~ubuntu-branches/ubuntu/trusty/ffc/trusty : revision 4

1

// This code conforms with the UFC specification version 1.0

2

// and was automatically generated by FFC version 0.7.0.

3

4

#ifndef __OPTIMIZATION_H

5

#define __OPTIMIZATION_H

6

7

#include <cmath>

8

#include <stdexcept>

9

#include <ufc.h>

10

11

/// This class defines the interface for a finite element.

12

13

class optimization_0_finite_element_0: public ufc::finite_element

14

{

15

public:

16

17

/// Constructor

18

optimization_0_finite_element_0() : ufc::finite_element()

19

{

20

// Do nothing

21

}

22

23

/// Destructor

24

virtual ~optimization_0_finite_element_0()

25

{

26

// Do nothing

27

}

28

29

/// Return a string identifying the finite element

30

virtual const char* signature() const

31

{

32

return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3)";

33

}

34

35

/// Return the cell shape

36

virtual ufc::shape cell_shape() const

37

{

38

return ufc::triangle;

39

}

40

41

/// Return the dimension of the finite element function space

42

virtual unsigned int space_dimension() const

43

{

44

return 10;

45

}

46

47

/// Return the rank of the value space

48

virtual unsigned int value_rank() const

49

{

50

return 0;

51

}

52

53

/// Return the dimension of the value space for axis i

54

virtual unsigned int value_dimension(unsigned int i) const

55

{

56

return 1;

57

}

58

59

/// Evaluate basis function i at given point in cell

60

virtual void evaluate_basis(unsigned int i,

61

double* values,

62

const double* coordinates,

63

const ufc::cell& c) const

64

{

65

// Extract vertex coordinates

66

const double * const * element_coordinates = c.coordinates;

67

68

// Compute Jacobian of affine map from reference cell

69

const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];

70

const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];

71

const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];

72

const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];

73

74

// Compute determinant of Jacobian

75

const double detJ = J_00*J_11 - J_01*J_10;

76

77

// Compute inverse of Jacobian

78

79

// Get coordinates and map to the reference (UFC) element

80

double x = (element_coordinates[0][1]*element_coordinates[2][0] -\

81

element_coordinates[0][0]*element_coordinates[2][1] +\

82

J_11*coordinates[0] - J_01*coordinates[1]) / detJ;

83

double y = (element_coordinates[1][1]*element_coordinates[0][0] -\

84

element_coordinates[1][0]*element_coordinates[0][1] -\

85

J_10*coordinates[0] + J_00*coordinates[1]) / detJ;

86

87

// Map coordinates to the reference square

88

if (std::abs(y - 1.0) < 1e-08)

89

x = -1.0;

90

else

91

x = 2.0 *x/(1.0 - y) - 1.0;

92

y = 2.0*y - 1.0;

93

94

// Reset values

95

*values = 0;

96

97

// Map degree of freedom to element degree of freedom

98

const unsigned int dof = i;

99

100

// Generate scalings

101

const double scalings_y_0 = 1;

102

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

103

const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);

104

const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);

105

106

// Compute psitilde_a

107

const double psitilde_a_0 = 1;

108

const double psitilde_a_1 = x;

109

const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;

110

const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;

111

112

// Compute psitilde_bs

113

const double psitilde_bs_0_0 = 1;

114

const double psitilde_bs_0_1 = 1.5*y + 0.5;

115

const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;

116

const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;

117

const double psitilde_bs_1_0 = 1;

118

const double psitilde_bs_1_1 = 2.5*y + 1.5;

119

const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;

120

const double psitilde_bs_2_0 = 1;

121

const double psitilde_bs_2_1 = 3.5*y + 2.5;

122

const double psitilde_bs_3_0 = 1;

123

124

// Compute basisvalues

125

const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;

126

const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;

127

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

128

const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;

129

const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;

130

const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;

131

const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;

132

const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;

133

const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;

134

const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;

135

136

// Table(s) of coefficients

137

static const double coefficients0[10][10] = \

138

{{0.0471404521, -0.0288675135, -0.0166666667, 0.0782460796, 0.0606091527, 0.0349927106, -0.0601337794, -0.0508223195, -0.0393667994, -0.0227284323},

139

{0.0471404521, 0.0288675135, -0.0166666667, 0.0782460796, -0.0606091527, 0.0349927106, 0.0601337794, -0.0508223195, 0.0393667994, -0.0227284323},

140

{0.0471404521, 0, 0.0333333333, 0, 0, 0.104978132, 0, 0, 0, 0.090913729},

141

{0.106066017, 0.259807621, -0.15, 0.117369119, 0.0606091527, -0.0787335989, 0, 0.101644639, -0.131222665, 0.090913729},

142

{0.106066017, 0, 0.3, 0, 0.151522882, 0.026244533, 0, 0, 0.131222665, -0.136370594},

143

{0.106066017, -0.259807621, -0.15, 0.117369119, -0.0606091527, -0.0787335989, 0, 0.101644639, 0.131222665, 0.090913729},

144

{0.106066017, 0, 0.3, 0, -0.151522882, 0.026244533, 0, 0, -0.131222665, -0.136370594},

145

{0.106066017, -0.259807621, -0.15, -0.0782460796, 0.090913729, 0.0962299542, 0.180401338, 0.0508223195, -0.0131222665, -0.0227284323},

146

{0.106066017, 0.259807621, -0.15, -0.0782460796, -0.090913729, 0.0962299542, -0.180401338, 0.0508223195, 0.0131222665, -0.0227284323},

147

{0.636396103, 0, 0, -0.234738239, 0, -0.26244533, 0, -0.203289278, 0, 0.090913729}};

148

149

// Extract relevant coefficients

150

const double coeff0_0 = coefficients0[dof][0];

151

const double coeff0_1 = coefficients0[dof][1];

152

const double coeff0_2 = coefficients0[dof][2];

153

const double coeff0_3 = coefficients0[dof][3];

154

const double coeff0_4 = coefficients0[dof][4];

155

const double coeff0_5 = coefficients0[dof][5];

156

const double coeff0_6 = coefficients0[dof][6];

157

const double coeff0_7 = coefficients0[dof][7];

158

const double coeff0_8 = coefficients0[dof][8];

159

const double coeff0_9 = coefficients0[dof][9];

160

161

// Compute value(s)

162

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9;

163

}

164

165

/// Evaluate all basis functions at given point in cell

166

virtual void evaluate_basis_all(double* values,

167

const double* coordinates,

168

const ufc::cell& c) const

169

{

170

throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");

171

}

172

173

/// Evaluate order n derivatives of basis function i at given point in cell

174

virtual void evaluate_basis_derivatives(unsigned int i,

175

unsigned int n,

176

double* values,

177

const double* coordinates,

178

const ufc::cell& c) const

179

{

180

// Extract vertex coordinates

181

const double * const * element_coordinates = c.coordinates;

182

183

// Compute Jacobian of affine map from reference cell

184

const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];

185

const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];

186

const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];

187

const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];

188

189

// Compute determinant of Jacobian

190

const double detJ = J_00*J_11 - J_01*J_10;

191

192

// Compute inverse of Jacobian

193

194

// Get coordinates and map to the reference (UFC) element

195

double x = (element_coordinates[0][1]*element_coordinates[2][0] -\

196

element_coordinates[0][0]*element_coordinates[2][1] +\

197

J_11*coordinates[0] - J_01*coordinates[1]) / detJ;

198

double y = (element_coordinates[1][1]*element_coordinates[0][0] -\

199

element_coordinates[1][0]*element_coordinates[0][1] -\

200

J_10*coordinates[0] + J_00*coordinates[1]) / detJ;

201

202

// Map coordinates to the reference square

203

if (std::abs(y - 1.0) < 1e-08)

204

x = -1.0;

205

else

206

x = 2.0 *x/(1.0 - y) - 1.0;

207

y = 2.0*y - 1.0;

208

209

// Compute number of derivatives

210

unsigned int num_derivatives = 1;

211

212

for (unsigned int j = 0; j < n; j++)

213

num_derivatives *= 2;

214

215

216

// Declare pointer to two dimensional array that holds combinations of derivatives and initialise

217

unsigned int **combinations = new unsigned int *[num_derivatives];

218

219

for (unsigned int j = 0; j < num_derivatives; j++)

220

{

221

combinations[j] = new unsigned int [n];

222

for (unsigned int k = 0; k < n; k++)

223

combinations[j][k] = 0;

224

}

225

226

// Generate combinations of derivatives

227

for (unsigned int row = 1; row < num_derivatives; row++)

228

{

229

for (unsigned int num = 0; num < row; num++)

230

{

231

for (unsigned int col = n-1; col+1 > 0; col--)

232

{

233

if (combinations[row][col] + 1 > 1)

234

combinations[row][col] = 0;

235

else

236

{

237

combinations[row][col] += 1;

238

break;

239

}

240

}

241

}

242

}

243

244

// Compute inverse of Jacobian

245

const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};

246

247

// Declare transformation matrix

248

// Declare pointer to two dimensional array and initialise

249

double **transform = new double *[num_derivatives];

250

251

for (unsigned int j = 0; j < num_derivatives; j++)

252

{

253

transform[j] = new double [num_derivatives];

254

for (unsigned int k = 0; k < num_derivatives; k++)

255

transform[j][k] = 1;

256

}

257

258

// Construct transformation matrix

259

for (unsigned int row = 0; row < num_derivatives; row++)

260

{

261

for (unsigned int col = 0; col < num_derivatives; col++)

262

{

263

for (unsigned int k = 0; k < n; k++)

264

transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];

265

}

266

}

267

268

// Reset values

269

for (unsigned int j = 0; j < 1*num_derivatives; j++)

270

values[j] = 0;

271

272

// Map degree of freedom to element degree of freedom

273

const unsigned int dof = i;

274

275

// Generate scalings

276

const double scalings_y_0 = 1;

277

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

278

const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);

279

const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);

280

281

// Compute psitilde_a

282

const double psitilde_a_0 = 1;

283

const double psitilde_a_1 = x;

284

const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;

285

const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;

286

287

// Compute psitilde_bs

288

const double psitilde_bs_0_0 = 1;

289

const double psitilde_bs_0_1 = 1.5*y + 0.5;

290

const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;

291

const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;

292

const double psitilde_bs_1_0 = 1;

293

const double psitilde_bs_1_1 = 2.5*y + 1.5;

294

const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;

295

const double psitilde_bs_2_0 = 1;

296

const double psitilde_bs_2_1 = 3.5*y + 2.5;

297

const double psitilde_bs_3_0 = 1;

298

299

// Compute basisvalues

300

const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;

301

const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;

302

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

303

const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;

304

const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;

305

const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;

306

const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;

307

const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;

308

const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;

309

const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;

310

311

// Table(s) of coefficients

312

static const double coefficients0[10][10] = \

313

{{0.0471404521, -0.0288675135, -0.0166666667, 0.0782460796, 0.0606091527, 0.0349927106, -0.0601337794, -0.0508223195, -0.0393667994, -0.0227284323},

314

{0.0471404521, 0.0288675135, -0.0166666667, 0.0782460796, -0.0606091527, 0.0349927106, 0.0601337794, -0.0508223195, 0.0393667994, -0.0227284323},

315

{0.0471404521, 0, 0.0333333333, 0, 0, 0.104978132, 0, 0, 0, 0.090913729},

316

{0.106066017, 0.259807621, -0.15, 0.117369119, 0.0606091527, -0.0787335989, 0, 0.101644639, -0.131222665, 0.090913729},

317

{0.106066017, 0, 0.3, 0, 0.151522882, 0.026244533, 0, 0, 0.131222665, -0.136370594},

318

{0.106066017, -0.259807621, -0.15, 0.117369119, -0.0606091527, -0.0787335989, 0, 0.101644639, 0.131222665, 0.090913729},

319

{0.106066017, 0, 0.3, 0, -0.151522882, 0.026244533, 0, 0, -0.131222665, -0.136370594},

320

{0.106066017, -0.259807621, -0.15, -0.0782460796, 0.090913729, 0.0962299542, 0.180401338, 0.0508223195, -0.0131222665, -0.0227284323},

321

{0.106066017, 0.259807621, -0.15, -0.0782460796, -0.090913729, 0.0962299542, -0.180401338, 0.0508223195, 0.0131222665, -0.0227284323},

322

{0.636396103, 0, 0, -0.234738239, 0, -0.26244533, 0, -0.203289278, 0, 0.090913729}};

323

324

// Interesting (new) part

325

// Tables of derivatives of the polynomial base (transpose)

326

static const double dmats0[10][10] = \

327

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

328

{4.89897949, 0, 0, 0, 0, 0, 0, 0, 0, 0},

329

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

330

{0, 9.48683298, 0, 0, 0, 0, 0, 0, 0, 0},

331

{4, 0, 7.07106781, 0, 0, 0, 0, 0, 0, 0},

332

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

333

{5.29150262, 0, -2.99332591, 13.662601, 0, 0.611010093, 0, 0, 0, 0},

334

{0, 4.38178046, 0, 0, 12.5219807, 0, 0, 0, 0, 0},

335

{3.46410162, 0, 7.83836718, 0, 0, 8.4, 0, 0, 0, 0},

336

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};

337

338

static const double dmats1[10][10] = \

339

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

340

{2.44948974, 0, 0, 0, 0, 0, 0, 0, 0, 0},

341

{4.24264069, 0, 0, 0, 0, 0, 0, 0, 0, 0},

342

{2.5819889, 4.74341649, -0.912870929, 0, 0, 0, 0, 0, 0, 0},

343

{2, 6.12372436, 3.53553391, 0, 0, 0, 0, 0, 0, 0},

344

{-2.30940108, 0, 8.16496581, 0, 0, 0, 0, 0, 0, 0},

345

{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.305505046, 0, 0, 0, 0},

346

{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0, 0, 0, 0},

347

{1.73205081, -5.09116882, 3.91918359, 0, 9.69948452, 4.2, 0, 0, 0, 0},

348

{5, 0, -2.82842712, 0, 0, 12.1243557, 0, 0, 0, 0}};

349

350

// Compute reference derivatives

351

// Declare pointer to array of derivatives on FIAT element

352

double *derivatives = new double [num_derivatives];

353

354

// Declare coefficients

355

double coeff0_0 = 0;

356

double coeff0_1 = 0;

357

double coeff0_2 = 0;

358

double coeff0_3 = 0;

359

double coeff0_4 = 0;

360

double coeff0_5 = 0;

361

double coeff0_6 = 0;

362

double coeff0_7 = 0;

363

double coeff0_8 = 0;

364

double coeff0_9 = 0;

365

366

// Declare new coefficients

367

double new_coeff0_0 = 0;

368

double new_coeff0_1 = 0;

369

double new_coeff0_2 = 0;

370

double new_coeff0_3 = 0;

371

double new_coeff0_4 = 0;

372

double new_coeff0_5 = 0;

373

double new_coeff0_6 = 0;

374

double new_coeff0_7 = 0;

375

double new_coeff0_8 = 0;

376

double new_coeff0_9 = 0;

377

378

// Loop possible derivatives

379

for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)

380

{

381

// Get values from coefficients array

382

new_coeff0_0 = coefficients0[dof][0];

383

new_coeff0_1 = coefficients0[dof][1];

384

new_coeff0_2 = coefficients0[dof][2];

385

new_coeff0_3 = coefficients0[dof][3];

386

new_coeff0_4 = coefficients0[dof][4];

387

new_coeff0_5 = coefficients0[dof][5];

388

new_coeff0_6 = coefficients0[dof][6];

389

new_coeff0_7 = coefficients0[dof][7];

390

new_coeff0_8 = coefficients0[dof][8];

391

new_coeff0_9 = coefficients0[dof][9];

392

393

// Loop derivative order

394

for (unsigned int j = 0; j < n; j++)

395

{

396

// Update old coefficients

397

coeff0_0 = new_coeff0_0;

398

coeff0_1 = new_coeff0_1;

399

coeff0_2 = new_coeff0_2;

400

coeff0_3 = new_coeff0_3;

401

coeff0_4 = new_coeff0_4;

402

coeff0_5 = new_coeff0_5;

403

coeff0_6 = new_coeff0_6;

404

coeff0_7 = new_coeff0_7;

405

coeff0_8 = new_coeff0_8;

406

coeff0_9 = new_coeff0_9;

407

408

if(combinations[deriv_num][j] == 0)

409

{

410

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0] + coeff0_6*dmats0[6][0] + coeff0_7*dmats0[7][0] + coeff0_8*dmats0[8][0] + coeff0_9*dmats0[9][0];

411

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1] + coeff0_6*dmats0[6][1] + coeff0_7*dmats0[7][1] + coeff0_8*dmats0[8][1] + coeff0_9*dmats0[9][1];

412

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2] + coeff0_6*dmats0[6][2] + coeff0_7*dmats0[7][2] + coeff0_8*dmats0[8][2] + coeff0_9*dmats0[9][2];

413

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3] + coeff0_6*dmats0[6][3] + coeff0_7*dmats0[7][3] + coeff0_8*dmats0[8][3] + coeff0_9*dmats0[9][3];

414

new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4] + coeff0_6*dmats0[6][4] + coeff0_7*dmats0[7][4] + coeff0_8*dmats0[8][4] + coeff0_9*dmats0[9][4];

415

new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5] + coeff0_6*dmats0[6][5] + coeff0_7*dmats0[7][5] + coeff0_8*dmats0[8][5] + coeff0_9*dmats0[9][5];

416

new_coeff0_6 = coeff0_0*dmats0[0][6] + coeff0_1*dmats0[1][6] + coeff0_2*dmats0[2][6] + coeff0_3*dmats0[3][6] + coeff0_4*dmats0[4][6] + coeff0_5*dmats0[5][6] + coeff0_6*dmats0[6][6] + coeff0_7*dmats0[7][6] + coeff0_8*dmats0[8][6] + coeff0_9*dmats0[9][6];

417

new_coeff0_7 = coeff0_0*dmats0[0][7] + coeff0_1*dmats0[1][7] + coeff0_2*dmats0[2][7] + coeff0_3*dmats0[3][7] + coeff0_4*dmats0[4][7] + coeff0_5*dmats0[5][7] + coeff0_6*dmats0[6][7] + coeff0_7*dmats0[7][7] + coeff0_8*dmats0[8][7] + coeff0_9*dmats0[9][7];

418

new_coeff0_8 = coeff0_0*dmats0[0][8] + coeff0_1*dmats0[1][8] + coeff0_2*dmats0[2][8] + coeff0_3*dmats0[3][8] + coeff0_4*dmats0[4][8] + coeff0_5*dmats0[5][8] + coeff0_6*dmats0[6][8] + coeff0_7*dmats0[7][8] + coeff0_8*dmats0[8][8] + coeff0_9*dmats0[9][8];

419

new_coeff0_9 = coeff0_0*dmats0[0][9] + coeff0_1*dmats0[1][9] + coeff0_2*dmats0[2][9] + coeff0_3*dmats0[3][9] + coeff0_4*dmats0[4][9] + coeff0_5*dmats0[5][9] + coeff0_6*dmats0[6][9] + coeff0_7*dmats0[7][9] + coeff0_8*dmats0[8][9] + coeff0_9*dmats0[9][9];

420

}

421

if(combinations[deriv_num][j] == 1)

422

{

423

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0] + coeff0_6*dmats1[6][0] + coeff0_7*dmats1[7][0] + coeff0_8*dmats1[8][0] + coeff0_9*dmats1[9][0];

424

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1] + coeff0_6*dmats1[6][1] + coeff0_7*dmats1[7][1] + coeff0_8*dmats1[8][1] + coeff0_9*dmats1[9][1];

425

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2] + coeff0_6*dmats1[6][2] + coeff0_7*dmats1[7][2] + coeff0_8*dmats1[8][2] + coeff0_9*dmats1[9][2];

426

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3] + coeff0_6*dmats1[6][3] + coeff0_7*dmats1[7][3] + coeff0_8*dmats1[8][3] + coeff0_9*dmats1[9][3];

427

new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4] + coeff0_6*dmats1[6][4] + coeff0_7*dmats1[7][4] + coeff0_8*dmats1[8][4] + coeff0_9*dmats1[9][4];

428

new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5] + coeff0_6*dmats1[6][5] + coeff0_7*dmats1[7][5] + coeff0_8*dmats1[8][5] + coeff0_9*dmats1[9][5];

429

new_coeff0_6 = coeff0_0*dmats1[0][6] + coeff0_1*dmats1[1][6] + coeff0_2*dmats1[2][6] + coeff0_3*dmats1[3][6] + coeff0_4*dmats1[4][6] + coeff0_5*dmats1[5][6] + coeff0_6*dmats1[6][6] + coeff0_7*dmats1[7][6] + coeff0_8*dmats1[8][6] + coeff0_9*dmats1[9][6];

430

new_coeff0_7 = coeff0_0*dmats1[0][7] + coeff0_1*dmats1[1][7] + coeff0_2*dmats1[2][7] + coeff0_3*dmats1[3][7] + coeff0_4*dmats1[4][7] + coeff0_5*dmats1[5][7] + coeff0_6*dmats1[6][7] + coeff0_7*dmats1[7][7] + coeff0_8*dmats1[8][7] + coeff0_9*dmats1[9][7];

431

new_coeff0_8 = coeff0_0*dmats1[0][8] + coeff0_1*dmats1[1][8] + coeff0_2*dmats1[2][8] + coeff0_3*dmats1[3][8] + coeff0_4*dmats1[4][8] + coeff0_5*dmats1[5][8] + coeff0_6*dmats1[6][8] + coeff0_7*dmats1[7][8] + coeff0_8*dmats1[8][8] + coeff0_9*dmats1[9][8];

432

new_coeff0_9 = coeff0_0*dmats1[0][9] + coeff0_1*dmats1[1][9] + coeff0_2*dmats1[2][9] + coeff0_3*dmats1[3][9] + coeff0_4*dmats1[4][9] + coeff0_5*dmats1[5][9] + coeff0_6*dmats1[6][9] + coeff0_7*dmats1[7][9] + coeff0_8*dmats1[8][9] + coeff0_9*dmats1[9][9];

433

}

434

435

}

436

// Compute derivatives on reference element as dot product of coefficients and basisvalues

437

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5 + new_coeff0_6*basisvalue6 + new_coeff0_7*basisvalue7 + new_coeff0_8*basisvalue8 + new_coeff0_9*basisvalue9;

438

}

439

440

// Transform derivatives back to physical element

441

for (unsigned int row = 0; row < num_derivatives; row++)

442

{

443

for (unsigned int col = 0; col < num_derivatives; col++)

444

{

445

values[row] += transform[row][col]*derivatives[col];

446

}

447

}

448

// Delete pointer to array of derivatives on FIAT element

449

delete [] derivatives;

450

451

// Delete pointer to array of combinations of derivatives and transform

452

for (unsigned int row = 0; row < num_derivatives; row++)

453

{

454

delete [] combinations[row];

455

delete [] transform[row];

456

}

457

458

delete [] combinations;

459

delete [] transform;

460

}

461

462

/// Evaluate order n derivatives of all basis functions at given point in cell

463

virtual void evaluate_basis_derivatives_all(unsigned int n,

464

double* values,

465

const double* coordinates,

466

const ufc::cell& c) const

467

{

468

throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");

469

}

470

471

/// Evaluate linear functional for dof i on the function f

472

virtual double evaluate_dof(unsigned int i,

473

const ufc::function& f,

474

const ufc::cell& c) const

475

{

476

// The reference points, direction and weights:

477

static const double X[10][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.666666667, 0.333333333}}, {{0.333333333, 0.666666667}}, {{0, 0.333333333}}, {{0, 0.666666667}}, {{0.333333333, 0}}, {{0.666666667, 0}}, {{0.333333333, 0.333333333}}};

478

static const double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};

479

static const double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};

480

481

const double * const * x = c.coordinates;

482

double result = 0.0;

483

// Iterate over the points:

484

// Evaluate basis functions for affine mapping

485

const double w0 = 1.0 - X[i][0][0] - X[i][0][1];

486

const double w1 = X[i][0][0];

487

const double w2 = X[i][0][1];

488

489

// Compute affine mapping y = F(X)

490

double y[2];

491

y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];

492

y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];

493

494

// Evaluate function at physical points

495

double values[1];

496

f.evaluate(values, y, c);

497

498

// Map function values using appropriate mapping

499

// Affine map: Do nothing

500

501

// Note that we do not map the weights (yet).

502

503

// Take directional components

504

for(int k = 0; k < 1; k++)

505

result += values[k]*D[i][0][k];

506

// Multiply by weights

507

result *= W[i][0];

508

509

return result;

510

}

511

512

/// Evaluate linear functionals for all dofs on the function f

513

virtual void evaluate_dofs(double* values,

514

const ufc::function& f,

515

const ufc::cell& c) const

516

{

517

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

518

}

519

520

/// Interpolate vertex values from dof values

521

virtual void interpolate_vertex_values(double* vertex_values,

522

const double* dof_values,

523

const ufc::cell& c) const

524

{

525

// Evaluate at vertices and use affine mapping

526

vertex_values[0] = dof_values[0];

527

vertex_values[1] = dof_values[1];

528

vertex_values[2] = dof_values[2];

529

}

530

531

/// Return the number of sub elements (for a mixed element)

532

virtual unsigned int num_sub_elements() const

533

{

534

return 1;

535

}

536

537

/// Create a new finite element for sub element i (for a mixed element)

538

virtual ufc::finite_element* create_sub_element(unsigned int i) const

539

{

540

return new optimization_0_finite_element_0();

541

}

542

543

};

544

545

/// This class defines the interface for a finite element.

546

547

class optimization_0_finite_element_1: public ufc::finite_element

548

{

549

public:

550

551

/// Constructor

552

optimization_0_finite_element_1() : ufc::finite_element()

553

{

554

// Do nothing

555

}

556

557

/// Destructor

558

virtual ~optimization_0_finite_element_1()

559

{

560

// Do nothing

561

}

562

563

/// Return a string identifying the finite element

564

virtual const char* signature() const

565

{

566

return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3)";

567

}

568

569

/// Return the cell shape

570

virtual ufc::shape cell_shape() const

571

{

572

return ufc::triangle;

573

}

574

575

/// Return the dimension of the finite element function space

576

virtual unsigned int space_dimension() const

577

{

578

return 10;

579

}

580

581

/// Return the rank of the value space

582

virtual unsigned int value_rank() const

583

{

584

return 0;

585

}

586

587

/// Return the dimension of the value space for axis i

588

virtual unsigned int value_dimension(unsigned int i) const

589

{

590

return 1;

591

}

592

593

/// Evaluate basis function i at given point in cell

594

virtual void evaluate_basis(unsigned int i,

595

double* values,

596

const double* coordinates,

597

const ufc::cell& c) const

598

{

599

// Extract vertex coordinates

600

const double * const * element_coordinates = c.coordinates;

601

602

// Compute Jacobian of affine map from reference cell

603

const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];

604

const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];

605

const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];

606

const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];

607

608

// Compute determinant of Jacobian

609

const double detJ = J_00*J_11 - J_01*J_10;

610

611

// Compute inverse of Jacobian

612

613

// Get coordinates and map to the reference (UFC) element

614

double x = (element_coordinates[0][1]*element_coordinates[2][0] -\

615

element_coordinates[0][0]*element_coordinates[2][1] +\

616

J_11*coordinates[0] - J_01*coordinates[1]) / detJ;

617

double y = (element_coordinates[1][1]*element_coordinates[0][0] -\

618

element_coordinates[1][0]*element_coordinates[0][1] -\

619

J_10*coordinates[0] + J_00*coordinates[1]) / detJ;

620

621

// Map coordinates to the reference square

622

if (std::abs(y - 1.0) < 1e-08)

623

x = -1.0;

624

else

625

x = 2.0 *x/(1.0 - y) - 1.0;

626

y = 2.0*y - 1.0;

627

628

// Reset values

629

*values = 0;

630

631

// Map degree of freedom to element degree of freedom

632

const unsigned int dof = i;

633

634

// Generate scalings

635

const double scalings_y_0 = 1;

636

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

637

const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);

638

const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);

639

640

// Compute psitilde_a

641

const double psitilde_a_0 = 1;

642

const double psitilde_a_1 = x;

643

const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;

644

const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;

645

646

// Compute psitilde_bs

647

const double psitilde_bs_0_0 = 1;

648

const double psitilde_bs_0_1 = 1.5*y + 0.5;

649

const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;

650

const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;

651

const double psitilde_bs_1_0 = 1;

652

const double psitilde_bs_1_1 = 2.5*y + 1.5;

653

const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;

654

const double psitilde_bs_2_0 = 1;

655

const double psitilde_bs_2_1 = 3.5*y + 2.5;

656

const double psitilde_bs_3_0 = 1;

657

658

// Compute basisvalues

659

const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;

660

const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;

661

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

662

const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;

663

const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;

664

const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;

665

const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;

666

const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;

667

const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;

668

const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;

669

670

// Table(s) of coefficients

671

static const double coefficients0[10][10] = \

672

{{0.0471404521, -0.0288675135, -0.0166666667, 0.0782460796, 0.0606091527, 0.0349927106, -0.0601337794, -0.0508223195, -0.0393667994, -0.0227284323},

673

{0.0471404521, 0.0288675135, -0.0166666667, 0.0782460796, -0.0606091527, 0.0349927106, 0.0601337794, -0.0508223195, 0.0393667994, -0.0227284323},

674

{0.0471404521, 0, 0.0333333333, 0, 0, 0.104978132, 0, 0, 0, 0.090913729},

675

{0.106066017, 0.259807621, -0.15, 0.117369119, 0.0606091527, -0.0787335989, 0, 0.101644639, -0.131222665, 0.090913729},

676

{0.106066017, 0, 0.3, 0, 0.151522882, 0.026244533, 0, 0, 0.131222665, -0.136370594},

677

{0.106066017, -0.259807621, -0.15, 0.117369119, -0.0606091527, -0.0787335989, 0, 0.101644639, 0.131222665, 0.090913729},

678

{0.106066017, 0, 0.3, 0, -0.151522882, 0.026244533, 0, 0, -0.131222665, -0.136370594},

679

{0.106066017, -0.259807621, -0.15, -0.0782460796, 0.090913729, 0.0962299542, 0.180401338, 0.0508223195, -0.0131222665, -0.0227284323},

680

{0.106066017, 0.259807621, -0.15, -0.0782460796, -0.090913729, 0.0962299542, -0.180401338, 0.0508223195, 0.0131222665, -0.0227284323},

681

{0.636396103, 0, 0, -0.234738239, 0, -0.26244533, 0, -0.203289278, 0, 0.090913729}};

682

683

// Extract relevant coefficients

684

const double coeff0_0 = coefficients0[dof][0];

685

const double coeff0_1 = coefficients0[dof][1];

686

const double coeff0_2 = coefficients0[dof][2];

687

const double coeff0_3 = coefficients0[dof][3];

688

const double coeff0_4 = coefficients0[dof][4];

689

const double coeff0_5 = coefficients0[dof][5];

690

const double coeff0_6 = coefficients0[dof][6];

691

const double coeff0_7 = coefficients0[dof][7];

692

const double coeff0_8 = coefficients0[dof][8];

693

const double coeff0_9 = coefficients0[dof][9];

694

695

// Compute value(s)

696

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9;

697

}

698

699

/// Evaluate all basis functions at given point in cell

700

virtual void evaluate_basis_all(double* values,

701

const double* coordinates,

702

const ufc::cell& c) const

703

{

704

throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");

705

}

706

707

/// Evaluate order n derivatives of basis function i at given point in cell

708

virtual void evaluate_basis_derivatives(unsigned int i,

709

unsigned int n,

710

double* values,

711

const double* coordinates,

712

const ufc::cell& c) const

713

{

714

// Extract vertex coordinates

715

const double * const * element_coordinates = c.coordinates;

716

717

// Compute Jacobian of affine map from reference cell

718

const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];

719

const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];

720

const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];

721

const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];

722

723

// Compute determinant of Jacobian

724

const double detJ = J_00*J_11 - J_01*J_10;

725

726

// Compute inverse of Jacobian

727

728

// Get coordinates and map to the reference (UFC) element

729

double x = (element_coordinates[0][1]*element_coordinates[2][0] -\

730

element_coordinates[0][0]*element_coordinates[2][1] +\

731

J_11*coordinates[0] - J_01*coordinates[1]) / detJ;

732

double y = (element_coordinates[1][1]*element_coordinates[0][0] -\

733

element_coordinates[1][0]*element_coordinates[0][1] -\

734

J_10*coordinates[0] + J_00*coordinates[1]) / detJ;

735

736

// Map coordinates to the reference square

737

if (std::abs(y - 1.0) < 1e-08)

738

x = -1.0;

739

else

740

x = 2.0 *x/(1.0 - y) - 1.0;

741

y = 2.0*y - 1.0;

742

743

// Compute number of derivatives

744

unsigned int num_derivatives = 1;

745

746

for (unsigned int j = 0; j < n; j++)

747

num_derivatives *= 2;

748

749

750

// Declare pointer to two dimensional array that holds combinations of derivatives and initialise

751

unsigned int **combinations = new unsigned int *[num_derivatives];

752

753

for (unsigned int j = 0; j < num_derivatives; j++)

754

{

755

combinations[j] = new unsigned int [n];

756

for (unsigned int k = 0; k < n; k++)

757

combinations[j][k] = 0;

758

}

759

760

// Generate combinations of derivatives

761

for (unsigned int row = 1; row < num_derivatives; row++)

762

{

763

for (unsigned int num = 0; num < row; num++)

764

{

765

for (unsigned int col = n-1; col+1 > 0; col--)

766

{

767

if (combinations[row][col] + 1 > 1)

768

combinations[row][col] = 0;

769

else

770

{

771

combinations[row][col] += 1;

772

break;

773

}

774

}

775

}

776

}

777

778

// Compute inverse of Jacobian

779

const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};

780

781

// Declare transformation matrix

782

// Declare pointer to two dimensional array and initialise

783

double **transform = new double *[num_derivatives];

784

785

for (unsigned int j = 0; j < num_derivatives; j++)

786

{

787

transform[j] = new double [num_derivatives];

788

for (unsigned int k = 0; k < num_derivatives; k++)

789

transform[j][k] = 1;

790

}

791

792

// Construct transformation matrix

793

for (unsigned int row = 0; row < num_derivatives; row++)

794

{

795

for (unsigned int col = 0; col < num_derivatives; col++)

796

{

797

for (unsigned int k = 0; k < n; k++)

798

transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];

799

}

800

}

801

802

// Reset values

803

for (unsigned int j = 0; j < 1*num_derivatives; j++)

804

values[j] = 0;

805

806

// Map degree of freedom to element degree of freedom

807

const unsigned int dof = i;

808

809

// Generate scalings

810

const double scalings_y_0 = 1;

811

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

812

const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);

813

const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);

814

815

// Compute psitilde_a

816

const double psitilde_a_0 = 1;

817

const double psitilde_a_1 = x;

818

const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;

819

const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;

820

821

// Compute psitilde_bs

822

const double psitilde_bs_0_0 = 1;

823

const double psitilde_bs_0_1 = 1.5*y + 0.5;

824

const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;

825

const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;

826

const double psitilde_bs_1_0 = 1;

827

const double psitilde_bs_1_1 = 2.5*y + 1.5;

828

const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;

829

const double psitilde_bs_2_0 = 1;

830

const double psitilde_bs_2_1 = 3.5*y + 2.5;

831

const double psitilde_bs_3_0 = 1;

832

833

// Compute basisvalues

834

const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;

835

const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;

836

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

837

const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;

838

const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;

839

const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;

840

const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;

841

const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;

842

const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;

843

const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;

844

845

// Table(s) of coefficients

846

static const double coefficients0[10][10] = \

847

{{0.0471404521, -0.0288675135, -0.0166666667, 0.0782460796, 0.0606091527, 0.0349927106, -0.0601337794, -0.0508223195, -0.0393667994, -0.0227284323},

848

{0.0471404521, 0.0288675135, -0.0166666667, 0.0782460796, -0.0606091527, 0.0349927106, 0.0601337794, -0.0508223195, 0.0393667994, -0.0227284323},

849

{0.0471404521, 0, 0.0333333333, 0, 0, 0.104978132, 0, 0, 0, 0.090913729},

850

{0.106066017, 0.259807621, -0.15, 0.117369119, 0.0606091527, -0.0787335989, 0, 0.101644639, -0.131222665, 0.090913729},

851

{0.106066017, 0, 0.3, 0, 0.151522882, 0.026244533, 0, 0, 0.131222665, -0.136370594},

852

{0.106066017, -0.259807621, -0.15, 0.117369119, -0.0606091527, -0.0787335989, 0, 0.101644639, 0.131222665, 0.090913729},

853

{0.106066017, 0, 0.3, 0, -0.151522882, 0.026244533, 0, 0, -0.131222665, -0.136370594},

854

{0.106066017, -0.259807621, -0.15, -0.0782460796, 0.090913729, 0.0962299542, 0.180401338, 0.0508223195, -0.0131222665, -0.0227284323},

855

{0.106066017, 0.259807621, -0.15, -0.0782460796, -0.090913729, 0.0962299542, -0.180401338, 0.0508223195, 0.0131222665, -0.0227284323},

856

{0.636396103, 0, 0, -0.234738239, 0, -0.26244533, 0, -0.203289278, 0, 0.090913729}};

857

858

// Interesting (new) part

859

// Tables of derivatives of the polynomial base (transpose)

860

static const double dmats0[10][10] = \

861

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

862

{4.89897949, 0, 0, 0, 0, 0, 0, 0, 0, 0},

863

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

864

{0, 9.48683298, 0, 0, 0, 0, 0, 0, 0, 0},

865

{4, 0, 7.07106781, 0, 0, 0, 0, 0, 0, 0},

866

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

867

{5.29150262, 0, -2.99332591, 13.662601, 0, 0.611010093, 0, 0, 0, 0},

868

{0, 4.38178046, 0, 0, 12.5219807, 0, 0, 0, 0, 0},

869

{3.46410162, 0, 7.83836718, 0, 0, 8.4, 0, 0, 0, 0},

870

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};

871

872

static const double dmats1[10][10] = \

873

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

874

{2.44948974, 0, 0, 0, 0, 0, 0, 0, 0, 0},

875

{4.24264069, 0, 0, 0, 0, 0, 0, 0, 0, 0},

876

{2.5819889, 4.74341649, -0.912870929, 0, 0, 0, 0, 0, 0, 0},

877

{2, 6.12372436, 3.53553391, 0, 0, 0, 0, 0, 0, 0},

878

{-2.30940108, 0, 8.16496581, 0, 0, 0, 0, 0, 0, 0},

879

{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.305505046, 0, 0, 0, 0},

880

{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0, 0, 0, 0},

881

{1.73205081, -5.09116882, 3.91918359, 0, 9.69948452, 4.2, 0, 0, 0, 0},

882

{5, 0, -2.82842712, 0, 0, 12.1243557, 0, 0, 0, 0}};

883

884

// Compute reference derivatives

885

// Declare pointer to array of derivatives on FIAT element

886

double *derivatives = new double [num_derivatives];

887

888

// Declare coefficients

889

double coeff0_0 = 0;

890

double coeff0_1 = 0;

891

double coeff0_2 = 0;

892

double coeff0_3 = 0;

893

double coeff0_4 = 0;

894

double coeff0_5 = 0;

895

double coeff0_6 = 0;

896

double coeff0_7 = 0;

897

double coeff0_8 = 0;

898

double coeff0_9 = 0;

899

900

// Declare new coefficients

901

double new_coeff0_0 = 0;

902

double new_coeff0_1 = 0;

903

double new_coeff0_2 = 0;

904

double new_coeff0_3 = 0;

905

double new_coeff0_4 = 0;

906

double new_coeff0_5 = 0;

907

double new_coeff0_6 = 0;

908

double new_coeff0_7 = 0;

909

double new_coeff0_8 = 0;

910

double new_coeff0_9 = 0;

911

912

// Loop possible derivatives

913

for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)

914

{

915

// Get values from coefficients array

916

new_coeff0_0 = coefficients0[dof][0];

917

new_coeff0_1 = coefficients0[dof][1];

918

new_coeff0_2 = coefficients0[dof][2];

919

new_coeff0_3 = coefficients0[dof][3];

920

new_coeff0_4 = coefficients0[dof][4];

921

new_coeff0_5 = coefficients0[dof][5];

922

new_coeff0_6 = coefficients0[dof][6];

923

new_coeff0_7 = coefficients0[dof][7];

924

new_coeff0_8 = coefficients0[dof][8];

925

new_coeff0_9 = coefficients0[dof][9];

926

927

// Loop derivative order

928

for (unsigned int j = 0; j < n; j++)

929

{

930

// Update old coefficients

931

coeff0_0 = new_coeff0_0;

932

coeff0_1 = new_coeff0_1;

933

coeff0_2 = new_coeff0_2;

934

coeff0_3 = new_coeff0_3;

935

coeff0_4 = new_coeff0_4;

936

coeff0_5 = new_coeff0_5;

937

coeff0_6 = new_coeff0_6;

938

coeff0_7 = new_coeff0_7;

939

coeff0_8 = new_coeff0_8;

940

coeff0_9 = new_coeff0_9;

941

942

if(combinations[deriv_num][j] == 0)

943

{

944

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0] + coeff0_6*dmats0[6][0] + coeff0_7*dmats0[7][0] + coeff0_8*dmats0[8][0] + coeff0_9*dmats0[9][0];

945

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1] + coeff0_6*dmats0[6][1] + coeff0_7*dmats0[7][1] + coeff0_8*dmats0[8][1] + coeff0_9*dmats0[9][1];

946

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2] + coeff0_6*dmats0[6][2] + coeff0_7*dmats0[7][2] + coeff0_8*dmats0[8][2] + coeff0_9*dmats0[9][2];

947

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3] + coeff0_6*dmats0[6][3] + coeff0_7*dmats0[7][3] + coeff0_8*dmats0[8][3] + coeff0_9*dmats0[9][3];

948

new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4] + coeff0_6*dmats0[6][4] + coeff0_7*dmats0[7][4] + coeff0_8*dmats0[8][4] + coeff0_9*dmats0[9][4];

949

new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5] + coeff0_6*dmats0[6][5] + coeff0_7*dmats0[7][5] + coeff0_8*dmats0[8][5] + coeff0_9*dmats0[9][5];

950

new_coeff0_6 = coeff0_0*dmats0[0][6] + coeff0_1*dmats0[1][6] + coeff0_2*dmats0[2][6] + coeff0_3*dmats0[3][6] + coeff0_4*dmats0[4][6] + coeff0_5*dmats0[5][6] + coeff0_6*dmats0[6][6] + coeff0_7*dmats0[7][6] + coeff0_8*dmats0[8][6] + coeff0_9*dmats0[9][6];

951

new_coeff0_7 = coeff0_0*dmats0[0][7] + coeff0_1*dmats0[1][7] + coeff0_2*dmats0[2][7] + coeff0_3*dmats0[3][7] + coeff0_4*dmats0[4][7] + coeff0_5*dmats0[5][7] + coeff0_6*dmats0[6][7] + coeff0_7*dmats0[7][7] + coeff0_8*dmats0[8][7] + coeff0_9*dmats0[9][7];

952

new_coeff0_8 = coeff0_0*dmats0[0][8] + coeff0_1*dmats0[1][8] + coeff0_2*dmats0[2][8] + coeff0_3*dmats0[3][8] + coeff0_4*dmats0[4][8] + coeff0_5*dmats0[5][8] + coeff0_6*dmats0[6][8] + coeff0_7*dmats0[7][8] + coeff0_8*dmats0[8][8] + coeff0_9*dmats0[9][8];

953

new_coeff0_9 = coeff0_0*dmats0[0][9] + coeff0_1*dmats0[1][9] + coeff0_2*dmats0[2][9] + coeff0_3*dmats0[3][9] + coeff0_4*dmats0[4][9] + coeff0_5*dmats0[5][9] + coeff0_6*dmats0[6][9] + coeff0_7*dmats0[7][9] + coeff0_8*dmats0[8][9] + coeff0_9*dmats0[9][9];

954

}

955

if(combinations[deriv_num][j] == 1)

956

{

957

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0] + coeff0_6*dmats1[6][0] + coeff0_7*dmats1[7][0] + coeff0_8*dmats1[8][0] + coeff0_9*dmats1[9][0];

958

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1] + coeff0_6*dmats1[6][1] + coeff0_7*dmats1[7][1] + coeff0_8*dmats1[8][1] + coeff0_9*dmats1[9][1];

959

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2] + coeff0_6*dmats1[6][2] + coeff0_7*dmats1[7][2] + coeff0_8*dmats1[8][2] + coeff0_9*dmats1[9][2];

960

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3] + coeff0_6*dmats1[6][3] + coeff0_7*dmats1[7][3] + coeff0_8*dmats1[8][3] + coeff0_9*dmats1[9][3];

961

new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4] + coeff0_6*dmats1[6][4] + coeff0_7*dmats1[7][4] + coeff0_8*dmats1[8][4] + coeff0_9*dmats1[9][4];

962

new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5] + coeff0_6*dmats1[6][5] + coeff0_7*dmats1[7][5] + coeff0_8*dmats1[8][5] + coeff0_9*dmats1[9][5];

963

new_coeff0_6 = coeff0_0*dmats1[0][6] + coeff0_1*dmats1[1][6] + coeff0_2*dmats1[2][6] + coeff0_3*dmats1[3][6] + coeff0_4*dmats1[4][6] + coeff0_5*dmats1[5][6] + coeff0_6*dmats1[6][6] + coeff0_7*dmats1[7][6] + coeff0_8*dmats1[8][6] + coeff0_9*dmats1[9][6];

964

new_coeff0_7 = coeff0_0*dmats1[0][7] + coeff0_1*dmats1[1][7] + coeff0_2*dmats1[2][7] + coeff0_3*dmats1[3][7] + coeff0_4*dmats1[4][7] + coeff0_5*dmats1[5][7] + coeff0_6*dmats1[6][7] + coeff0_7*dmats1[7][7] + coeff0_8*dmats1[8][7] + coeff0_9*dmats1[9][7];

965

new_coeff0_8 = coeff0_0*dmats1[0][8] + coeff0_1*dmats1[1][8] + coeff0_2*dmats1[2][8] + coeff0_3*dmats1[3][8] + coeff0_4*dmats1[4][8] + coeff0_5*dmats1[5][8] + coeff0_6*dmats1[6][8] + coeff0_7*dmats1[7][8] + coeff0_8*dmats1[8][8] + coeff0_9*dmats1[9][8];

966

new_coeff0_9 = coeff0_0*dmats1[0][9] + coeff0_1*dmats1[1][9] + coeff0_2*dmats1[2][9] + coeff0_3*dmats1[3][9] + coeff0_4*dmats1[4][9] + coeff0_5*dmats1[5][9] + coeff0_6*dmats1[6][9] + coeff0_7*dmats1[7][9] + coeff0_8*dmats1[8][9] + coeff0_9*dmats1[9][9];

967

}

968

969

}

970

// Compute derivatives on reference element as dot product of coefficients and basisvalues

971

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5 + new_coeff0_6*basisvalue6 + new_coeff0_7*basisvalue7 + new_coeff0_8*basisvalue8 + new_coeff0_9*basisvalue9;

972

}

973

974

// Transform derivatives back to physical element

975

for (unsigned int row = 0; row < num_derivatives; row++)

976

{

977

for (unsigned int col = 0; col < num_derivatives; col++)

978

{

979

values[row] += transform[row][col]*derivatives[col];

980

}

981

}

982

// Delete pointer to array of derivatives on FIAT element

983

delete [] derivatives;

984

985

// Delete pointer to array of combinations of derivatives and transform

986

for (unsigned int row = 0; row < num_derivatives; row++)

987

{

988

delete [] combinations[row];

989

delete [] transform[row];

990

}

991

992

delete [] combinations;

993

delete [] transform;

994

}

995

996

/// Evaluate order n derivatives of all basis functions at given point in cell

997

virtual void evaluate_basis_derivatives_all(unsigned int n,

998

double* values,

999

const double* coordinates,

1000

const ufc::cell& c) const

1001

{

1002

throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");

1003

}

1004

1005

/// Evaluate linear functional for dof i on the function f

1006

virtual double evaluate_dof(unsigned int i,

1007

const ufc::function& f,

1008

const ufc::cell& c) const

1009

{

1010

// The reference points, direction and weights:

1011

static const double X[10][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.666666667, 0.333333333}}, {{0.333333333, 0.666666667}}, {{0, 0.333333333}}, {{0, 0.666666667}}, {{0.333333333, 0}}, {{0.666666667, 0}}, {{0.333333333, 0.333333333}}};

1012

static const double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};

1013

static const double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};

1014

1015

const double * const * x = c.coordinates;

1016

double result = 0.0;

1017

// Iterate over the points:

1018

// Evaluate basis functions for affine mapping

1019

const double w0 = 1.0 - X[i][0][0] - X[i][0][1];

1020

const double w1 = X[i][0][0];

1021

const double w2 = X[i][0][1];

1022

1023

// Compute affine mapping y = F(X)

1024

double y[2];

1025

y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];

1026

y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];

1027

1028

// Evaluate function at physical points

1029

double values[1];

1030

f.evaluate(values, y, c);

1031

1032

// Map function values using appropriate mapping

1033

// Affine map: Do nothing

1034

1035

// Note that we do not map the weights (yet).

1036

1037

// Take directional components

1038

for(int k = 0; k < 1; k++)

1039

result += values[k]*D[i][0][k];

1040

// Multiply by weights

1041

result *= W[i][0];

1042

1043

return result;

1044

}

1045

1046

/// Evaluate linear functionals for all dofs on the function f

1047

virtual void evaluate_dofs(double* values,

1048

const ufc::function& f,

1049

const ufc::cell& c) const

1050

{

1051

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

1052

}

1053

1054

/// Interpolate vertex values from dof values

1055

virtual void interpolate_vertex_values(double* vertex_values,

1056

const double* dof_values,

1057

const ufc::cell& c) const

1058

{

1059

// Evaluate at vertices and use affine mapping

1060

vertex_values[0] = dof_values[0];

1061

vertex_values[1] = dof_values[1];

1062

vertex_values[2] = dof_values[2];

1063

}

1064

1065

/// Return the number of sub elements (for a mixed element)

1066

virtual unsigned int num_sub_elements() const

1067

{

1068

return 1;

1069

}

1070

1071

/// Create a new finite element for sub element i (for a mixed element)

1072

virtual ufc::finite_element* create_sub_element(unsigned int i) const

1073

{

1074

return new optimization_0_finite_element_1();

1075

}

1076

1077

};

1078

1079

/// This class defines the interface for a local-to-global mapping of

1080

/// degrees of freedom (dofs).

1081

1082

class optimization_0_dof_map_0: public ufc::dof_map

1083

{

1084

private:

1085

1086

unsigned int __global_dimension;

1087

1088

public:

1089

1090

/// Constructor

1091

optimization_0_dof_map_0() : ufc::dof_map()

1092

{

1093

__global_dimension = 0;

1094

}

1095

1096

/// Destructor

1097

virtual ~optimization_0_dof_map_0()

1098

{

1099

// Do nothing

1100

}

1101

1102

/// Return a string identifying the dof map

1103

virtual const char* signature() const

1104

{

1105

return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3)";

1106

}

1107

1108

/// Return true iff mesh entities of topological dimension d are needed

1109

virtual bool needs_mesh_entities(unsigned int d) const

1110

{

1111

switch ( d )

1112

{

1113

case 0:

1114

return true;

1115

break;

1116

case 1:

1117

return true;

1118

break;

1119

case 2:

1120

return true;

1121

break;

1122

}

1123

return false;

1124

}

1125

1126

/// Initialize dof map for mesh (return true iff init_cell() is needed)

1127

virtual bool init_mesh(const ufc::mesh& m)

1128

{

1129

__global_dimension = m.num_entities[0] + 2*m.num_entities[1] + m.num_entities[2];

1130

return false;

1131

}

1132

1133

/// Initialize dof map for given cell

1134

virtual void init_cell(const ufc::mesh& m,

1135

const ufc::cell& c)

1136

{

1137

// Do nothing

1138

}

1139

1140

/// Finish initialization of dof map for cells

1141

virtual void init_cell_finalize()

1142

{

1143

// Do nothing

1144

}

1145

1146

/// Return the dimension of the global finite element function space

1147

virtual unsigned int global_dimension() const

1148

{

1149

return __global_dimension;

1150

}

1151

1152

/// Return the dimension of the local finite element function space for a cell

1153

virtual unsigned int local_dimension(const ufc::cell& c) const

1154

{

1155

return 10;

1156

}

1157

1158

/// Return the maximum dimension of the local finite element function space

1159

virtual unsigned int max_local_dimension() const

1160

{

1161

return 10;

1162

}

1163

1164

// Return the geometric dimension of the coordinates this dof map provides

1165

virtual unsigned int geometric_dimension() const

1166

{

1167

return 2;

1168

}

1169

1170

/// Return the number of dofs on each cell facet

1171

virtual unsigned int num_facet_dofs() const

1172

{

1173

return 4;

1174

}

1175

1176

/// Return the number of dofs associated with each cell entity of dimension d

1177

virtual unsigned int num_entity_dofs(unsigned int d) const

1178

{

1179

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

1180

}

1181

1182

/// Tabulate the local-to-global mapping of dofs on a cell

1183

virtual void tabulate_dofs(unsigned int* dofs,

1184

const ufc::mesh& m,

1185

const ufc::cell& c) const

1186

{

1187

dofs[0] = c.entity_indices[0][0];

1188

dofs[1] = c.entity_indices[0][1];

1189

dofs[2] = c.entity_indices[0][2];

1190

unsigned int offset = m.num_entities[0];

1191

dofs[3] = offset + 2*c.entity_indices[1][0];

1192

dofs[4] = offset + 2*c.entity_indices[1][0] + 1;

1193

dofs[5] = offset + 2*c.entity_indices[1][1];

1194

dofs[6] = offset + 2*c.entity_indices[1][1] + 1;

1195

dofs[7] = offset + 2*c.entity_indices[1][2];

1196

dofs[8] = offset + 2*c.entity_indices[1][2] + 1;

1197

offset = offset + 2*m.num_entities[1];

1198

dofs[9] = offset + c.entity_indices[2][0];

1199

}

1200

1201

/// Tabulate the local-to-local mapping from facet dofs to cell dofs

1202

virtual void tabulate_facet_dofs(unsigned int* dofs,

1203

unsigned int facet) const

1204

{

1205

switch ( facet )

1206

{

1207

case 0:

1208

dofs[0] = 1;

1209

dofs[1] = 2;

1210

dofs[2] = 3;

1211

dofs[3] = 4;

1212

break;

1213

case 1:

1214

dofs[0] = 0;

1215

dofs[1] = 2;

1216

dofs[2] = 5;

1217

dofs[3] = 6;

1218

break;

1219

case 2:

1220

dofs[0] = 0;

1221

dofs[1] = 1;

1222

dofs[2] = 7;

1223

dofs[3] = 8;

1224

break;

1225

}

1226

}

1227

1228

/// Tabulate the local-to-local mapping of dofs on entity (d, i)

1229

virtual void tabulate_entity_dofs(unsigned int* dofs,

1230

unsigned int d, unsigned int i) const

1231

{

1232

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

1233

}

1234

1235

/// Tabulate the coordinates of all dofs on a cell

1236

virtual void tabulate_coordinates(double** coordinates,

1237

const ufc::cell& c) const

1238

{

1239

const double * const * x = c.coordinates;

1240

coordinates[0][0] = x[0][0];

1241

coordinates[0][1] = x[0][1];

1242

coordinates[1][0] = x[1][0];

1243

coordinates[1][1] = x[1][1];

1244

coordinates[2][0] = x[2][0];

1245

coordinates[2][1] = x[2][1];

1246

coordinates[3][0] = 0.666666667*x[1][0] + 0.333333333*x[2][0];

1247

coordinates[3][1] = 0.666666667*x[1][1] + 0.333333333*x[2][1];

1248

coordinates[4][0] = 0.333333333*x[1][0] + 0.666666667*x[2][0];

1249

coordinates[4][1] = 0.333333333*x[1][1] + 0.666666667*x[2][1];

1250

coordinates[5][0] = 0.666666667*x[0][0] + 0.333333333*x[2][0];

1251

coordinates[5][1] = 0.666666667*x[0][1] + 0.333333333*x[2][1];

1252

coordinates[6][0] = 0.333333333*x[0][0] + 0.666666667*x[2][0];

1253

coordinates[6][1] = 0.333333333*x[0][1] + 0.666666667*x[2][1];

1254

coordinates[7][0] = 0.666666667*x[0][0] + 0.333333333*x[1][0];

1255

coordinates[7][1] = 0.666666667*x[0][1] + 0.333333333*x[1][1];

1256

coordinates[8][0] = 0.333333333*x[0][0] + 0.666666667*x[1][0];

1257

coordinates[8][1] = 0.333333333*x[0][1] + 0.666666667*x[1][1];

1258

coordinates[9][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];

1259

coordinates[9][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];

1260

}

1261

1262

/// Return the number of sub dof maps (for a mixed element)

1263

virtual unsigned int num_sub_dof_maps() const

1264

{

1265

return 1;

1266

}

1267

1268

/// Create a new dof_map for sub dof map i (for a mixed element)

1269

virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const

1270

{

1271

return new optimization_0_dof_map_0();

1272

}

1273

1274

};

1275

1276

/// This class defines the interface for a local-to-global mapping of

1277

/// degrees of freedom (dofs).

1278

1279

class optimization_0_dof_map_1: public ufc::dof_map

1280

{

1281

private:

1282

1283

unsigned int __global_dimension;

1284

1285

public:

1286

1287

/// Constructor

1288

optimization_0_dof_map_1() : ufc::dof_map()

1289

{

1290

__global_dimension = 0;

1291

}

1292

1293

/// Destructor

1294

virtual ~optimization_0_dof_map_1()

1295

{

1296

// Do nothing

1297

}

1298

1299

/// Return a string identifying the dof map

1300

virtual const char* signature() const

1301

{

1302

return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3)";

1303

}

1304

1305

/// Return true iff mesh entities of topological dimension d are needed

1306

virtual bool needs_mesh_entities(unsigned int d) const

1307

{

1308

switch ( d )

1309

{

1310

case 0:

1311

return true;

1312

break;

1313

case 1:

1314

return true;

1315

break;

1316

case 2:

1317

return true;

1318

break;

1319

}

1320

return false;

1321

}

1322

1323

/// Initialize dof map for mesh (return true iff init_cell() is needed)

1324

virtual bool init_mesh(const ufc::mesh& m)

1325

{

1326

__global_dimension = m.num_entities[0] + 2*m.num_entities[1] + m.num_entities[2];

1327

return false;

1328

}

1329

1330

/// Initialize dof map for given cell

1331

virtual void init_cell(const ufc::mesh& m,

1332

const ufc::cell& c)

1333

{

1334

// Do nothing

1335

}

1336

1337

/// Finish initialization of dof map for cells

1338

virtual void init_cell_finalize()

1339

{

1340

// Do nothing

1341

}

1342

1343

/// Return the dimension of the global finite element function space

1344

virtual unsigned int global_dimension() const

1345

{

1346

return __global_dimension;

1347

}

1348

1349

/// Return the dimension of the local finite element function space for a cell

1350

virtual unsigned int local_dimension(const ufc::cell& c) const

1351

{

1352

return 10;

1353

}

1354

1355

/// Return the maximum dimension of the local finite element function space

1356

virtual unsigned int max_local_dimension() const

1357

{

1358

return 10;

1359

}

1360

1361

// Return the geometric dimension of the coordinates this dof map provides

1362

virtual unsigned int geometric_dimension() const

1363

{

1364

return 2;

1365

}

1366

1367

/// Return the number of dofs on each cell facet

1368

virtual unsigned int num_facet_dofs() const

1369

{

1370

return 4;

1371

}

1372

1373

/// Return the number of dofs associated with each cell entity of dimension d

1374

virtual unsigned int num_entity_dofs(unsigned int d) const

1375

{

1376

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

1377

}

1378

1379

/// Tabulate the local-to-global mapping of dofs on a cell

1380

virtual void tabulate_dofs(unsigned int* dofs,

1381

const ufc::mesh& m,

1382

const ufc::cell& c) const

1383

{

1384

dofs[0] = c.entity_indices[0][0];

1385

dofs[1] = c.entity_indices[0][1];

1386

dofs[2] = c.entity_indices[0][2];

1387

unsigned int offset = m.num_entities[0];

1388

dofs[3] = offset + 2*c.entity_indices[1][0];

1389

dofs[4] = offset + 2*c.entity_indices[1][0] + 1;

1390

dofs[5] = offset + 2*c.entity_indices[1][1];

1391

dofs[6] = offset + 2*c.entity_indices[1][1] + 1;

1392

dofs[7] = offset + 2*c.entity_indices[1][2];

1393

dofs[8] = offset + 2*c.entity_indices[1][2] + 1;

1394

offset = offset + 2*m.num_entities[1];

1395

dofs[9] = offset + c.entity_indices[2][0];

1396

}

1397

1398

/// Tabulate the local-to-local mapping from facet dofs to cell dofs

1399

virtual void tabulate_facet_dofs(unsigned int* dofs,

1400

unsigned int facet) const

1401

{

1402

switch ( facet )

1403

{

1404

case 0:

1405

dofs[0] = 1;

1406

dofs[1] = 2;

1407

dofs[2] = 3;

1408

dofs[3] = 4;

1409

break;

1410

case 1:

1411

dofs[0] = 0;

1412

dofs[1] = 2;

1413

dofs[2] = 5;

1414

dofs[3] = 6;

1415

break;

1416

case 2:

1417

dofs[0] = 0;

1418

dofs[1] = 1;

1419

dofs[2] = 7;

1420

dofs[3] = 8;

1421

break;

1422

}

1423

}

1424

1425

/// Tabulate the local-to-local mapping of dofs on entity (d, i)

1426

virtual void tabulate_entity_dofs(unsigned int* dofs,

1427

unsigned int d, unsigned int i) const

1428

{

1429

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

1430

}

1431

1432

/// Tabulate the coordinates of all dofs on a cell

1433

virtual void tabulate_coordinates(double** coordinates,

1434

const ufc::cell& c) const

1435

{

1436

const double * const * x = c.coordinates;

1437

coordinates[0][0] = x[0][0];

1438

coordinates[0][1] = x[0][1];

1439

coordinates[1][0] = x[1][0];

1440

coordinates[1][1] = x[1][1];

1441

coordinates[2][0] = x[2][0];

1442

coordinates[2][1] = x[2][1];

1443

coordinates[3][0] = 0.666666667*x[1][0] + 0.333333333*x[2][0];

1444

coordinates[3][1] = 0.666666667*x[1][1] + 0.333333333*x[2][1];

1445

coordinates[4][0] = 0.333333333*x[1][0] + 0.666666667*x[2][0];

1446

coordinates[4][1] = 0.333333333*x[1][1] + 0.666666667*x[2][1];

1447

coordinates[5][0] = 0.666666667*x[0][0] + 0.333333333*x[2][0];

1448

coordinates[5][1] = 0.666666667*x[0][1] + 0.333333333*x[2][1];

1449

coordinates[6][0] = 0.333333333*x[0][0] + 0.666666667*x[2][0];

1450

coordinates[6][1] = 0.333333333*x[0][1] + 0.666666667*x[2][1];

1451

coordinates[7][0] = 0.666666667*x[0][0] + 0.333333333*x[1][0];

1452

coordinates[7][1] = 0.666666667*x[0][1] + 0.333333333*x[1][1];

1453

coordinates[8][0] = 0.333333333*x[0][0] + 0.666666667*x[1][0];

1454

coordinates[8][1] = 0.333333333*x[0][1] + 0.666666667*x[1][1];

1455

coordinates[9][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];

1456

coordinates[9][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];

1457

}

1458

1459

/// Return the number of sub dof maps (for a mixed element)

1460

virtual unsigned int num_sub_dof_maps() const

1461

{

1462

return 1;

1463

}

1464

1465

/// Create a new dof_map for sub dof map i (for a mixed element)

1466

virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const

1467

{

1468

return new optimization_0_dof_map_1();

1469

}

1470

1471

};

1472

1473

/// This class defines the interface for the tabulation of the cell

1474

/// tensor corresponding to the local contribution to a form from

1475

/// the integral over a cell.

1476

1477

class optimization_0_cell_integral_0_quadrature: public ufc::cell_integral

1478

{

1479

public:

1480

1481

/// Constructor

1482

optimization_0_cell_integral_0_quadrature() : ufc::cell_integral()

1483

{

1484

// Do nothing

1485

}

1486

1487

/// Destructor

1488

virtual ~optimization_0_cell_integral_0_quadrature()

1489

{

1490

// Do nothing

1491

}

1492

1493

/// Tabulate the tensor for the contribution from a local cell

1494

virtual void tabulate_tensor(double* A,

1495

const double * const * w,

1496

const ufc::cell& c) const

1497

{

1498

// Extract vertex coordinates

1499

const double * const * x = c.coordinates;

1500

1501

// Compute Jacobian of affine map from reference cell

1502

const double J_00 = x[1][0] - x[0][0];

1503

const double J_01 = x[2][0] - x[0][0];

1504

const double J_10 = x[1][1] - x[0][1];

1505

const double J_11 = x[2][1] - x[0][1];

1506

1507

// Compute determinant of Jacobian

1508

double detJ = J_00*J_11 - J_01*J_10;

1509

1510

// Compute inverse of Jacobian

1511

const double Jinv_00 = J_11 / detJ;

1512

const double Jinv_01 = -J_01 / detJ;

1513

const double Jinv_10 = -J_10 / detJ;

1514

const double Jinv_11 = J_00 / detJ;

1515

1516

// Set scale factor

1517

const double det = std::abs(detJ);

1518

1519

1520

// Array of quadrature weights

1521

static const double W9[9] = {0.0558144205, 0.0636780851, 0.0193963833, 0.0893030728, 0.101884936, 0.0310342133, 0.0558144205, 0.0636780851, 0.0193963833};

1522

// Quadrature points on the UFC reference element: (0.102717655, 0.0885879595), (0.0665540678, 0.409466864), (0.0239311323, 0.787659462), (0.45570602, 0.0885879595), (0.295266568, 0.409466864), (0.106170269, 0.787659462), (0.808694386, 0.0885879595), (0.523979068, 0.409466864), (0.188409406, 0.787659462)

1523

1524

// Value of basis functions at quadrature points.

1525

static const double FE0_D01[9][10] = \

1526

{{-2.55056976, 0, 0.308654023, -0.319792072, -0.216541666, 3.6540445, -1.41212877, -1.7805847, 0.319792072, 1.99712637},

1527

{0.00933175346, 0, -0.421749753, -0.239695812, 0.436302203, -2.60173084, 3.01414884, -0.642076032, 0.239695812, 0.20577383},

1528

{0.216460246, 0, 2.28656512, -0.0999586575, 0.40124864, -0.831016526, -1.67200884, -0.014048866, 0.0999586575, -0.387199774},

1529

{0.297836494, 0, 0.308654023, 0.752840597, -0.960685296, 0.06149462, -0.667985137, -3.55635828, -0.752840597, 4.51704358},

1530

{0.480437438, 0, -0.421749753, -0.151737883, 1.93565109, -1.57348764, 1.51479995, -1.02522379, 0.151737883, -0.9104273},

1531

{-0.19664128, 0, 2.28656512, -0.325592509, 1.78013625, 0.960972609, -3.05089645, 0.173418808, 0.325592509, -1.95355506},

1532

{-0.217978481, 0, 0.308654023, 5.18969449, -1.70482893, -0.166834036, 0.0761584935, 1.39631059, -5.18969449, 0.308518342},

1533

{-0.460810883, 0, -0.421749753, 1.34857405, 3.43499997, 0.86710957, 0.015451066, 1.41633647, -1.34857405, -4.85133644},

1534

{-0.792351247, 0, 2.28656512, -0.368617919, 3.15902386, 2.93557019, -4.42978406, 0.726103365, 0.368617919, -3.88512722}};

1535

1536

static const double FE0_D10[9][10] = \

1537

{{-2.55056976, 0.217978481, 0, -0.152958037, -0.292700159, -1.53564999, 0.292700159, 3.40910979, -1.07651851, 1.68860803},

1538

{0.00933175346, 0.460810883, 0, -1.10680538, 0.420851137, -3.95030489, -0.420851137, 0.706498019, -1.17664066, 5.05711027},

1539

{0.216460246, 0.792351247, 0, -3.03552884, 4.8310327, -0.462398607, -4.8310327, -0.382666785, -0.626144708, 3.49792745},

1540

{0.297836494, -0.297836494, 0, 0.691345977, -0.292700159, -0.691345977, 0.292700159, -2.80351769, 2.80351769, 0},

1541

{0.480437438, -0.480437438, 0, 1.42174975, 0.420851137, -1.42174975, -0.420851137, -1.17696167, 1.17696167, 0},

1542

{-0.19664128, 0.19664128, 0, -1.28656512, 4.8310327, 1.28656512, -4.8310327, -0.152173702, 0.152173702, 0},

1543

{-0.217978481, 2.55056976, 0, 1.53564999, -0.292700159, 0.152958037, 0.292700159, 1.07651851, -3.40910979, -1.68860803},

1544

{-0.460810883, -0.00933175346, 0, 3.95030489, 0.420851137, 1.10680538, -0.420851137, 1.17664066, -0.706498019, -5.05711027},

1545

{-0.792351247, -0.216460246, 0, 0.462398607, 4.8310327, 3.03552884, -4.8310327, 0.626144708, 0.382666785, -3.49792745}};

1546

1547

1548

// Compute element tensor using UFL quadrature representation

1549

// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)

1550

// Total number of operations to compute element tensor: 16200

1551

1552

// Loop quadrature points for integral

1553

// Number of operations to compute element tensor for following IP loop = 16200

1554

for (unsigned int ip = 0; ip < 9; ip++)

1555

{

1556

1557

// Number of operations for primary indices: 1800

1558

for (unsigned int j = 0; j < 10; j++)

1559

{

1560

for (unsigned int k = 0; k < 10; k++)

1561

{

1562

// Number of operations to compute entry: 18

1563

A[j*10 + k] += ((Jinv_01*FE0_D10[ip][j] + Jinv_11*FE0_D01[ip][j])*(Jinv_01*FE0_D10[ip][k] + Jinv_11*FE0_D01[ip][k]) + (Jinv_00*FE0_D10[ip][j] + Jinv_10*FE0_D01[ip][j])*(Jinv_00*FE0_D10[ip][k] + Jinv_10*FE0_D01[ip][k]))*W9[ip]*det;

1564

}// end loop over 'k'

1565

}// end loop over 'j'

1566

}// end loop over 'ip'

1567

}

1568

1569

};

1570

1571

/// This class defines the interface for the tabulation of the cell

1572

/// tensor corresponding to the local contribution to a form from

1573

/// the integral over a cell.

1574

1575

class optimization_0_cell_integral_0: public ufc::cell_integral

1576

{

1577

private:

1578

1579

optimization_0_cell_integral_0_quadrature integral_0_quadrature;

1580

1581

public:

1582

1583

/// Constructor

1584

optimization_0_cell_integral_0() : ufc::cell_integral()

1585

{

1586

// Do nothing

1587

}

1588

1589

/// Destructor

1590

virtual ~optimization_0_cell_integral_0()

1591

{

1592

// Do nothing

1593

}

1594

1595

/// Tabulate the tensor for the contribution from a local cell

1596

virtual void tabulate_tensor(double* A,

1597

const double * const * w,

1598

const ufc::cell& c) const

1599

{

1600

// Reset values of the element tensor block

1601

for (unsigned int j = 0; j < 100; j++)

1602

A[j] = 0;

1603

1604

// Add all contributions to element tensor

1605

integral_0_quadrature.tabulate_tensor(A, w, c);

1606

}

1607

1608

};

1609

1610

/// This class defines the interface for the assembly of the global

1611

/// tensor corresponding to a form with r + n arguments, that is, a

1612

/// mapping

1613

///

1614

/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R

1615

///

1616

/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r

1617

/// global tensor A is defined by

1618

///

1619

/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),

1620

///

1621

/// where each argument Vj represents the application to the

1622

/// sequence of basis functions of Vj and w1, w2, ..., wn are given

1623

/// fixed functions (coefficients).

1624

1625

class optimization_form_0: public ufc::form

1626

{

1627

public:

1628

1629

/// Constructor

1630

optimization_form_0() : ufc::form()

1631

{

1632

// Do nothing

1633

}

1634

1635

/// Destructor

1636

virtual ~optimization_form_0()

1637

{

1638

// Do nothing

1639

}

1640

1641

/// Return a string identifying the form

1642

virtual const char* signature() const

1643

{

1644

return "Form([Integral(IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})), Measure('cell', 0, None))])";

1645

}

1646

1647

/// Return the rank of the global tensor (r)

1648

virtual unsigned int rank() const

1649

{

1650

return 2;

1651

}

1652

1653

/// Return the number of coefficients (n)

1654

virtual unsigned int num_coefficients() const

1655

{

1656

return 0;

1657

}

1658

1659

/// Return the number of cell integrals

1660

virtual unsigned int num_cell_integrals() const

1661

{

1662

return 1;

1663

}

1664

1665

/// Return the number of exterior facet integrals

1666

virtual unsigned int num_exterior_facet_integrals() const

1667

{

1668

return 0;

1669

}

1670

1671

/// Return the number of interior facet integrals

1672

virtual unsigned int num_interior_facet_integrals() const

1673

{

1674

return 0;

1675

}

1676

1677

/// Create a new finite element for argument function i

1678

virtual ufc::finite_element* create_finite_element(unsigned int i) const

1679

{

1680

switch ( i )

1681

{

1682

case 0:

1683

return new optimization_0_finite_element_0();

1684

break;

1685

case 1:

1686

return new optimization_0_finite_element_1();

1687

break;

1688

}

1689

return 0;

1690

}

1691

1692

/// Create a new dof map for argument function i

1693

virtual ufc::dof_map* create_dof_map(unsigned int i) const

1694

{

1695

switch ( i )

1696

{

1697

case 0:

1698

return new optimization_0_dof_map_0();

1699

break;

1700

case 1:

1701

return new optimization_0_dof_map_1();

1702

break;

1703

}

1704

return 0;

1705

}

1706

1707

/// Create a new cell integral on sub domain i

1708

virtual ufc::cell_integral* create_cell_integral(unsigned int i) const

1709

{

1710

return new optimization_0_cell_integral_0();

1711

}

1712

1713

/// Create a new exterior facet integral on sub domain i

1714

virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const

1715

{

1716

return 0;

1717

}

1718

1719

/// Create a new interior facet integral on sub domain i

1720

virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const

1721

{

1722

return 0;

1723

}

1724

1725

};

1726

1727

/// This class defines the interface for a finite element.

1728

1729

class optimization_1_finite_element_0: public ufc::finite_element

1730

{

1731

public:

1732

1733

/// Constructor

1734

optimization_1_finite_element_0() : ufc::finite_element()

1735

{

1736

// Do nothing

1737

}

1738

1739

/// Destructor

1740

virtual ~optimization_1_finite_element_0()

1741

{

1742

// Do nothing

1743

}

1744

1745

/// Return a string identifying the finite element

1746

virtual const char* signature() const

1747

{

1748

return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3)";

1749

}

1750

1751

/// Return the cell shape

1752

virtual ufc::shape cell_shape() const

1753

{

1754

return ufc::triangle;

1755

}

1756

1757

/// Return the dimension of the finite element function space

1758

virtual unsigned int space_dimension() const

1759

{

1760

return 10;

1761

}

1762

1763

/// Return the rank of the value space

1764

virtual unsigned int value_rank() const

1765

{

1766

return 0;

1767

}

1768

1769

/// Return the dimension of the value space for axis i

1770

virtual unsigned int value_dimension(unsigned int i) const

1771

{

1772

return 1;

1773

}

1774

1775

/// Evaluate basis function i at given point in cell

1776

virtual void evaluate_basis(unsigned int i,

1777

double* values,

1778

const double* coordinates,

1779

const ufc::cell& c) const

1780

{

1781

// Extract vertex coordinates

1782

const double * const * element_coordinates = c.coordinates;

1783

1784

// Compute Jacobian of affine map from reference cell

1785

const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];

1786

const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];

1787

const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];

1788

const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];

1789

1790

// Compute determinant of Jacobian

1791

const double detJ = J_00*J_11 - J_01*J_10;

1792

1793

// Compute inverse of Jacobian

1794

1795

// Get coordinates and map to the reference (UFC) element

1796

double x = (element_coordinates[0][1]*element_coordinates[2][0] -\

1797

element_coordinates[0][0]*element_coordinates[2][1] +\

1798

J_11*coordinates[0] - J_01*coordinates[1]) / detJ;

1799

double y = (element_coordinates[1][1]*element_coordinates[0][0] -\

1800

element_coordinates[1][0]*element_coordinates[0][1] -\

1801

J_10*coordinates[0] + J_00*coordinates[1]) / detJ;

1802

1803

// Map coordinates to the reference square

1804

if (std::abs(y - 1.0) < 1e-08)

1805

x = -1.0;

1806

else

1807

x = 2.0 *x/(1.0 - y) - 1.0;

1808

y = 2.0*y - 1.0;

1809

1810

// Reset values

1811

*values = 0;

1812

1813

// Map degree of freedom to element degree of freedom

1814

const unsigned int dof = i;

1815

1816

// Generate scalings

1817

const double scalings_y_0 = 1;

1818

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

1819

const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);

1820

const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);

1821

1822

// Compute psitilde_a

1823

const double psitilde_a_0 = 1;

1824

const double psitilde_a_1 = x;

1825

const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;

1826

const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;

1827

1828

// Compute psitilde_bs

1829

const double psitilde_bs_0_0 = 1;

1830

const double psitilde_bs_0_1 = 1.5*y + 0.5;

1831

const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;

1832

const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;

1833

const double psitilde_bs_1_0 = 1;

1834

const double psitilde_bs_1_1 = 2.5*y + 1.5;

1835

const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;

1836

const double psitilde_bs_2_0 = 1;

1837

const double psitilde_bs_2_1 = 3.5*y + 2.5;

1838

const double psitilde_bs_3_0 = 1;

1839

1840

// Compute basisvalues

1841

const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;

1842

const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;

1843

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

1844

const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;

1845

const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;

1846

const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;

1847

const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;

1848

const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;

1849

const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;

1850

const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;

1851

1852

// Table(s) of coefficients

1853

static const double coefficients0[10][10] = \

1854

{{0.0471404521, -0.0288675135, -0.0166666667, 0.0782460796, 0.0606091527, 0.0349927106, -0.0601337794, -0.0508223195, -0.0393667994, -0.0227284323},

1855

{0.0471404521, 0.0288675135, -0.0166666667, 0.0782460796, -0.0606091527, 0.0349927106, 0.0601337794, -0.0508223195, 0.0393667994, -0.0227284323},

1856

{0.0471404521, 0, 0.0333333333, 0, 0, 0.104978132, 0, 0, 0, 0.090913729},

1857

{0.106066017, 0.259807621, -0.15, 0.117369119, 0.0606091527, -0.0787335989, 0, 0.101644639, -0.131222665, 0.090913729},

1858

{0.106066017, 0, 0.3, 0, 0.151522882, 0.026244533, 0, 0, 0.131222665, -0.136370594},

1859

{0.106066017, -0.259807621, -0.15, 0.117369119, -0.0606091527, -0.0787335989, 0, 0.101644639, 0.131222665, 0.090913729},

1860

{0.106066017, 0, 0.3, 0, -0.151522882, 0.026244533, 0, 0, -0.131222665, -0.136370594},

1861

{0.106066017, -0.259807621, -0.15, -0.0782460796, 0.090913729, 0.0962299542, 0.180401338, 0.0508223195, -0.0131222665, -0.0227284323},

1862

{0.106066017, 0.259807621, -0.15, -0.0782460796, -0.090913729, 0.0962299542, -0.180401338, 0.0508223195, 0.0131222665, -0.0227284323},

1863

{0.636396103, 0, 0, -0.234738239, 0, -0.26244533, 0, -0.203289278, 0, 0.090913729}};

1864

1865

// Extract relevant coefficients

1866

const double coeff0_0 = coefficients0[dof][0];

1867

const double coeff0_1 = coefficients0[dof][1];

1868

const double coeff0_2 = coefficients0[dof][2];

1869

const double coeff0_3 = coefficients0[dof][3];

1870

const double coeff0_4 = coefficients0[dof][4];

1871

const double coeff0_5 = coefficients0[dof][5];

1872

const double coeff0_6 = coefficients0[dof][6];

1873

const double coeff0_7 = coefficients0[dof][7];

1874

const double coeff0_8 = coefficients0[dof][8];

1875

const double coeff0_9 = coefficients0[dof][9];

1876

1877

// Compute value(s)

1878

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9;

1879

}

1880

1881

/// Evaluate all basis functions at given point in cell

1882

virtual void evaluate_basis_all(double* values,

1883

const double* coordinates,

1884

const ufc::cell& c) const

1885

{

1886

throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");

1887

}

1888

1889

/// Evaluate order n derivatives of basis function i at given point in cell

1890

virtual void evaluate_basis_derivatives(unsigned int i,

1891

unsigned int n,

1892

double* values,

1893

const double* coordinates,

1894

const ufc::cell& c) const

1895

{

1896

// Extract vertex coordinates

1897

const double * const * element_coordinates = c.coordinates;

1898

1899

// Compute Jacobian of affine map from reference cell

1900

const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];

1901

const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];

1902

const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];

1903

const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];

1904

1905

// Compute determinant of Jacobian

1906

const double detJ = J_00*J_11 - J_01*J_10;

1907

1908

// Compute inverse of Jacobian

1909

1910

// Get coordinates and map to the reference (UFC) element

1911

double x = (element_coordinates[0][1]*element_coordinates[2][0] -\

1912

element_coordinates[0][0]*element_coordinates[2][1] +\

1913

J_11*coordinates[0] - J_01*coordinates[1]) / detJ;

1914

double y = (element_coordinates[1][1]*element_coordinates[0][0] -\

1915

element_coordinates[1][0]*element_coordinates[0][1] -\

1916

J_10*coordinates[0] + J_00*coordinates[1]) / detJ;

1917

1918

// Map coordinates to the reference square

1919

if (std::abs(y - 1.0) < 1e-08)

1920

x = -1.0;

1921

else

1922

x = 2.0 *x/(1.0 - y) - 1.0;

1923

y = 2.0*y - 1.0;

1924

1925

// Compute number of derivatives

1926

unsigned int num_derivatives = 1;

1927

1928

for (unsigned int j = 0; j < n; j++)

1929

num_derivatives *= 2;

1930

1931

1932

// Declare pointer to two dimensional array that holds combinations of derivatives and initialise

1933

unsigned int **combinations = new unsigned int *[num_derivatives];

1934

1935

for (unsigned int j = 0; j < num_derivatives; j++)

1936

{

1937

combinations[j] = new unsigned int [n];

1938

for (unsigned int k = 0; k < n; k++)

1939

combinations[j][k] = 0;

1940

}

1941

1942

// Generate combinations of derivatives

1943

for (unsigned int row = 1; row < num_derivatives; row++)

1944

{

1945

for (unsigned int num = 0; num < row; num++)

1946

{

1947

for (unsigned int col = n-1; col+1 > 0; col--)

1948

{

1949

if (combinations[row][col] + 1 > 1)

1950

combinations[row][col] = 0;

1951

else

1952

{

1953

combinations[row][col] += 1;

1954

break;

1955

}

1956

}

1957

}

1958

}

1959

1960

// Compute inverse of Jacobian

1961

const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};

1962

1963

// Declare transformation matrix

1964

// Declare pointer to two dimensional array and initialise

1965

double **transform = new double *[num_derivatives];

1966

1967

for (unsigned int j = 0; j < num_derivatives; j++)

1968

{

1969

transform[j] = new double [num_derivatives];

1970

for (unsigned int k = 0; k < num_derivatives; k++)

1971

transform[j][k] = 1;

1972

}

1973

1974

// Construct transformation matrix

1975

for (unsigned int row = 0; row < num_derivatives; row++)

1976

{

1977

for (unsigned int col = 0; col < num_derivatives; col++)

1978

{

1979

for (unsigned int k = 0; k < n; k++)

1980

transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];

1981

}

1982

}

1983

1984

// Reset values

1985

for (unsigned int j = 0; j < 1*num_derivatives; j++)

1986

values[j] = 0;

1987

1988

// Map degree of freedom to element degree of freedom

1989

const unsigned int dof = i;

1990

1991

// Generate scalings

1992

const double scalings_y_0 = 1;

1993

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

1994

const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);

1995

const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);

1996

1997

// Compute psitilde_a

1998

const double psitilde_a_0 = 1;

1999

const double psitilde_a_1 = x;

2000

const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;

2001

const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;

2002

2003

// Compute psitilde_bs

2004

const double psitilde_bs_0_0 = 1;

2005

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2006

const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;

2007

const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;

2008

const double psitilde_bs_1_0 = 1;

2009

const double psitilde_bs_1_1 = 2.5*y + 1.5;

2010

const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;

2011

const double psitilde_bs_2_0 = 1;

2012

const double psitilde_bs_2_1 = 3.5*y + 2.5;

2013

const double psitilde_bs_3_0 = 1;

2014

2015

// Compute basisvalues

2016

const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;

2017

const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;

2018

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2019

const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;

2020

const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;

2021

const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;

2022

const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;

2023

const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;

2024

const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;

2025

const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;

2026

2027

// Table(s) of coefficients

2028

static const double coefficients0[10][10] = \

2029

{{0.0471404521, -0.0288675135, -0.0166666667, 0.0782460796, 0.0606091527, 0.0349927106, -0.0601337794, -0.0508223195, -0.0393667994, -0.0227284323},

2030

{0.0471404521, 0.0288675135, -0.0166666667, 0.0782460796, -0.0606091527, 0.0349927106, 0.0601337794, -0.0508223195, 0.0393667994, -0.0227284323},

2031

{0.0471404521, 0, 0.0333333333, 0, 0, 0.104978132, 0, 0, 0, 0.090913729},

2032

{0.106066017, 0.259807621, -0.15, 0.117369119, 0.0606091527, -0.0787335989, 0, 0.101644639, -0.131222665, 0.090913729},

2033

{0.106066017, 0, 0.3, 0, 0.151522882, 0.026244533, 0, 0, 0.131222665, -0.136370594},

2034

{0.106066017, -0.259807621, -0.15, 0.117369119, -0.0606091527, -0.0787335989, 0, 0.101644639, 0.131222665, 0.090913729},

2035

{0.106066017, 0, 0.3, 0, -0.151522882, 0.026244533, 0, 0, -0.131222665, -0.136370594},

2036

{0.106066017, -0.259807621, -0.15, -0.0782460796, 0.090913729, 0.0962299542, 0.180401338, 0.0508223195, -0.0131222665, -0.0227284323},

2037

{0.106066017, 0.259807621, -0.15, -0.0782460796, -0.090913729, 0.0962299542, -0.180401338, 0.0508223195, 0.0131222665, -0.0227284323},

2038

{0.636396103, 0, 0, -0.234738239, 0, -0.26244533, 0, -0.203289278, 0, 0.090913729}};

2039

2040

// Interesting (new) part

2041

// Tables of derivatives of the polynomial base (transpose)

2042

static const double dmats0[10][10] = \

2043

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2044

{4.89897949, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2045

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2046

{0, 9.48683298, 0, 0, 0, 0, 0, 0, 0, 0},

2047

{4, 0, 7.07106781, 0, 0, 0, 0, 0, 0, 0},

2048

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2049

{5.29150262, 0, -2.99332591, 13.662601, 0, 0.611010093, 0, 0, 0, 0},

2050

{0, 4.38178046, 0, 0, 12.5219807, 0, 0, 0, 0, 0},

2051

{3.46410162, 0, 7.83836718, 0, 0, 8.4, 0, 0, 0, 0},

2052

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};

2053

2054

static const double dmats1[10][10] = \

2055

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2056

{2.44948974, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2057

{4.24264069, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2058

{2.5819889, 4.74341649, -0.912870929, 0, 0, 0, 0, 0, 0, 0},

2059

{2, 6.12372436, 3.53553391, 0, 0, 0, 0, 0, 0, 0},

2060

{-2.30940108, 0, 8.16496581, 0, 0, 0, 0, 0, 0, 0},

2061

{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.305505046, 0, 0, 0, 0},

2062

{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0, 0, 0, 0},

2063

{1.73205081, -5.09116882, 3.91918359, 0, 9.69948452, 4.2, 0, 0, 0, 0},

2064

{5, 0, -2.82842712, 0, 0, 12.1243557, 0, 0, 0, 0}};

2065

2066

// Compute reference derivatives

2067

// Declare pointer to array of derivatives on FIAT element

2068

double *derivatives = new double [num_derivatives];

2069

2070

// Declare coefficients

2071

double coeff0_0 = 0;

2072

double coeff0_1 = 0;

2073

double coeff0_2 = 0;

2074

double coeff0_3 = 0;

2075

double coeff0_4 = 0;

2076

double coeff0_5 = 0;

2077

double coeff0_6 = 0;

2078

double coeff0_7 = 0;

2079

double coeff0_8 = 0;

2080

double coeff0_9 = 0;

2081

2082

// Declare new coefficients

2083

double new_coeff0_0 = 0;

2084

double new_coeff0_1 = 0;

2085

double new_coeff0_2 = 0;

2086

double new_coeff0_3 = 0;

2087

double new_coeff0_4 = 0;

2088

double new_coeff0_5 = 0;

2089

double new_coeff0_6 = 0;

2090

double new_coeff0_7 = 0;

2091

double new_coeff0_8 = 0;

2092

double new_coeff0_9 = 0;

2093

2094

// Loop possible derivatives

2095

for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)

2096

{

2097

// Get values from coefficients array

2098

new_coeff0_0 = coefficients0[dof][0];

2099

new_coeff0_1 = coefficients0[dof][1];

2100

new_coeff0_2 = coefficients0[dof][2];

2101

new_coeff0_3 = coefficients0[dof][3];

2102

new_coeff0_4 = coefficients0[dof][4];

2103

new_coeff0_5 = coefficients0[dof][5];

2104

new_coeff0_6 = coefficients0[dof][6];

2105

new_coeff0_7 = coefficients0[dof][7];

2106

new_coeff0_8 = coefficients0[dof][8];

2107

new_coeff0_9 = coefficients0[dof][9];

2108

2109

// Loop derivative order

2110

for (unsigned int j = 0; j < n; j++)

2111

{

2112

// Update old coefficients

2113

coeff0_0 = new_coeff0_0;

2114

coeff0_1 = new_coeff0_1;

2115

coeff0_2 = new_coeff0_2;

2116

coeff0_3 = new_coeff0_3;

2117

coeff0_4 = new_coeff0_4;

2118

coeff0_5 = new_coeff0_5;

2119

coeff0_6 = new_coeff0_6;

2120

coeff0_7 = new_coeff0_7;

2121

coeff0_8 = new_coeff0_8;

2122

coeff0_9 = new_coeff0_9;

2123

2124

if(combinations[deriv_num][j] == 0)

2125

{

2126

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0] + coeff0_6*dmats0[6][0] + coeff0_7*dmats0[7][0] + coeff0_8*dmats0[8][0] + coeff0_9*dmats0[9][0];

2127

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1] + coeff0_6*dmats0[6][1] + coeff0_7*dmats0[7][1] + coeff0_8*dmats0[8][1] + coeff0_9*dmats0[9][1];

2128

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2] + coeff0_6*dmats0[6][2] + coeff0_7*dmats0[7][2] + coeff0_8*dmats0[8][2] + coeff0_9*dmats0[9][2];

2129

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3] + coeff0_6*dmats0[6][3] + coeff0_7*dmats0[7][3] + coeff0_8*dmats0[8][3] + coeff0_9*dmats0[9][3];

2130

new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4] + coeff0_6*dmats0[6][4] + coeff0_7*dmats0[7][4] + coeff0_8*dmats0[8][4] + coeff0_9*dmats0[9][4];

2131

new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5] + coeff0_6*dmats0[6][5] + coeff0_7*dmats0[7][5] + coeff0_8*dmats0[8][5] + coeff0_9*dmats0[9][5];

2132

new_coeff0_6 = coeff0_0*dmats0[0][6] + coeff0_1*dmats0[1][6] + coeff0_2*dmats0[2][6] + coeff0_3*dmats0[3][6] + coeff0_4*dmats0[4][6] + coeff0_5*dmats0[5][6] + coeff0_6*dmats0[6][6] + coeff0_7*dmats0[7][6] + coeff0_8*dmats0[8][6] + coeff0_9*dmats0[9][6];

2133

new_coeff0_7 = coeff0_0*dmats0[0][7] + coeff0_1*dmats0[1][7] + coeff0_2*dmats0[2][7] + coeff0_3*dmats0[3][7] + coeff0_4*dmats0[4][7] + coeff0_5*dmats0[5][7] + coeff0_6*dmats0[6][7] + coeff0_7*dmats0[7][7] + coeff0_8*dmats0[8][7] + coeff0_9*dmats0[9][7];

2134

new_coeff0_8 = coeff0_0*dmats0[0][8] + coeff0_1*dmats0[1][8] + coeff0_2*dmats0[2][8] + coeff0_3*dmats0[3][8] + coeff0_4*dmats0[4][8] + coeff0_5*dmats0[5][8] + coeff0_6*dmats0[6][8] + coeff0_7*dmats0[7][8] + coeff0_8*dmats0[8][8] + coeff0_9*dmats0[9][8];

2135

new_coeff0_9 = coeff0_0*dmats0[0][9] + coeff0_1*dmats0[1][9] + coeff0_2*dmats0[2][9] + coeff0_3*dmats0[3][9] + coeff0_4*dmats0[4][9] + coeff0_5*dmats0[5][9] + coeff0_6*dmats0[6][9] + coeff0_7*dmats0[7][9] + coeff0_8*dmats0[8][9] + coeff0_9*dmats0[9][9];

2136

}

2137

if(combinations[deriv_num][j] == 1)

2138

{

2139

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0] + coeff0_6*dmats1[6][0] + coeff0_7*dmats1[7][0] + coeff0_8*dmats1[8][0] + coeff0_9*dmats1[9][0];

2140

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1] + coeff0_6*dmats1[6][1] + coeff0_7*dmats1[7][1] + coeff0_8*dmats1[8][1] + coeff0_9*dmats1[9][1];

2141

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2] + coeff0_6*dmats1[6][2] + coeff0_7*dmats1[7][2] + coeff0_8*dmats1[8][2] + coeff0_9*dmats1[9][2];

2142

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3] + coeff0_6*dmats1[6][3] + coeff0_7*dmats1[7][3] + coeff0_8*dmats1[8][3] + coeff0_9*dmats1[9][3];

2143

new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4] + coeff0_6*dmats1[6][4] + coeff0_7*dmats1[7][4] + coeff0_8*dmats1[8][4] + coeff0_9*dmats1[9][4];

2144

new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5] + coeff0_6*dmats1[6][5] + coeff0_7*dmats1[7][5] + coeff0_8*dmats1[8][5] + coeff0_9*dmats1[9][5];

2145

new_coeff0_6 = coeff0_0*dmats1[0][6] + coeff0_1*dmats1[1][6] + coeff0_2*dmats1[2][6] + coeff0_3*dmats1[3][6] + coeff0_4*dmats1[4][6] + coeff0_5*dmats1[5][6] + coeff0_6*dmats1[6][6] + coeff0_7*dmats1[7][6] + coeff0_8*dmats1[8][6] + coeff0_9*dmats1[9][6];

2146

new_coeff0_7 = coeff0_0*dmats1[0][7] + coeff0_1*dmats1[1][7] + coeff0_2*dmats1[2][7] + coeff0_3*dmats1[3][7] + coeff0_4*dmats1[4][7] + coeff0_5*dmats1[5][7] + coeff0_6*dmats1[6][7] + coeff0_7*dmats1[7][7] + coeff0_8*dmats1[8][7] + coeff0_9*dmats1[9][7];

2147

new_coeff0_8 = coeff0_0*dmats1[0][8] + coeff0_1*dmats1[1][8] + coeff0_2*dmats1[2][8] + coeff0_3*dmats1[3][8] + coeff0_4*dmats1[4][8] + coeff0_5*dmats1[5][8] + coeff0_6*dmats1[6][8] + coeff0_7*dmats1[7][8] + coeff0_8*dmats1[8][8] + coeff0_9*dmats1[9][8];

2148

new_coeff0_9 = coeff0_0*dmats1[0][9] + coeff0_1*dmats1[1][9] + coeff0_2*dmats1[2][9] + coeff0_3*dmats1[3][9] + coeff0_4*dmats1[4][9] + coeff0_5*dmats1[5][9] + coeff0_6*dmats1[6][9] + coeff0_7*dmats1[7][9] + coeff0_8*dmats1[8][9] + coeff0_9*dmats1[9][9];

2149

}

2150

2151

}

2152

// Compute derivatives on reference element as dot product of coefficients and basisvalues

2153

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5 + new_coeff0_6*basisvalue6 + new_coeff0_7*basisvalue7 + new_coeff0_8*basisvalue8 + new_coeff0_9*basisvalue9;

2154

}

2155

2156

// Transform derivatives back to physical element

2157

for (unsigned int row = 0; row < num_derivatives; row++)

2158

{

2159

for (unsigned int col = 0; col < num_derivatives; col++)

2160

{

2161

values[row] += transform[row][col]*derivatives[col];

2162

}

2163

}

2164

// Delete pointer to array of derivatives on FIAT element

2165

delete [] derivatives;

2166

2167

// Delete pointer to array of combinations of derivatives and transform

2168

for (unsigned int row = 0; row < num_derivatives; row++)

2169

{

2170

delete [] combinations[row];

2171

delete [] transform[row];

2172

}

2173

2174

delete [] combinations;

2175

delete [] transform;

2176

}

2177

2178

/// Evaluate order n derivatives of all basis functions at given point in cell

2179

virtual void evaluate_basis_derivatives_all(unsigned int n,

2180

double* values,

2181

const double* coordinates,

2182

const ufc::cell& c) const

2183

{

2184

throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");

2185

}

2186

2187

/// Evaluate linear functional for dof i on the function f

2188

virtual double evaluate_dof(unsigned int i,

2189

const ufc::function& f,

2190

const ufc::cell& c) const

2191

{

2192

// The reference points, direction and weights:

2193

static const double X[10][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.666666667, 0.333333333}}, {{0.333333333, 0.666666667}}, {{0, 0.333333333}}, {{0, 0.666666667}}, {{0.333333333, 0}}, {{0.666666667, 0}}, {{0.333333333, 0.333333333}}};

2194

static const double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};

2195

static const double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};

2196

2197

const double * const * x = c.coordinates;

2198

double result = 0.0;

2199

// Iterate over the points:

2200

// Evaluate basis functions for affine mapping

2201

const double w0 = 1.0 - X[i][0][0] - X[i][0][1];

2202

const double w1 = X[i][0][0];

2203

const double w2 = X[i][0][1];

2204

2205

// Compute affine mapping y = F(X)

2206

double y[2];

2207

y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];

2208

y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];

2209

2210

// Evaluate function at physical points

2211

double values[1];

2212

f.evaluate(values, y, c);

2213

2214

// Map function values using appropriate mapping

2215

// Affine map: Do nothing

2216

2217

// Note that we do not map the weights (yet).

2218

2219

// Take directional components

2220

for(int k = 0; k < 1; k++)

2221

result += values[k]*D[i][0][k];

2222

// Multiply by weights

2223

result *= W[i][0];

2224

2225

return result;

2226

}

2227

2228

/// Evaluate linear functionals for all dofs on the function f

2229

virtual void evaluate_dofs(double* values,

2230

const ufc::function& f,

2231

const ufc::cell& c) const

2232

{

2233

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2234

}

2235

2236

/// Interpolate vertex values from dof values

2237

virtual void interpolate_vertex_values(double* vertex_values,

2238

const double* dof_values,

2239

const ufc::cell& c) const

2240

{

2241

// Evaluate at vertices and use affine mapping

2242

vertex_values[0] = dof_values[0];

2243

vertex_values[1] = dof_values[1];

2244

vertex_values[2] = dof_values[2];

2245

}

2246

2247

/// Return the number of sub elements (for a mixed element)

2248

virtual unsigned int num_sub_elements() const

2249

{

2250

return 1;

2251

}

2252

2253

/// Create a new finite element for sub element i (for a mixed element)

2254

virtual ufc::finite_element* create_sub_element(unsigned int i) const

2255

{

2256

return new optimization_1_finite_element_0();

2257

}

2258

2259

};

2260

2261

/// This class defines the interface for a finite element.

2262

2263

class optimization_1_finite_element_1: public ufc::finite_element

2264

{

2265

public:

2266

2267

/// Constructor

2268

optimization_1_finite_element_1() : ufc::finite_element()

2269

{

2270

// Do nothing

2271

}

2272

2273

/// Destructor

2274

virtual ~optimization_1_finite_element_1()

2275

{

2276

// Do nothing

2277

}

2278

2279

/// Return a string identifying the finite element

2280

virtual const char* signature() const

2281

{

2282

return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3)";

2283

}

2284

2285

/// Return the cell shape

2286

virtual ufc::shape cell_shape() const

2287

{

2288

return ufc::triangle;

2289

}

2290

2291

/// Return the dimension of the finite element function space

2292

virtual unsigned int space_dimension() const

2293

{

2294

return 10;

2295

}

2296

2297

/// Return the rank of the value space

2298

virtual unsigned int value_rank() const

2299

{

2300

return 0;

2301

}

2302

2303

/// Return the dimension of the value space for axis i

2304

virtual unsigned int value_dimension(unsigned int i) const

2305

{

2306

return 1;

2307

}

2308

2309

/// Evaluate basis function i at given point in cell

2310

virtual void evaluate_basis(unsigned int i,

2311

double* values,

2312

const double* coordinates,

2313

const ufc::cell& c) const

2314

{

2315

// Extract vertex coordinates

2316

const double * const * element_coordinates = c.coordinates;

2317

2318

// Compute Jacobian of affine map from reference cell

2319

const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];

2320

const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];

2321

const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];

2322

const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];

2323

2324

// Compute determinant of Jacobian

2325

const double detJ = J_00*J_11 - J_01*J_10;

2326

2327

// Compute inverse of Jacobian

2328

2329

// Get coordinates and map to the reference (UFC) element

2330

double x = (element_coordinates[0][1]*element_coordinates[2][0] -\

2331

element_coordinates[0][0]*element_coordinates[2][1] +\

2332

J_11*coordinates[0] - J_01*coordinates[1]) / detJ;

2333

double y = (element_coordinates[1][1]*element_coordinates[0][0] -\

2334

element_coordinates[1][0]*element_coordinates[0][1] -\

2335

J_10*coordinates[0] + J_00*coordinates[1]) / detJ;

2336

2337

// Map coordinates to the reference square

2338

if (std::abs(y - 1.0) < 1e-08)

2339

x = -1.0;

2340

else

2341

x = 2.0 *x/(1.0 - y) - 1.0;

2342

y = 2.0*y - 1.0;

2343

2344

// Reset values

2345

*values = 0;

2346

2347

// Map degree of freedom to element degree of freedom

2348

const unsigned int dof = i;

2349

2350

// Generate scalings

2351

const double scalings_y_0 = 1;

2352

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

2353

const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);

2354

const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);

2355

2356

// Compute psitilde_a

2357

const double psitilde_a_0 = 1;

2358

const double psitilde_a_1 = x;

2359

const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;

2360

const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;

2361

2362

// Compute psitilde_bs

2363

const double psitilde_bs_0_0 = 1;

2364

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2365

const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;

2366

const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;

2367

const double psitilde_bs_1_0 = 1;

2368

const double psitilde_bs_1_1 = 2.5*y + 1.5;

2369

const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;

2370

const double psitilde_bs_2_0 = 1;

2371

const double psitilde_bs_2_1 = 3.5*y + 2.5;

2372

const double psitilde_bs_3_0 = 1;

2373

2374

// Compute basisvalues

2375

const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;

2376

const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;

2377

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2378

const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;

2379

const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;

2380

const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;

2381

const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;

2382

const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;

2383

const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;

2384

const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;

2385

2386

// Table(s) of coefficients

2387

static const double coefficients0[10][10] = \

2388

{{0.0471404521, -0.0288675135, -0.0166666667, 0.0782460796, 0.0606091527, 0.0349927106, -0.0601337794, -0.0508223195, -0.0393667994, -0.0227284323},

2389

{0.0471404521, 0.0288675135, -0.0166666667, 0.0782460796, -0.0606091527, 0.0349927106, 0.0601337794, -0.0508223195, 0.0393667994, -0.0227284323},

2390

{0.0471404521, 0, 0.0333333333, 0, 0, 0.104978132, 0, 0, 0, 0.090913729},

2391

{0.106066017, 0.259807621, -0.15, 0.117369119, 0.0606091527, -0.0787335989, 0, 0.101644639, -0.131222665, 0.090913729},

2392

{0.106066017, 0, 0.3, 0, 0.151522882, 0.026244533, 0, 0, 0.131222665, -0.136370594},

2393

{0.106066017, -0.259807621, -0.15, 0.117369119, -0.0606091527, -0.0787335989, 0, 0.101644639, 0.131222665, 0.090913729},

2394

{0.106066017, 0, 0.3, 0, -0.151522882, 0.026244533, 0, 0, -0.131222665, -0.136370594},

2395

{0.106066017, -0.259807621, -0.15, -0.0782460796, 0.090913729, 0.0962299542, 0.180401338, 0.0508223195, -0.0131222665, -0.0227284323},

2396

{0.106066017, 0.259807621, -0.15, -0.0782460796, -0.090913729, 0.0962299542, -0.180401338, 0.0508223195, 0.0131222665, -0.0227284323},

2397

{0.636396103, 0, 0, -0.234738239, 0, -0.26244533, 0, -0.203289278, 0, 0.090913729}};

2398

2399

// Extract relevant coefficients

2400

const double coeff0_0 = coefficients0[dof][0];

2401

const double coeff0_1 = coefficients0[dof][1];

2402

const double coeff0_2 = coefficients0[dof][2];

2403

const double coeff0_3 = coefficients0[dof][3];

2404

const double coeff0_4 = coefficients0[dof][4];

2405

const double coeff0_5 = coefficients0[dof][5];

2406

const double coeff0_6 = coefficients0[dof][6];

2407

const double coeff0_7 = coefficients0[dof][7];

2408

const double coeff0_8 = coefficients0[dof][8];

2409

const double coeff0_9 = coefficients0[dof][9];

2410

2411

// Compute value(s)

2412

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5 + coeff0_6*basisvalue6 + coeff0_7*basisvalue7 + coeff0_8*basisvalue8 + coeff0_9*basisvalue9;

2413

}

2414

2415

/// Evaluate all basis functions at given point in cell

2416

virtual void evaluate_basis_all(double* values,

2417

const double* coordinates,

2418

const ufc::cell& c) const

2419

{

2420

throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");

2421

}

2422

2423

/// Evaluate order n derivatives of basis function i at given point in cell

2424

virtual void evaluate_basis_derivatives(unsigned int i,

2425

unsigned int n,

2426

double* values,

2427

const double* coordinates,

2428

const ufc::cell& c) const

2429

{

2430

// Extract vertex coordinates

2431

const double * const * element_coordinates = c.coordinates;

2432

2433

// Compute Jacobian of affine map from reference cell

2434

const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];

2435

const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];

2436

const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];

2437

const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];

2438

2439

// Compute determinant of Jacobian

2440

const double detJ = J_00*J_11 - J_01*J_10;

2441

2442

// Compute inverse of Jacobian

2443

2444

// Get coordinates and map to the reference (UFC) element

2445

double x = (element_coordinates[0][1]*element_coordinates[2][0] -\

2446

element_coordinates[0][0]*element_coordinates[2][1] +\

2447

J_11*coordinates[0] - J_01*coordinates[1]) / detJ;

2448

double y = (element_coordinates[1][1]*element_coordinates[0][0] -\

2449

element_coordinates[1][0]*element_coordinates[0][1] -\

2450

J_10*coordinates[0] + J_00*coordinates[1]) / detJ;

2451

2452

// Map coordinates to the reference square

2453

if (std::abs(y - 1.0) < 1e-08)

2454

x = -1.0;

2455

else

2456

x = 2.0 *x/(1.0 - y) - 1.0;

2457

y = 2.0*y - 1.0;

2458

2459

// Compute number of derivatives

2460

unsigned int num_derivatives = 1;

2461

2462

for (unsigned int j = 0; j < n; j++)

2463

num_derivatives *= 2;

2464

2465

2466

// Declare pointer to two dimensional array that holds combinations of derivatives and initialise

2467

unsigned int **combinations = new unsigned int *[num_derivatives];

2468

2469

for (unsigned int j = 0; j < num_derivatives; j++)

2470

{

2471

combinations[j] = new unsigned int [n];

2472

for (unsigned int k = 0; k < n; k++)

2473

combinations[j][k] = 0;

2474

}

2475

2476

// Generate combinations of derivatives

2477

for (unsigned int row = 1; row < num_derivatives; row++)

2478

{

2479

for (unsigned int num = 0; num < row; num++)

2480

{

2481

for (unsigned int col = n-1; col+1 > 0; col--)

2482

{

2483

if (combinations[row][col] + 1 > 1)

2484

combinations[row][col] = 0;

2485

else

2486

{

2487

combinations[row][col] += 1;

2488

break;

2489

}

2490

}

2491

}

2492

}

2493

2494

// Compute inverse of Jacobian

2495

const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};

2496

2497

// Declare transformation matrix

2498

// Declare pointer to two dimensional array and initialise

2499

double **transform = new double *[num_derivatives];

2500

2501

for (unsigned int j = 0; j < num_derivatives; j++)

2502

{

2503

transform[j] = new double [num_derivatives];

2504

for (unsigned int k = 0; k < num_derivatives; k++)

2505

transform[j][k] = 1;

2506

}

2507

2508

// Construct transformation matrix

2509

for (unsigned int row = 0; row < num_derivatives; row++)

2510

{

2511

for (unsigned int col = 0; col < num_derivatives; col++)

2512

{

2513

for (unsigned int k = 0; k < n; k++)

2514

transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];

2515

}

2516

}

2517

2518

// Reset values

2519

for (unsigned int j = 0; j < 1*num_derivatives; j++)

2520

values[j] = 0;

2521

2522

// Map degree of freedom to element degree of freedom

2523

const unsigned int dof = i;

2524

2525

// Generate scalings

2526

const double scalings_y_0 = 1;

2527

const double scalings_y_1 = scalings_y_0*(0.5 - 0.5*y);

2528

const double scalings_y_2 = scalings_y_1*(0.5 - 0.5*y);

2529

const double scalings_y_3 = scalings_y_2*(0.5 - 0.5*y);

2530

2531

// Compute psitilde_a

2532

const double psitilde_a_0 = 1;

2533

const double psitilde_a_1 = x;

2534

const double psitilde_a_2 = 1.5*x*psitilde_a_1 - 0.5*psitilde_a_0;

2535

const double psitilde_a_3 = 1.66666667*x*psitilde_a_2 - 0.666666667*psitilde_a_1;

2536

2537

// Compute psitilde_bs

2538

const double psitilde_bs_0_0 = 1;

2539

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2540

const double psitilde_bs_0_2 = 0.111111111*psitilde_bs_0_1 + 1.66666667*y*psitilde_bs_0_1 - 0.555555556*psitilde_bs_0_0;

2541

const double psitilde_bs_0_3 = 0.05*psitilde_bs_0_2 + 1.75*y*psitilde_bs_0_2 - 0.7*psitilde_bs_0_1;

2542

const double psitilde_bs_1_0 = 1;

2543

const double psitilde_bs_1_1 = 2.5*y + 1.5;

2544

const double psitilde_bs_1_2 = 0.54*psitilde_bs_1_1 + 2.1*y*psitilde_bs_1_1 - 0.56*psitilde_bs_1_0;

2545

const double psitilde_bs_2_0 = 1;

2546

const double psitilde_bs_2_1 = 3.5*y + 2.5;

2547

const double psitilde_bs_3_0 = 1;

2548

2549

// Compute basisvalues

2550

const double basisvalue0 = 0.707106781*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;

2551

const double basisvalue1 = 1.73205081*psitilde_a_1*scalings_y_1*psitilde_bs_1_0;

2552

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2553

const double basisvalue3 = 2.73861279*psitilde_a_2*scalings_y_2*psitilde_bs_2_0;

2554

const double basisvalue4 = 2.12132034*psitilde_a_1*scalings_y_1*psitilde_bs_1_1;

2555

const double basisvalue5 = 1.22474487*psitilde_a_0*scalings_y_0*psitilde_bs_0_2;

2556

const double basisvalue6 = 3.74165739*psitilde_a_3*scalings_y_3*psitilde_bs_3_0;

2557

const double basisvalue7 = 3.16227766*psitilde_a_2*scalings_y_2*psitilde_bs_2_1;

2558

const double basisvalue8 = 2.44948974*psitilde_a_1*scalings_y_1*psitilde_bs_1_2;

2559

const double basisvalue9 = 1.41421356*psitilde_a_0*scalings_y_0*psitilde_bs_0_3;

2560

2561

// Table(s) of coefficients

2562

static const double coefficients0[10][10] = \

2563

{{0.0471404521, -0.0288675135, -0.0166666667, 0.0782460796, 0.0606091527, 0.0349927106, -0.0601337794, -0.0508223195, -0.0393667994, -0.0227284323},

2564

{0.0471404521, 0.0288675135, -0.0166666667, 0.0782460796, -0.0606091527, 0.0349927106, 0.0601337794, -0.0508223195, 0.0393667994, -0.0227284323},

2565

{0.0471404521, 0, 0.0333333333, 0, 0, 0.104978132, 0, 0, 0, 0.090913729},

2566

{0.106066017, 0.259807621, -0.15, 0.117369119, 0.0606091527, -0.0787335989, 0, 0.101644639, -0.131222665, 0.090913729},

2567

{0.106066017, 0, 0.3, 0, 0.151522882, 0.026244533, 0, 0, 0.131222665, -0.136370594},

2568

{0.106066017, -0.259807621, -0.15, 0.117369119, -0.0606091527, -0.0787335989, 0, 0.101644639, 0.131222665, 0.090913729},

2569

{0.106066017, 0, 0.3, 0, -0.151522882, 0.026244533, 0, 0, -0.131222665, -0.136370594},

2570

{0.106066017, -0.259807621, -0.15, -0.0782460796, 0.090913729, 0.0962299542, 0.180401338, 0.0508223195, -0.0131222665, -0.0227284323},

2571

{0.106066017, 0.259807621, -0.15, -0.0782460796, -0.090913729, 0.0962299542, -0.180401338, 0.0508223195, 0.0131222665, -0.0227284323},

2572

{0.636396103, 0, 0, -0.234738239, 0, -0.26244533, 0, -0.203289278, 0, 0.090913729}};

2573

2574

// Interesting (new) part

2575

// Tables of derivatives of the polynomial base (transpose)

2576

static const double dmats0[10][10] = \

2577

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2578

{4.89897949, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2579

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2580

{0, 9.48683298, 0, 0, 0, 0, 0, 0, 0, 0},

2581

{4, 0, 7.07106781, 0, 0, 0, 0, 0, 0, 0},

2582

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2583

{5.29150262, 0, -2.99332591, 13.662601, 0, 0.611010093, 0, 0, 0, 0},

2584

{0, 4.38178046, 0, 0, 12.5219807, 0, 0, 0, 0, 0},

2585

{3.46410162, 0, 7.83836718, 0, 0, 8.4, 0, 0, 0, 0},

2586

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};

2587

2588

static const double dmats1[10][10] = \

2589

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2590

{2.44948974, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2591

{4.24264069, 0, 0, 0, 0, 0, 0, 0, 0, 0},

2592

{2.5819889, 4.74341649, -0.912870929, 0, 0, 0, 0, 0, 0, 0},

2593

{2, 6.12372436, 3.53553391, 0, 0, 0, 0, 0, 0, 0},

2594

{-2.30940108, 0, 8.16496581, 0, 0, 0, 0, 0, 0, 0},

2595

{2.64575131, 5.18459256, -1.49666295, 6.83130051, -1.05830052, 0.305505046, 0, 0, 0, 0},

2596

{2.23606798, 2.19089023, 2.52982213, 8.08290377, 6.26099034, -1.80739223, 0, 0, 0, 0},

2597

{1.73205081, -5.09116882, 3.91918359, 0, 9.69948452, 4.2, 0, 0, 0, 0},

2598

{5, 0, -2.82842712, 0, 0, 12.1243557, 0, 0, 0, 0}};

2599

2600

// Compute reference derivatives

2601

// Declare pointer to array of derivatives on FIAT element

2602

double *derivatives = new double [num_derivatives];

2603

2604

// Declare coefficients

2605

double coeff0_0 = 0;

2606

double coeff0_1 = 0;

2607

double coeff0_2 = 0;

2608

double coeff0_3 = 0;

2609

double coeff0_4 = 0;

2610

double coeff0_5 = 0;

2611

double coeff0_6 = 0;

2612

double coeff0_7 = 0;

2613

double coeff0_8 = 0;

2614

double coeff0_9 = 0;

2615

2616

// Declare new coefficients

2617

double new_coeff0_0 = 0;

2618

double new_coeff0_1 = 0;

2619

double new_coeff0_2 = 0;

2620

double new_coeff0_3 = 0;

2621

double new_coeff0_4 = 0;

2622

double new_coeff0_5 = 0;

2623

double new_coeff0_6 = 0;

2624

double new_coeff0_7 = 0;

2625

double new_coeff0_8 = 0;

2626

double new_coeff0_9 = 0;

2627

2628

// Loop possible derivatives

2629

for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)

2630

{

2631

// Get values from coefficients array

2632

new_coeff0_0 = coefficients0[dof][0];

2633

new_coeff0_1 = coefficients0[dof][1];

2634

new_coeff0_2 = coefficients0[dof][2];

2635

new_coeff0_3 = coefficients0[dof][3];

2636

new_coeff0_4 = coefficients0[dof][4];

2637

new_coeff0_5 = coefficients0[dof][5];

2638

new_coeff0_6 = coefficients0[dof][6];

2639

new_coeff0_7 = coefficients0[dof][7];

2640

new_coeff0_8 = coefficients0[dof][8];

2641

new_coeff0_9 = coefficients0[dof][9];

2642

2643

// Loop derivative order

2644

for (unsigned int j = 0; j < n; j++)

2645

{

2646

// Update old coefficients

2647

coeff0_0 = new_coeff0_0;

2648

coeff0_1 = new_coeff0_1;

2649

coeff0_2 = new_coeff0_2;

2650

coeff0_3 = new_coeff0_3;

2651

coeff0_4 = new_coeff0_4;

2652

coeff0_5 = new_coeff0_5;

2653

coeff0_6 = new_coeff0_6;

2654

coeff0_7 = new_coeff0_7;

2655

coeff0_8 = new_coeff0_8;

2656

coeff0_9 = new_coeff0_9;

2657

2658

if(combinations[deriv_num][j] == 0)

2659

{

2660

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0] + coeff0_6*dmats0[6][0] + coeff0_7*dmats0[7][0] + coeff0_8*dmats0[8][0] + coeff0_9*dmats0[9][0];

2661

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1] + coeff0_6*dmats0[6][1] + coeff0_7*dmats0[7][1] + coeff0_8*dmats0[8][1] + coeff0_9*dmats0[9][1];

2662

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2] + coeff0_6*dmats0[6][2] + coeff0_7*dmats0[7][2] + coeff0_8*dmats0[8][2] + coeff0_9*dmats0[9][2];

2663

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3] + coeff0_6*dmats0[6][3] + coeff0_7*dmats0[7][3] + coeff0_8*dmats0[8][3] + coeff0_9*dmats0[9][3];

2664

new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4] + coeff0_6*dmats0[6][4] + coeff0_7*dmats0[7][4] + coeff0_8*dmats0[8][4] + coeff0_9*dmats0[9][4];

2665

new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5] + coeff0_6*dmats0[6][5] + coeff0_7*dmats0[7][5] + coeff0_8*dmats0[8][5] + coeff0_9*dmats0[9][5];

2666

new_coeff0_6 = coeff0_0*dmats0[0][6] + coeff0_1*dmats0[1][6] + coeff0_2*dmats0[2][6] + coeff0_3*dmats0[3][6] + coeff0_4*dmats0[4][6] + coeff0_5*dmats0[5][6] + coeff0_6*dmats0[6][6] + coeff0_7*dmats0[7][6] + coeff0_8*dmats0[8][6] + coeff0_9*dmats0[9][6];

2667

new_coeff0_7 = coeff0_0*dmats0[0][7] + coeff0_1*dmats0[1][7] + coeff0_2*dmats0[2][7] + coeff0_3*dmats0[3][7] + coeff0_4*dmats0[4][7] + coeff0_5*dmats0[5][7] + coeff0_6*dmats0[6][7] + coeff0_7*dmats0[7][7] + coeff0_8*dmats0[8][7] + coeff0_9*dmats0[9][7];

2668

new_coeff0_8 = coeff0_0*dmats0[0][8] + coeff0_1*dmats0[1][8] + coeff0_2*dmats0[2][8] + coeff0_3*dmats0[3][8] + coeff0_4*dmats0[4][8] + coeff0_5*dmats0[5][8] + coeff0_6*dmats0[6][8] + coeff0_7*dmats0[7][8] + coeff0_8*dmats0[8][8] + coeff0_9*dmats0[9][8];

2669

new_coeff0_9 = coeff0_0*dmats0[0][9] + coeff0_1*dmats0[1][9] + coeff0_2*dmats0[2][9] + coeff0_3*dmats0[3][9] + coeff0_4*dmats0[4][9] + coeff0_5*dmats0[5][9] + coeff0_6*dmats0[6][9] + coeff0_7*dmats0[7][9] + coeff0_8*dmats0[8][9] + coeff0_9*dmats0[9][9];

2670

}

2671

if(combinations[deriv_num][j] == 1)

2672

{

2673

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0] + coeff0_6*dmats1[6][0] + coeff0_7*dmats1[7][0] + coeff0_8*dmats1[8][0] + coeff0_9*dmats1[9][0];

2674

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1] + coeff0_6*dmats1[6][1] + coeff0_7*dmats1[7][1] + coeff0_8*dmats1[8][1] + coeff0_9*dmats1[9][1];

2675

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2] + coeff0_6*dmats1[6][2] + coeff0_7*dmats1[7][2] + coeff0_8*dmats1[8][2] + coeff0_9*dmats1[9][2];

2676

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3] + coeff0_6*dmats1[6][3] + coeff0_7*dmats1[7][3] + coeff0_8*dmats1[8][3] + coeff0_9*dmats1[9][3];

2677

new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4] + coeff0_6*dmats1[6][4] + coeff0_7*dmats1[7][4] + coeff0_8*dmats1[8][4] + coeff0_9*dmats1[9][4];

2678

new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5] + coeff0_6*dmats1[6][5] + coeff0_7*dmats1[7][5] + coeff0_8*dmats1[8][5] + coeff0_9*dmats1[9][5];

2679

new_coeff0_6 = coeff0_0*dmats1[0][6] + coeff0_1*dmats1[1][6] + coeff0_2*dmats1[2][6] + coeff0_3*dmats1[3][6] + coeff0_4*dmats1[4][6] + coeff0_5*dmats1[5][6] + coeff0_6*dmats1[6][6] + coeff0_7*dmats1[7][6] + coeff0_8*dmats1[8][6] + coeff0_9*dmats1[9][6];

2680

new_coeff0_7 = coeff0_0*dmats1[0][7] + coeff0_1*dmats1[1][7] + coeff0_2*dmats1[2][7] + coeff0_3*dmats1[3][7] + coeff0_4*dmats1[4][7] + coeff0_5*dmats1[5][7] + coeff0_6*dmats1[6][7] + coeff0_7*dmats1[7][7] + coeff0_8*dmats1[8][7] + coeff0_9*dmats1[9][7];

2681

new_coeff0_8 = coeff0_0*dmats1[0][8] + coeff0_1*dmats1[1][8] + coeff0_2*dmats1[2][8] + coeff0_3*dmats1[3][8] + coeff0_4*dmats1[4][8] + coeff0_5*dmats1[5][8] + coeff0_6*dmats1[6][8] + coeff0_7*dmats1[7][8] + coeff0_8*dmats1[8][8] + coeff0_9*dmats1[9][8];

2682

new_coeff0_9 = coeff0_0*dmats1[0][9] + coeff0_1*dmats1[1][9] + coeff0_2*dmats1[2][9] + coeff0_3*dmats1[3][9] + coeff0_4*dmats1[4][9] + coeff0_5*dmats1[5][9] + coeff0_6*dmats1[6][9] + coeff0_7*dmats1[7][9] + coeff0_8*dmats1[8][9] + coeff0_9*dmats1[9][9];

2683

}

2684

2685

}

2686

// Compute derivatives on reference element as dot product of coefficients and basisvalues

2687

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5 + new_coeff0_6*basisvalue6 + new_coeff0_7*basisvalue7 + new_coeff0_8*basisvalue8 + new_coeff0_9*basisvalue9;

2688

}

2689

2690

// Transform derivatives back to physical element

2691

for (unsigned int row = 0; row < num_derivatives; row++)

2692

{

2693

for (unsigned int col = 0; col < num_derivatives; col++)

2694

{

2695

values[row] += transform[row][col]*derivatives[col];

2696

}

2697

}

2698

// Delete pointer to array of derivatives on FIAT element

2699

delete [] derivatives;

2700

2701

// Delete pointer to array of combinations of derivatives and transform

2702

for (unsigned int row = 0; row < num_derivatives; row++)

2703

{

2704

delete [] combinations[row];

2705

delete [] transform[row];

2706

}

2707

2708

delete [] combinations;

2709

delete [] transform;

2710

}

2711

2712

/// Evaluate order n derivatives of all basis functions at given point in cell

2713

virtual void evaluate_basis_derivatives_all(unsigned int n,

2714

double* values,

2715

const double* coordinates,

2716

const ufc::cell& c) const

2717

{

2718

throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");

2719

}

2720

2721

/// Evaluate linear functional for dof i on the function f

2722

virtual double evaluate_dof(unsigned int i,

2723

const ufc::function& f,

2724

const ufc::cell& c) const

2725

{

2726

// The reference points, direction and weights:

2727

static const double X[10][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.666666667, 0.333333333}}, {{0.333333333, 0.666666667}}, {{0, 0.333333333}}, {{0, 0.666666667}}, {{0.333333333, 0}}, {{0.666666667, 0}}, {{0.333333333, 0.333333333}}};

2728

static const double W[10][1] = {{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}};

2729

static const double D[10][1][1] = {{{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}, {{1}}};

2730

2731

const double * const * x = c.coordinates;

2732

double result = 0.0;

2733

// Iterate over the points:

2734

// Evaluate basis functions for affine mapping

2735

const double w0 = 1.0 - X[i][0][0] - X[i][0][1];

2736

const double w1 = X[i][0][0];

2737

const double w2 = X[i][0][1];

2738

2739

// Compute affine mapping y = F(X)

2740

double y[2];

2741

y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];

2742

y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];

2743

2744

// Evaluate function at physical points

2745

double values[1];

2746

f.evaluate(values, y, c);

2747

2748

// Map function values using appropriate mapping

2749

// Affine map: Do nothing

2750

2751

// Note that we do not map the weights (yet).

2752

2753

// Take directional components

2754

for(int k = 0; k < 1; k++)

2755

result += values[k]*D[i][0][k];

2756

// Multiply by weights

2757

result *= W[i][0];

2758

2759

return result;

2760

}

2761

2762

/// Evaluate linear functionals for all dofs on the function f

2763

virtual void evaluate_dofs(double* values,

2764

const ufc::function& f,

2765

const ufc::cell& c) const

2766

{

2767

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2768

}

2769

2770

/// Interpolate vertex values from dof values

2771

virtual void interpolate_vertex_values(double* vertex_values,

2772

const double* dof_values,

2773

const ufc::cell& c) const

2774

{

2775

// Evaluate at vertices and use affine mapping

2776

vertex_values[0] = dof_values[0];

2777

vertex_values[1] = dof_values[1];

2778

vertex_values[2] = dof_values[2];

2779

}

2780

2781

/// Return the number of sub elements (for a mixed element)

2782

virtual unsigned int num_sub_elements() const

2783

{

2784

return 1;

2785

}

2786

2787

/// Create a new finite element for sub element i (for a mixed element)

2788

virtual ufc::finite_element* create_sub_element(unsigned int i) const

2789

{

2790

return new optimization_1_finite_element_1();

2791

}

2792

2793

};

2794

2795

/// This class defines the interface for a local-to-global mapping of

2796

/// degrees of freedom (dofs).

2797

2798

class optimization_1_dof_map_0: public ufc::dof_map

2799

{

2800

private:

2801

2802

unsigned int __global_dimension;

2803

2804

public:

2805

2806

/// Constructor

2807

optimization_1_dof_map_0() : ufc::dof_map()

2808

{

2809

__global_dimension = 0;

2810

}

2811

2812

/// Destructor

2813

virtual ~optimization_1_dof_map_0()

2814

{

2815

// Do nothing

2816

}

2817

2818

/// Return a string identifying the dof map

2819

virtual const char* signature() const

2820

{

2821

return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3)";

2822

}

2823

2824

/// Return true iff mesh entities of topological dimension d are needed

2825

virtual bool needs_mesh_entities(unsigned int d) const

2826

{

2827

switch ( d )

2828

{

2829

case 0:

2830

return true;

2831

break;

2832

case 1:

2833

return true;

2834

break;

2835

case 2:

2836

return true;

2837

break;

2838

}

2839

return false;

2840

}

2841

2842

/// Initialize dof map for mesh (return true iff init_cell() is needed)

2843

virtual bool init_mesh(const ufc::mesh& m)

2844

{

2845

__global_dimension = m.num_entities[0] + 2*m.num_entities[1] + m.num_entities[2];

2846

return false;

2847

}

2848

2849

/// Initialize dof map for given cell

2850

virtual void init_cell(const ufc::mesh& m,

2851

const ufc::cell& c)

2852

{

2853

// Do nothing

2854

}

2855

2856

/// Finish initialization of dof map for cells

2857

virtual void init_cell_finalize()

2858

{

2859

// Do nothing

2860

}

2861

2862

/// Return the dimension of the global finite element function space

2863

virtual unsigned int global_dimension() const

2864

{

2865

return __global_dimension;

2866

}

2867

2868

/// Return the dimension of the local finite element function space for a cell

2869

virtual unsigned int local_dimension(const ufc::cell& c) const

2870

{

2871

return 10;

2872

}

2873

2874

/// Return the maximum dimension of the local finite element function space

2875

virtual unsigned int max_local_dimension() const

2876

{

2877

return 10;

2878

}

2879

2880

// Return the geometric dimension of the coordinates this dof map provides

2881

virtual unsigned int geometric_dimension() const

2882

{

2883

return 2;

2884

}

2885

2886

/// Return the number of dofs on each cell facet

2887

virtual unsigned int num_facet_dofs() const

2888

{

2889

return 4;

2890

}

2891

2892

/// Return the number of dofs associated with each cell entity of dimension d

2893

virtual unsigned int num_entity_dofs(unsigned int d) const

2894

{

2895

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2896

}

2897

2898

/// Tabulate the local-to-global mapping of dofs on a cell

2899

virtual void tabulate_dofs(unsigned int* dofs,

2900

const ufc::mesh& m,

2901

const ufc::cell& c) const

2902

{

2903

dofs[0] = c.entity_indices[0][0];

2904

dofs[1] = c.entity_indices[0][1];

2905

dofs[2] = c.entity_indices[0][2];

2906

unsigned int offset = m.num_entities[0];

2907

dofs[3] = offset + 2*c.entity_indices[1][0];

2908

dofs[4] = offset + 2*c.entity_indices[1][0] + 1;

2909

dofs[5] = offset + 2*c.entity_indices[1][1];

2910

dofs[6] = offset + 2*c.entity_indices[1][1] + 1;

2911

dofs[7] = offset + 2*c.entity_indices[1][2];

2912

dofs[8] = offset + 2*c.entity_indices[1][2] + 1;

2913

offset = offset + 2*m.num_entities[1];

2914

dofs[9] = offset + c.entity_indices[2][0];

2915

}

2916

2917

/// Tabulate the local-to-local mapping from facet dofs to cell dofs

2918

virtual void tabulate_facet_dofs(unsigned int* dofs,

2919

unsigned int facet) const

2920

{

2921

switch ( facet )

2922

{

2923

case 0:

2924

dofs[0] = 1;

2925

dofs[1] = 2;

2926

dofs[2] = 3;

2927

dofs[3] = 4;

2928

break;

2929

case 1:

2930

dofs[0] = 0;

2931

dofs[1] = 2;

2932

dofs[2] = 5;

2933

dofs[3] = 6;

2934

break;

2935

case 2:

2936

dofs[0] = 0;

2937

dofs[1] = 1;

2938

dofs[2] = 7;

2939

dofs[3] = 8;

2940

break;

2941

}

2942

}

2943

2944

/// Tabulate the local-to-local mapping of dofs on entity (d, i)

2945

virtual void tabulate_entity_dofs(unsigned int* dofs,

2946

unsigned int d, unsigned int i) const

2947

{

2948

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

2949

}

2950

2951

/// Tabulate the coordinates of all dofs on a cell

2952

virtual void tabulate_coordinates(double** coordinates,

2953

const ufc::cell& c) const

2954

{

2955

const double * const * x = c.coordinates;

2956

coordinates[0][0] = x[0][0];

2957

coordinates[0][1] = x[0][1];

2958

coordinates[1][0] = x[1][0];

2959

coordinates[1][1] = x[1][1];

2960

coordinates[2][0] = x[2][0];

2961

coordinates[2][1] = x[2][1];

2962

coordinates[3][0] = 0.666666667*x[1][0] + 0.333333333*x[2][0];

2963

coordinates[3][1] = 0.666666667*x[1][1] + 0.333333333*x[2][1];

2964

coordinates[4][0] = 0.333333333*x[1][0] + 0.666666667*x[2][0];

2965

coordinates[4][1] = 0.333333333*x[1][1] + 0.666666667*x[2][1];

2966

coordinates[5][0] = 0.666666667*x[0][0] + 0.333333333*x[2][0];

2967

coordinates[5][1] = 0.666666667*x[0][1] + 0.333333333*x[2][1];

2968

coordinates[6][0] = 0.333333333*x[0][0] + 0.666666667*x[2][0];

2969

coordinates[6][1] = 0.333333333*x[0][1] + 0.666666667*x[2][1];

2970

coordinates[7][0] = 0.666666667*x[0][0] + 0.333333333*x[1][0];

2971

coordinates[7][1] = 0.666666667*x[0][1] + 0.333333333*x[1][1];

2972

coordinates[8][0] = 0.333333333*x[0][0] + 0.666666667*x[1][0];

2973

coordinates[8][1] = 0.333333333*x[0][1] + 0.666666667*x[1][1];

2974

coordinates[9][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];

2975

coordinates[9][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];

2976

}

2977

2978

/// Return the number of sub dof maps (for a mixed element)

2979

virtual unsigned int num_sub_dof_maps() const

2980

{

2981

return 1;

2982

}

2983

2984

/// Create a new dof_map for sub dof map i (for a mixed element)

2985

virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const

2986

{

2987

return new optimization_1_dof_map_0();

2988

}

2989

2990

};

2991

2992

/// This class defines the interface for a local-to-global mapping of

2993

/// degrees of freedom (dofs).

2994

2995

class optimization_1_dof_map_1: public ufc::dof_map

2996

{

2997

private:

2998

2999

unsigned int __global_dimension;

3000

3001

public:

3002

3003

/// Constructor

3004

optimization_1_dof_map_1() : ufc::dof_map()

3005

{

3006

__global_dimension = 0;

3007

}

3008

3009

/// Destructor

3010

virtual ~optimization_1_dof_map_1()

3011

{

3012

// Do nothing

3013

}

3014

3015

/// Return a string identifying the dof map

3016

virtual const char* signature() const

3017

{

3018

return "FFC dof map for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3)";

3019

}

3020

3021

/// Return true iff mesh entities of topological dimension d are needed

3022

virtual bool needs_mesh_entities(unsigned int d) const

3023

{

3024

switch ( d )

3025

{

3026

case 0:

3027

return true;

3028

break;

3029

case 1:

3030

return true;

3031

break;

3032

case 2:

3033

return true;

3034

break;

3035

}

3036

return false;

3037

}

3038

3039

/// Initialize dof map for mesh (return true iff init_cell() is needed)

3040

virtual bool init_mesh(const ufc::mesh& m)

3041

{

3042

__global_dimension = m.num_entities[0] + 2*m.num_entities[1] + m.num_entities[2];

3043

return false;

3044

}

3045

3046

/// Initialize dof map for given cell

3047

virtual void init_cell(const ufc::mesh& m,

3048

const ufc::cell& c)

3049

{

3050

// Do nothing

3051

}

3052

3053

/// Finish initialization of dof map for cells

3054

virtual void init_cell_finalize()

3055

{

3056

// Do nothing

3057

}

3058

3059

/// Return the dimension of the global finite element function space

3060

virtual unsigned int global_dimension() const

3061

{

3062

return __global_dimension;

3063

}

3064

3065

/// Return the dimension of the local finite element function space for a cell

3066

virtual unsigned int local_dimension(const ufc::cell& c) const

3067

{

3068

return 10;

3069

}

3070

3071

/// Return the maximum dimension of the local finite element function space

3072

virtual unsigned int max_local_dimension() const

3073

{

3074

return 10;

3075

}

3076

3077

// Return the geometric dimension of the coordinates this dof map provides

3078

virtual unsigned int geometric_dimension() const

3079

{

3080

return 2;

3081

}

3082

3083

/// Return the number of dofs on each cell facet

3084

virtual unsigned int num_facet_dofs() const

3085

{

3086

return 4;

3087

}

3088

3089

/// Return the number of dofs associated with each cell entity of dimension d

3090

virtual unsigned int num_entity_dofs(unsigned int d) const

3091

{

3092

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

3093

}

3094

3095

/// Tabulate the local-to-global mapping of dofs on a cell

3096

virtual void tabulate_dofs(unsigned int* dofs,

3097

const ufc::mesh& m,

3098

const ufc::cell& c) const

3099

{

3100

dofs[0] = c.entity_indices[0][0];

3101

dofs[1] = c.entity_indices[0][1];

3102

dofs[2] = c.entity_indices[0][2];

3103

unsigned int offset = m.num_entities[0];

3104

dofs[3] = offset + 2*c.entity_indices[1][0];

3105

dofs[4] = offset + 2*c.entity_indices[1][0] + 1;

3106

dofs[5] = offset + 2*c.entity_indices[1][1];

3107

dofs[6] = offset + 2*c.entity_indices[1][1] + 1;

3108

dofs[7] = offset + 2*c.entity_indices[1][2];

3109

dofs[8] = offset + 2*c.entity_indices[1][2] + 1;

3110

offset = offset + 2*m.num_entities[1];

3111

dofs[9] = offset + c.entity_indices[2][0];

3112

}

3113

3114

/// Tabulate the local-to-local mapping from facet dofs to cell dofs

3115

virtual void tabulate_facet_dofs(unsigned int* dofs,

3116

unsigned int facet) const

3117

{

3118

switch ( facet )

3119

{

3120

case 0:

3121

dofs[0] = 1;

3122

dofs[1] = 2;

3123

dofs[2] = 3;

3124

dofs[3] = 4;

3125

break;

3126

case 1:

3127

dofs[0] = 0;

3128

dofs[1] = 2;

3129

dofs[2] = 5;

3130

dofs[3] = 6;

3131

break;

3132

case 2:

3133

dofs[0] = 0;

3134

dofs[1] = 1;

3135

dofs[2] = 7;

3136

dofs[3] = 8;

3137

break;

3138

}

3139

}

3140

3141

/// Tabulate the local-to-local mapping of dofs on entity (d, i)

3142

virtual void tabulate_entity_dofs(unsigned int* dofs,

3143

unsigned int d, unsigned int i) const

3144

{

3145

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

3146

}

3147

3148

/// Tabulate the coordinates of all dofs on a cell

3149

virtual void tabulate_coordinates(double** coordinates,

3150

const ufc::cell& c) const

3151

{

3152

const double * const * x = c.coordinates;

3153

coordinates[0][0] = x[0][0];

3154

coordinates[0][1] = x[0][1];

3155

coordinates[1][0] = x[1][0];

3156

coordinates[1][1] = x[1][1];

3157

coordinates[2][0] = x[2][0];

3158

coordinates[2][1] = x[2][1];

3159

coordinates[3][0] = 0.666666667*x[1][0] + 0.333333333*x[2][0];

3160

coordinates[3][1] = 0.666666667*x[1][1] + 0.333333333*x[2][1];

3161

coordinates[4][0] = 0.333333333*x[1][0] + 0.666666667*x[2][0];

3162

coordinates[4][1] = 0.333333333*x[1][1] + 0.666666667*x[2][1];

3163

coordinates[5][0] = 0.666666667*x[0][0] + 0.333333333*x[2][0];

3164

coordinates[5][1] = 0.666666667*x[0][1] + 0.333333333*x[2][1];

3165

coordinates[6][0] = 0.333333333*x[0][0] + 0.666666667*x[2][0];

3166

coordinates[6][1] = 0.333333333*x[0][1] + 0.666666667*x[2][1];

3167

coordinates[7][0] = 0.666666667*x[0][0] + 0.333333333*x[1][0];

3168

coordinates[7][1] = 0.666666667*x[0][1] + 0.333333333*x[1][1];

3169

coordinates[8][0] = 0.333333333*x[0][0] + 0.666666667*x[1][0];

3170

coordinates[8][1] = 0.333333333*x[0][1] + 0.666666667*x[1][1];

3171

coordinates[9][0] = 0.333333333*x[0][0] + 0.333333333*x[1][0] + 0.333333333*x[2][0];

3172

coordinates[9][1] = 0.333333333*x[0][1] + 0.333333333*x[1][1] + 0.333333333*x[2][1];

3173

}

3174

3175

/// Return the number of sub dof maps (for a mixed element)

3176

virtual unsigned int num_sub_dof_maps() const

3177

{

3178

return 1;

3179

}

3180

3181

/// Create a new dof_map for sub dof map i (for a mixed element)

3182

virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const

3183

{

3184

return new optimization_1_dof_map_1();

3185

}

3186

3187

};

3188

3189

/// This class defines the interface for the tabulation of the cell

3190

/// tensor corresponding to the local contribution to a form from

3191

/// the integral over a cell.

3192

3193

class optimization_1_cell_integral_0_quadrature: public ufc::cell_integral

3194

{

3195

public:

3196

3197

/// Constructor

3198

optimization_1_cell_integral_0_quadrature() : ufc::cell_integral()

3199

{

3200

// Do nothing

3201

}

3202

3203

/// Destructor

3204

virtual ~optimization_1_cell_integral_0_quadrature()

3205

{

3206

// Do nothing

3207

}

3208

3209

/// Tabulate the tensor for the contribution from a local cell

3210

virtual void tabulate_tensor(double* A,

3211

const double * const * w,

3212

const ufc::cell& c) const

3213

{

3214

// Extract vertex coordinates

3215

const double * const * x = c.coordinates;

3216

3217

// Compute Jacobian of affine map from reference cell

3218

const double J_00 = x[1][0] - x[0][0];

3219

const double J_01 = x[2][0] - x[0][0];

3220

const double J_10 = x[1][1] - x[0][1];

3221

const double J_11 = x[2][1] - x[0][1];

3222

3223

// Compute determinant of Jacobian

3224

double detJ = J_00*J_11 - J_01*J_10;

3225

3226

// Compute inverse of Jacobian

3227

3228

// Set scale factor

3229

const double det = std::abs(detJ);

3230

3231

3232

// Array of quadrature weights

3233

static const double W16[16] = {0.0235683682, 0.0353880679, 0.0225840493, 0.00542322591, 0.0441850885, 0.0663442161, 0.0423397245, 0.0101672596, 0.0441850885, 0.0663442161, 0.0423397245, 0.0101672596, 0.0235683682, 0.0353880679, 0.0225840493, 0.00542322591};

3234

// Quadrature points on the UFC reference element: (0.0654669946, 0.0571041961), (0.0502101232, 0.276843014), (0.0289120842, 0.583590432), (0.00970378513, 0.860240136), (0.311164552, 0.0571041961), (0.23864866, 0.276843014), (0.137419104, 0.583590432), (0.0461220799, 0.860240136), (0.631731252, 0.0571041961), (0.484508327, 0.276843014), (0.278990463, 0.583590432), (0.0936377844, 0.860240136), (0.877428809, 0.0571041961), (0.672946863, 0.276843014), (0.387497483, 0.583590432), (0.130056079, 0.860240136)

3235

3236

// Value of basis functions at quadrature points.

3237

static const double FE0[16][10] = \

3238

{{0.452785097, 0.0474429619, 0.0432681417, -0.0135189305, -0.0139409921, 0.368034723, -0.186845726, 0.421932693, -0.207723774, 0.0885658061},

3239

{0.00645879502, 0.0394349906, 0.0274339474, -0.0531293004, -0.0106006539, 0.854147931, -0.142076465, 0.154914052, -0.129146102, 0.252562805},

3240

{-0.0263670054, 0.0252592508, -0.054598899, -0.0693419892, 0.0570043159, 0.165357063, 0.7640068, 0.00819207629, -0.0460423009, 0.176530688},

3241

{0.0638397524, 0.00928416146, 0.394831554, -0.0364705916, 0.0593783939, -0.307024415, 0.795825649, -0.00346333405, -0.00551383498, 0.0293126643},

3242

{-0.0296351119, 0.0110353632, 0.0432681417, -0.00531782109, -0.06626152, 0.145321523, -0.134525198, 0.791866619, -0.0588298934, 0.303077897},

3243

{-0.0600402894, 0.0435224405, 0.0274339474, -0.0844512385, -0.0503848963, 0.273746478, -0.102292222, 0.235979334, -0.14779975, 0.864286197},

3244

{0.0264492636, 0.0641186652, -0.054598899, -0.212107011, 0.27094145, -0.119446619, 0.550069666, -0.0281263133, -0.101399595, 0.604099392},

3245

{0.0578762154, 0.0369909804, 0.394831554, -0.153838064, 0.28222544, -0.260654105, 0.572978603, -0.0139750622, -0.0167453887, 0.100309826},

3246

{0.0110353632, -0.0296351119, 0.0432681417, 0.145321523, -0.134525198, -0.00531782109, -0.06626152, -0.0588298934, 0.791866619, 0.303077897},

3247

{0.0435224405, -0.0600402894, 0.0274339474, 0.273746478, -0.102292222, -0.0844512385, -0.0503848963, -0.14779975, 0.235979334, 0.864286197},

3248

{0.0641186652, 0.0264492636, -0.054598899, -0.119446619, 0.550069666, -0.212107011, 0.27094145, -0.101399595, -0.0281263133, 0.604099392},

3249

{0.0369909804, 0.0578762154, 0.394831554, -0.260654105, 0.572978603, -0.153838064, 0.28222544, -0.0167453887, -0.0139750622, 0.100309826},

3250

{0.0474429619, 0.452785097, 0.0432681417, 0.368034723, -0.186845726, -0.0135189305, -0.0139409921, -0.207723774, 0.421932693, 0.0885658061},

3251

{0.0394349906, 0.00645879502, 0.0274339474, 0.854147931, -0.142076465, -0.0531293004, -0.0106006539, -0.129146102, 0.154914052, 0.252562805},

3252

{0.0252592508, -0.0263670054, -0.054598899, 0.165357063, 0.7640068, -0.0693419892, 0.0570043159, -0.0460423009, 0.00819207629, 0.176530688},

3253

{0.00928416146, 0.0638397524, 0.394831554, -0.307024415, 0.795825649, -0.0364705916, 0.0593783939, -0.00551383498, -0.00346333405, 0.0293126643}};

3254

3255

3256

// Compute element tensor using UFL quadrature representation

3257

// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)

3258

// Total number of operations to compute element tensor: 960

3259

3260

// Loop quadrature points for integral

3261

// Number of operations to compute element tensor for following IP loop = 960

3262

for (unsigned int ip = 0; ip < 16; ip++)

3263

{

3264

3265

// Function declarations

3266

double F0 = 0;

3267

3268

// Total number of operations to compute function values = 20

3269

for (unsigned int r = 0; r < 10; r++)

3270

{

3271

F0 += FE0[ip][r]*w[0][r];

3272

}// end loop over 'r'

3273

3274

// Number of operations for primary indices: 40

3275

for (unsigned int j = 0; j < 10; j++)

3276

{

3277

// Number of operations to compute entry: 4

3278

A[j] += FE0[ip][j]*F0*W16[ip]*det;

3279

}// end loop over 'j'

3280

}// end loop over 'ip'

3281

}

3282

3283

};

3284

3285

/// This class defines the interface for the tabulation of the cell

3286

/// tensor corresponding to the local contribution to a form from

3287

/// the integral over a cell.

3288

3289

class optimization_1_cell_integral_0: public ufc::cell_integral

3290

{

3291

private:

3292

3293

optimization_1_cell_integral_0_quadrature integral_0_quadrature;

3294

3295

public:

3296

3297

/// Constructor

3298

optimization_1_cell_integral_0() : ufc::cell_integral()

3299

{

3300

// Do nothing

3301

}

3302

3303

/// Destructor

3304

virtual ~optimization_1_cell_integral_0()

3305

{

3306

// Do nothing

3307

}

3308

3309

/// Tabulate the tensor for the contribution from a local cell

3310

virtual void tabulate_tensor(double* A,

3311

const double * const * w,

3312

const ufc::cell& c) const

3313

{

3314

// Reset values of the element tensor block

3315

for (unsigned int j = 0; j < 10; j++)

3316

A[j] = 0;

3317

3318

// Add all contributions to element tensor

3319

integral_0_quadrature.tabulate_tensor(A, w, c);

3320

}

3321

3322

};

3323

3324

/// This class defines the interface for the assembly of the global

3325

/// tensor corresponding to a form with r + n arguments, that is, a

3326

/// mapping

3327

///

3328

/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R

3329

///

3330

/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r

3331

/// global tensor A is defined by

3332

///

3333

/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),

3334

///

3335

/// where each argument Vj represents the application to the

3336

/// sequence of basis functions of Vj and w1, w2, ..., wn are given

3337

/// fixed functions (coefficients).

3338

3339

class optimization_form_1: public ufc::form

3340

{

3341

public:

3342

3343

/// Constructor

3344

optimization_form_1() : ufc::form()

3345

{

3346

// Do nothing

3347

}

3348

3349

/// Destructor

3350

virtual ~optimization_form_1()

3351

{

3352

// Do nothing

3353

}

3354

3355

/// Return a string identifying the form

3356

virtual const char* signature() const

3357

{

3358

return "Form([Integral(Product(BasisFunction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3), 0), Function(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 3), 0)), Measure('cell', 0, None))])";

3359

}

3360

3361

/// Return the rank of the global tensor (r)

3362

virtual unsigned int rank() const

3363

{

3364

return 1;

3365

}

3366

3367

/// Return the number of coefficients (n)

3368

virtual unsigned int num_coefficients() const

3369

{

3370

return 1;

3371

}

3372

3373

/// Return the number of cell integrals

3374

virtual unsigned int num_cell_integrals() const

3375

{

3376

return 1;

3377

}

3378

3379

/// Return the number of exterior facet integrals

3380

virtual unsigned int num_exterior_facet_integrals() const

3381

{

3382

return 0;

3383

}

3384

3385

/// Return the number of interior facet integrals

3386

virtual unsigned int num_interior_facet_integrals() const

3387

{

3388

return 0;

3389

}

3390

3391

/// Create a new finite element for argument function i

3392

virtual ufc::finite_element* create_finite_element(unsigned int i) const

3393

{

3394

switch ( i )

3395

{

3396

case 0:

3397

return new optimization_1_finite_element_0();

3398

break;

3399

case 1:

3400

return new optimization_1_finite_element_1();

3401

break;

3402

}

3403

return 0;

3404

}

3405

3406

/// Create a new dof map for argument function i

3407

virtual ufc::dof_map* create_dof_map(unsigned int i) const

3408

{

3409

switch ( i )

3410

{

3411

case 0:

3412

return new optimization_1_dof_map_0();

3413

break;

3414

case 1:

3415

return new optimization_1_dof_map_1();

3416

break;

3417

}

3418

return 0;

3419

}

3420

3421

/// Create a new cell integral on sub domain i

3422

virtual ufc::cell_integral* create_cell_integral(unsigned int i) const

3423

{

3424

return new optimization_1_cell_integral_0();

3425

}

3426

3427

/// Create a new exterior facet integral on sub domain i

3428

virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const

3429

{

3430

return 0;

3431

}

3432

3433

/// Create a new interior facet integral on sub domain i

3434

virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const

3435

{

3436

return 0;

3437

}

3438

3439

};

3440

3441

#endif