~ubuntu-branches/ubuntu/maverick/dolfin/maverick

{{0.0471404520791031, -0.0288675134594812, -0.0166666666666666, 0.0782460796435951, 0.0606091526731327, 0.0349927106111883, -0.0601337794302955, -0.0508223195384204, -0.0393667994375868, -0.0227284322524248},

143

{0.0471404520791032, 0.0288675134594812, -0.0166666666666666, 0.0782460796435952, -0.0606091526731327, 0.0349927106111883, 0.0601337794302955, -0.0508223195384204, 0.0393667994375868, -0.0227284322524248},

144

{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},

145

{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},

146

{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},

147

{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},

148

{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},

149

{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},

150

{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},

151

{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};

152

153

// Extract relevant coefficients

154

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

155

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

156

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

157

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

158

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

159

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

160

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

161

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

162

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

163

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

164

165

// Compute value(s)

166

*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;

167

}

168

169

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

170

virtual void evaluate_basis_all(double* values,

171

const double* coordinates,

172

const ufc::cell& c) const

173

{

174

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

175

}

176

177

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

178

virtual void evaluate_basis_derivatives(unsigned int i,

179

unsigned int n,

180

double* values,

181

const double* coordinates,

182

const ufc::cell& c) const

183

{

184

// Extract vertex coordinates

185

const double * const * element_coordinates = c.coordinates;

186

187

// Compute Jacobian of affine map from reference cell

188

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

189

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

190

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

191

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

192

193

// Compute determinant of Jacobian

194

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

195

196

// Compute inverse of Jacobian

197

198

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

199

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

200

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

201

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

202

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

203

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

204

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

205

206

// Map coordinates to the reference square

207

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

208

x = -1.0;

209

else

210

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

211

y = 2.0*y - 1.0;

212

213

// Compute number of derivatives

214

unsigned int num_derivatives = 1;

215

216

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

217

num_derivatives *= 2;

218

219

220

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

221

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

222

223

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

224

{

225

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

226

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

227

combinations[j][k] = 0;

228

}

229

230

// Generate combinations of derivatives

231

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

232

{

233

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

234

{

235

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

236

{

237

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

238

combinations[row][col] = 0;

239

else

240

{

241

combinations[row][col] += 1;

242

break;

243

}

244

}

245

}

246

}

247

248

// Compute inverse of Jacobian

249

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

250

251

// Declare transformation matrix

252

// Declare pointer to two dimensional array and initialise

253

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

254

255

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

256

{

257

transform[j] = new double [num_derivatives];

258

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

259

transform[j][k] = 1;

260

}

261

262

// Construct transformation matrix

263

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

264

{

265

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

266

{

267

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

268

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

269

}

270

}

271

272

// Reset values

273

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

274

values[j] = 0;

275

276

// Map degree of freedom to element degree of freedom

277

const unsigned int dof = i;

278

279

// Generate scalings

280

const double scalings_y_0 = 1;

281

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

282

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

283

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

284

285

// Compute psitilde_a

286

const double psitilde_a_0 = 1;

287

const double psitilde_a_1 = x;

288

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

289

const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;

290

291

// Compute psitilde_bs

292

const double psitilde_bs_0_0 = 1;

293

const double psitilde_bs_0_1 = 1.5*y + 0.5;

294

const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;

295

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;

296

const double psitilde_bs_1_0 = 1;

297

const double psitilde_bs_1_1 = 2.5*y + 1.5;

298

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;

299

const double psitilde_bs_2_0 = 1;

300

const double psitilde_bs_2_1 = 3.5*y + 2.5;

301

const double psitilde_bs_3_0 = 1;

302

303

// Compute basisvalues

304

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

305

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

306

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

307

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

308

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

309

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

310

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

311

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

312

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

313

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

314

315

// Table(s) of coefficients

316

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

317

318

319

{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},

320

{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},

321

{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},

322

{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},

323

{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},

324

{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},

325

{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},

326

{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};

327

328

// Interesting (new) part

329

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

330

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

331

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

332

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

333

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

334

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

335

{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0},

336

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

337

{5.29150262212918, 0, -2.99332590941915, 13.6626010212795, 0, 0.611010092660779, 0, 0, 0, 0},

338

{0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0},

339

{3.46410161513775, 0, 7.83836717690617, 0, 0, 8.4, 0, 0, 0, 0},

340

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

341

342

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

343

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

344

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

345

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

346

{2.58198889747161, 4.74341649025257, -0.912870929175276, 0, 0, 0, 0, 0, 0, 0},

347

{2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0},

348

{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0},

349

{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330391, 0, 0, 0, 0},

350

{2.23606797749979, 2.19089023002067, 2.5298221281347, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0, 0, 0, 0},

351

{1.73205080756888, -5.09116882454314, 3.91918358845309, 0, 9.69948452238571, 4.2, 0, 0, 0, 0},

352

{5, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0}};

353

354

// Compute reference derivatives

355

// Declare pointer to array of derivatives on FIAT element

356

double *derivatives = new double [num_derivatives];

357

358

// Declare coefficients

359

double coeff0_0 = 0;

360

double coeff0_1 = 0;

361

double coeff0_2 = 0;

362

double coeff0_3 = 0;

363

double coeff0_4 = 0;

364

double coeff0_5 = 0;

365

double coeff0_6 = 0;

366

double coeff0_7 = 0;

367

double coeff0_8 = 0;

368

double coeff0_9 = 0;

369

370

// Declare new coefficients

371

double new_coeff0_0 = 0;

372

double new_coeff0_1 = 0;

373

double new_coeff0_2 = 0;

374

double new_coeff0_3 = 0;

375

double new_coeff0_4 = 0;

376

double new_coeff0_5 = 0;

377

double new_coeff0_6 = 0;

378

double new_coeff0_7 = 0;

379

double new_coeff0_8 = 0;

380

double new_coeff0_9 = 0;

381

382

// Loop possible derivatives

383

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

384

{

385

// Get values from coefficients array

386

new_coeff0_0 = coefficients0[dof][0];

387

new_coeff0_1 = coefficients0[dof][1];

388

new_coeff0_2 = coefficients0[dof][2];

389

new_coeff0_3 = coefficients0[dof][3];

390

new_coeff0_4 = coefficients0[dof][4];

391

new_coeff0_5 = coefficients0[dof][5];

392

new_coeff0_6 = coefficients0[dof][6];

393

new_coeff0_7 = coefficients0[dof][7];

394

new_coeff0_8 = coefficients0[dof][8];

395

new_coeff0_9 = coefficients0[dof][9];

396

397

// Loop derivative order

398

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

399

{

400

// Update old coefficients

401

coeff0_0 = new_coeff0_0;

402

coeff0_1 = new_coeff0_1;

403

coeff0_2 = new_coeff0_2;

404

coeff0_3 = new_coeff0_3;

405

coeff0_4 = new_coeff0_4;

406

coeff0_5 = new_coeff0_5;

407

coeff0_6 = new_coeff0_6;

408

coeff0_7 = new_coeff0_7;

409

coeff0_8 = new_coeff0_8;

410

coeff0_9 = new_coeff0_9;

411

412

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

413

{

414

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];

415

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];

416

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];

417

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];

418

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];

419

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];

420

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];

421

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];

422

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];

423

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];

424

}

425

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

426

{

427

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];

428

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];

429

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];

430

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];

431

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];

432

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];

433

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];

434

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];

435

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];

436

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];

437

}

438

439

}

440

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

441

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;

442

}

443

444

// Transform derivatives back to physical element

445

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

446

{

447

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

448

{

449

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

450

}

451

}

452

// Delete pointer to array of derivatives on FIAT element

453

delete [] derivatives;

454

455

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

456

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

457

{

458

delete [] combinations[row];

459

delete [] transform[row];

460

}

461

462

delete [] combinations;

463

delete [] transform;

464

}

465

466

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

467

virtual void evaluate_basis_derivatives_all(unsigned int n,

468

double* values,

469

const double* coordinates,

470

const ufc::cell& c) const

471

{

472

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

473

}

474

475

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

476

virtual double evaluate_dof(unsigned int i,

477

const ufc::function& f,

478

const ufc::cell& c) const

479

{

480

// The reference points, direction and weights:

481

const static double X[10][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}, {{0.666666666666667, 0.333333333333333}}, {{0.333333333333333, 0.666666666666667}}, {{0, 0.333333333333333}}, {{0, 0.666666666666667}}, {{0.333333333333333, 0}}, {{0.666666666666667, 0}}, {{0.333333333333333, 0.333333333333333}}};

482

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

483

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

484

485

const double * const * x = c.coordinates;

486

double result = 0.0;

487

// Iterate over the points:

488

// Evaluate basis functions for affine mapping

489

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

490

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

491

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

492

493

// Compute affine mapping y = F(X)

494

double y[2];

495

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

496

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

497

498

// Evaluate function at physical points

499

double values[1];

500

f.evaluate(values, y, c);

501

502

// Map function values using appropriate mapping

503

// Affine map: Do nothing

504

505

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

506

507

// Take directional components

508

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

509

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

510

// Multiply by weights

511

result *= W[i][0];

512

513

return result;

514

}

515

516

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

517

virtual void evaluate_dofs(double* values,

518

const ufc::function& f,

519

const ufc::cell& c) const

520

{

521

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

522

}

523

524

/// Interpolate vertex values from dof values

525

virtual void interpolate_vertex_values(double* vertex_values,

526

const double* dof_values,

527

const ufc::cell& c) const

528

{

529

// Evaluate at vertices and use affine mapping

530

vertex_values[0] = dof_values[0];

531

vertex_values[1] = dof_values[1];

532

vertex_values[2] = dof_values[2];

533

}

534

535

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

536

virtual unsigned int num_sub_elements() const

537

{

538

return 1;

539

}

540

541

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

542

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

543

{

544

return new UFC_PoissonP3BilinearForm_finite_element_0();

545

}

546

547

};

548

549

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

550

551

class UFC_PoissonP3BilinearForm_finite_element_1: public ufc::finite_element

552

{

553

public:

554

555

/// Constructor

556

UFC_PoissonP3BilinearForm_finite_element_1() : ufc::finite_element()

557

{

558

// Do nothing

559

}

560

561

/// Destructor

562

virtual ~UFC_PoissonP3BilinearForm_finite_element_1()

563

{

564

// Do nothing

565

}

566

567

/// Return a string identifying the finite element

568

virtual const char* signature() const

569

{

570

return "FiniteElement('Lagrange', 'triangle', 3)";

571

}

572

573

/// Return the cell shape

574

virtual ufc::shape cell_shape() const

575

{

576

return ufc::triangle;

577

}

578

579

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

580

virtual unsigned int space_dimension() const

581

{

582

return 10;

583

}

584

585

/// Return the rank of the value space

586

virtual unsigned int value_rank() const

587

{

588

return 0;

589

}

590

591

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

592

virtual unsigned int value_dimension(unsigned int i) const

593

{

594

return 1;

595

}

596

597

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

598

virtual void evaluate_basis(unsigned int i,

599

double* values,

600

const double* coordinates,

601

const ufc::cell& c) const

602

{

603

// Extract vertex coordinates

604

const double * const * element_coordinates = c.coordinates;

605

606

// Compute Jacobian of affine map from reference cell

607

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

608

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

609

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

610

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

611

612

// Compute determinant of Jacobian

613

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

614

615

// Compute inverse of Jacobian

616

617

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

618

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

619

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

620

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

621

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

622

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

623

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

624

625

// Map coordinates to the reference square

626

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

627

x = -1.0;

628

else

629

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

630

y = 2.0*y - 1.0;

631

632

// Reset values

633

*values = 0;

634

635

// Map degree of freedom to element degree of freedom

636

const unsigned int dof = i;

637

638

// Generate scalings

639

const double scalings_y_0 = 1;

640

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

641

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

642

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

643

644

// Compute psitilde_a

645

const double psitilde_a_0 = 1;

646

const double psitilde_a_1 = x;

647

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

648

const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;

649

650

// Compute psitilde_bs

651

const double psitilde_bs_0_0 = 1;

652

const double psitilde_bs_0_1 = 1.5*y + 0.5;

653

const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;

654

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;

655

const double psitilde_bs_1_0 = 1;

656

const double psitilde_bs_1_1 = 2.5*y + 1.5;

657

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;

658

const double psitilde_bs_2_0 = 1;

659

const double psitilde_bs_2_1 = 3.5*y + 2.5;

660

const double psitilde_bs_3_0 = 1;

661

662

// Compute basisvalues

663

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

664

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

665

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

666

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

667

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

668

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

669

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

670

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

671

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

672

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

673

674

// Table(s) of coefficients

675

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

676

677

678

{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},

679

{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},

680

{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},

681

{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},

682

{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},

683

{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},

684

{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},

685

{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};

686

687

// Extract relevant coefficients

688

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

689

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

690

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

691

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

692

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

693

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

694

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

695

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

696

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

697

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

698

699

// Compute value(s)

700

701

}

702

703

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

704

virtual void evaluate_basis_all(double* values,

705

const double* coordinates,

706

const ufc::cell& c) const

707

{

708

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

709

}

710

711

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

712

virtual void evaluate_basis_derivatives(unsigned int i,

713

unsigned int n,

714

double* values,

715

const double* coordinates,

716

const ufc::cell& c) const

717

{

718

// Extract vertex coordinates

719

const double * const * element_coordinates = c.coordinates;

720

721

// Compute Jacobian of affine map from reference cell

722

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

723

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

724

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

725

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

726

727

// Compute determinant of Jacobian

728

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

729

730

// Compute inverse of Jacobian

731

732

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

733

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

734

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

735

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

736

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

737

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

738

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

739

740

// Map coordinates to the reference square

741

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

742

x = -1.0;

743

else

744

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

745

y = 2.0*y - 1.0;

746

747

// Compute number of derivatives

748

unsigned int num_derivatives = 1;

749

750

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

751

num_derivatives *= 2;

752

753

754

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

755

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

756

757

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

758

{

759

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

760

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

761

combinations[j][k] = 0;

762

}

763

764

// Generate combinations of derivatives

765

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

766

{

767

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

768

{

769

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

770

{

771

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

772

combinations[row][col] = 0;

773

else

774

{

775

combinations[row][col] += 1;

776

break;

777

}

778

}

779

}

780

}

781

782

// Compute inverse of Jacobian

783

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

784

785

// Declare transformation matrix

786

// Declare pointer to two dimensional array and initialise

787

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

788

789

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

790

{

791

transform[j] = new double [num_derivatives];

792

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

793

transform[j][k] = 1;

794

}

795

796

// Construct transformation matrix

797

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

798

{

799

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

800

{

801

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

802

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

803

}

804

}

805

806

// Reset values

807

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

808

values[j] = 0;

809

810

// Map degree of freedom to element degree of freedom

811

const unsigned int dof = i;

812

813

// Generate scalings

814

const double scalings_y_0 = 1;

815

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

816

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

817

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

818

819

// Compute psitilde_a

820

const double psitilde_a_0 = 1;

821

const double psitilde_a_1 = x;

822

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

823

const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;

824

825

// Compute psitilde_bs

826

const double psitilde_bs_0_0 = 1;

827

const double psitilde_bs_0_1 = 1.5*y + 0.5;

828

const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;

829

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;

830

const double psitilde_bs_1_0 = 1;

831

const double psitilde_bs_1_1 = 2.5*y + 1.5;

832

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;

833

const double psitilde_bs_2_0 = 1;

834

const double psitilde_bs_2_1 = 3.5*y + 2.5;

835

const double psitilde_bs_3_0 = 1;

836

837

// Compute basisvalues

838

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

839

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

840

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

841

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

842

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

843

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

844

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

845

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

846

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

847

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

848

849

// Table(s) of coefficients

850

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

851

852

853

{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},

854

{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},

855

{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},

856

{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},

857

{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},

858

{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},

859

{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},

860

{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};

861

862

// Interesting (new) part

863

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

864

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

865

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

866

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

867

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

868

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

869

{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0},

870

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

871

{5.29150262212918, 0, -2.99332590941915, 13.6626010212795, 0, 0.611010092660779, 0, 0, 0, 0},

872

{0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0},

873

{3.46410161513775, 0, 7.83836717690617, 0, 0, 8.4, 0, 0, 0, 0},

874

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

875

876

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

877

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

878

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

879

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

880

{2.58198889747161, 4.74341649025257, -0.912870929175276, 0, 0, 0, 0, 0, 0, 0},

881

{2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0},

882

{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0},

883

{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330391, 0, 0, 0, 0},

884

{2.23606797749979, 2.19089023002067, 2.5298221281347, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0, 0, 0, 0},

885

{1.73205080756888, -5.09116882454314, 3.91918358845309, 0, 9.69948452238571, 4.2, 0, 0, 0, 0},

886

{5, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0}};

887

888

// Compute reference derivatives

889

// Declare pointer to array of derivatives on FIAT element

890

double *derivatives = new double [num_derivatives];

891

892

// Declare coefficients

893

double coeff0_0 = 0;

894

double coeff0_1 = 0;

895

double coeff0_2 = 0;

896

double coeff0_3 = 0;

897

double coeff0_4 = 0;

898

double coeff0_5 = 0;

899

double coeff0_6 = 0;

900

double coeff0_7 = 0;

901

double coeff0_8 = 0;

902

double coeff0_9 = 0;

903

904

// Declare new coefficients

905

double new_coeff0_0 = 0;

906

double new_coeff0_1 = 0;

907

double new_coeff0_2 = 0;

908

double new_coeff0_3 = 0;

909

double new_coeff0_4 = 0;

910

double new_coeff0_5 = 0;

911

double new_coeff0_6 = 0;

912

double new_coeff0_7 = 0;

913

double new_coeff0_8 = 0;

914

double new_coeff0_9 = 0;

915

916

// Loop possible derivatives

917

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

918

{

919

// Get values from coefficients array

920

new_coeff0_0 = coefficients0[dof][0];

921

new_coeff0_1 = coefficients0[dof][1];

922

new_coeff0_2 = coefficients0[dof][2];

923

new_coeff0_3 = coefficients0[dof][3];

924

new_coeff0_4 = coefficients0[dof][4];

925

new_coeff0_5 = coefficients0[dof][5];

926

new_coeff0_6 = coefficients0[dof][6];

927

new_coeff0_7 = coefficients0[dof][7];

928

new_coeff0_8 = coefficients0[dof][8];

929

new_coeff0_9 = coefficients0[dof][9];

930

931

// Loop derivative order

932

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

933

{

934

// Update old coefficients

935

coeff0_0 = new_coeff0_0;

936

coeff0_1 = new_coeff0_1;

937

coeff0_2 = new_coeff0_2;

938

coeff0_3 = new_coeff0_3;

939

coeff0_4 = new_coeff0_4;

940

coeff0_5 = new_coeff0_5;

941

coeff0_6 = new_coeff0_6;

942

coeff0_7 = new_coeff0_7;

943

coeff0_8 = new_coeff0_8;

944

coeff0_9 = new_coeff0_9;

945

946

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

947

{

948

949

950

951

952

953

954

955

956

957

958

}

959

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

960

{

961

962

963

964

965

966

967

968

969

970

971

}

972

973

}

974

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

975

976

}

977

978

// Transform derivatives back to physical element

979

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

980

{

981

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

982

{

983

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

984

}

985

}

986

// Delete pointer to array of derivatives on FIAT element

987

delete [] derivatives;

988

989

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

990

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

991

{

992

delete [] combinations[row];

993

delete [] transform[row];

994

}

995

996

delete [] combinations;

997

delete [] transform;

998

}

999

1000

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

1001

virtual void evaluate_basis_derivatives_all(unsigned int n,

1002

double* values,

1003

const double* coordinates,

1004

const ufc::cell& c) const

1005

{

1006

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

1007

}

1008

1009

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

1010

virtual double evaluate_dof(unsigned int i,

1011

const ufc::function& f,

1012

const ufc::cell& c) const

1013

{

1014

// The reference points, direction and weights:

1015

1016

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

1017

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

1018

1019

const double * const * x = c.coordinates;

1020

double result = 0.0;

1021

// Iterate over the points:

1022

// Evaluate basis functions for affine mapping

1023

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

1024

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

1025

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

1026

1027

// Compute affine mapping y = F(X)

1028

double y[2];

1029

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

1030

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

1031

1032

// Evaluate function at physical points

1033

double values[1];

1034

f.evaluate(values, y, c);

1035

1036

// Map function values using appropriate mapping

1037

// Affine map: Do nothing

1038

1039

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

1040

1041

// Take directional components

1042

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

1043

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

1044

// Multiply by weights

1045

result *= W[i][0];

1046

1047

return result;

1048

}

1049

1050

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

1051

virtual void evaluate_dofs(double* values,

1052

const ufc::function& f,

1053

const ufc::cell& c) const

1054

{

1055

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

1056

}

1057

1058

/// Interpolate vertex values from dof values

1059

virtual void interpolate_vertex_values(double* vertex_values,

1060

const double* dof_values,

1061

const ufc::cell& c) const

1062

{

1063

// Evaluate at vertices and use affine mapping

1064

vertex_values[0] = dof_values[0];

1065

vertex_values[1] = dof_values[1];

1066

vertex_values[2] = dof_values[2];

1067

}

1068

1069

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

1070

virtual unsigned int num_sub_elements() const

1071

{

1072

return 1;

1073

}

1074

1075

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

1076

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

1077

{

1078

return new UFC_PoissonP3BilinearForm_finite_element_1();

1079

}

1080

1081

};

1082

1083

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

1084

/// degrees of freedom (dofs).

1085

1086

class UFC_PoissonP3BilinearForm_dof_map_0: public ufc::dof_map

1087

{

1088

private:

1089

1090

unsigned int __global_dimension;

1091

1092

public:

1093

1094

/// Constructor

1095

UFC_PoissonP3BilinearForm_dof_map_0() : ufc::dof_map()

1096

{

1097

__global_dimension = 0;

1098

}

1099

1100

/// Destructor

1101

virtual ~UFC_PoissonP3BilinearForm_dof_map_0()

1102

{

1103

// Do nothing

1104

}

1105

1106

/// Return a string identifying the dof map

1107

virtual const char* signature() const

1108

{

1109

return "FFC dof map for FiniteElement('Lagrange', 'triangle', 3)";

1110

}

1111

1112

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

1113

virtual bool needs_mesh_entities(unsigned int d) const

1114

{

1115

switch ( d )

1116

{

1117

case 0:

1118

return true;

1119

break;

1120

case 1:

1121

return true;

1122

break;

1123

case 2:

1124

return true;

1125

break;

1126

}

1127

return false;

1128

}

1129

1130

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

1131

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

1132

{

1133

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

1134

return false;

1135

}

1136

1137

/// Initialize dof map for given cell

1138

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

1139

const ufc::cell& c)

1140

{

1141

// Do nothing

1142

}

1143

1144

/// Finish initialization of dof map for cells

1145

virtual void init_cell_finalize()

1146

{

1147

// Do nothing

1148

}

1149

1150

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

1151

virtual unsigned int global_dimension() const

1152

{

1153

return __global_dimension;

1154

}

1155

1156

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

1157

virtual unsigned int local_dimension() const

1158

{

1159

return 10;

1160

}

1161

1162

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

1163

virtual unsigned int geometric_dimension() const

1164

{

1165

return 2;

1166

}

1167

1168

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

1169

virtual unsigned int num_facet_dofs() const

1170

{

1171

return 4;

1172

}

1173

1174

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

1175

virtual unsigned int num_entity_dofs(unsigned int d) const

1176

{

1177

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

1178

}

1179

1180

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

1181

virtual void tabulate_dofs(unsigned int* dofs,

1182

const ufc::mesh& m,

1183

const ufc::cell& c) const

1184

{

1185

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

1186

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

1187

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

1188

unsigned int offset = m.num_entities[0];

1189

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

1190

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

1191

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

1192

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

1193

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

1194

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

1195

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

1196

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

1197

}

1198

1199

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

1200

virtual void tabulate_facet_dofs(unsigned int* dofs,

1201

unsigned int facet) const

1202

{

1203

switch ( facet )

1204

{

1205

case 0:

1206

dofs[0] = 1;

1207

dofs[1] = 2;

1208

dofs[2] = 3;

1209

dofs[3] = 4;

1210

break;

1211

case 1:

1212

dofs[0] = 0;

1213

dofs[1] = 2;

1214

dofs[2] = 5;

1215

dofs[3] = 6;

1216

break;

1217

case 2:

1218

dofs[0] = 0;

1219

dofs[1] = 1;

1220

dofs[2] = 7;

1221

dofs[3] = 8;

1222

break;

1223

}

1224

}

1225

1226

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

1227

virtual void tabulate_entity_dofs(unsigned int* dofs,

1228

unsigned int d, unsigned int i) const

1229

{

1230

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

1231

}

1232

1233

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

1234

virtual void tabulate_coordinates(double** coordinates,

1235

const ufc::cell& c) const

1236

{

1237

const double * const * x = c.coordinates;

1238

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

1239

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

1240

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

1241

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

1242

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

1243

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

1244

coordinates[3][0] = 0.666666666666667*x[1][0] + 0.333333333333333*x[2][0];

1245

coordinates[3][1] = 0.666666666666667*x[1][1] + 0.333333333333333*x[2][1];

1246

coordinates[4][0] = 0.333333333333333*x[1][0] + 0.666666666666667*x[2][0];

1247

coordinates[4][1] = 0.333333333333333*x[1][1] + 0.666666666666667*x[2][1];

1248

coordinates[5][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[2][0];

1249

coordinates[5][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[2][1];

1250

coordinates[6][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[2][0];

1251

coordinates[6][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[2][1];

1252

coordinates[7][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[1][0];

1253

coordinates[7][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[1][1];

1254

coordinates[8][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[1][0];

1255

coordinates[8][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[1][1];

1256

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

1257

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

1258

}

1259

1260

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

1261

virtual unsigned int num_sub_dof_maps() const

1262

{

1263

return 1;

1264

}

1265

1266

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

1267

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

1268

{

1269

return new UFC_PoissonP3BilinearForm_dof_map_0();

1270

}

1271

1272

};

1273

1274

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

1275

/// degrees of freedom (dofs).

1276

1277

class UFC_PoissonP3BilinearForm_dof_map_1: public ufc::dof_map

1278

{

1279

private:

1280

1281

unsigned int __global_dimension;

1282

1283

public:

1284

1285

/// Constructor

1286

UFC_PoissonP3BilinearForm_dof_map_1() : ufc::dof_map()

1287

{

1288

__global_dimension = 0;

1289

}

1290

1291

/// Destructor

1292

virtual ~UFC_PoissonP3BilinearForm_dof_map_1()

1293

{

1294

// Do nothing

1295

}

1296

1297

/// Return a string identifying the dof map

1298

virtual const char* signature() const

1299

{

1300

return "FFC dof map for FiniteElement('Lagrange', 'triangle', 3)";

1301

}

1302

1303

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

1304

virtual bool needs_mesh_entities(unsigned int d) const

1305

{

1306

switch ( d )

1307

{

1308

case 0:

1309

return true;

1310

break;

1311

case 1:

1312

return true;

1313

break;

1314

case 2:

1315

return true;

1316

break;

1317

}

1318

return false;

1319

}

1320

1321

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

1322

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

1323

{

1324

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

1325

return false;

1326

}

1327

1328

/// Initialize dof map for given cell

1329

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

1330

const ufc::cell& c)

1331

{

1332

// Do nothing

1333

}

1334

1335

/// Finish initialization of dof map for cells

1336

virtual void init_cell_finalize()

1337

{

1338

// Do nothing

1339

}

1340

1341

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

1342

virtual unsigned int global_dimension() const

1343

{

1344

return __global_dimension;

1345

}

1346

1347

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

1348

virtual unsigned int local_dimension() const

1349

{

1350

return 10;

1351

}

1352

1353

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

1354

virtual unsigned int geometric_dimension() const

1355

{

1356

return 2;

1357

}

1358

1359

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

1360

virtual unsigned int num_facet_dofs() const

1361

{

1362

return 4;

1363

}

1364

1365

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

1366

virtual unsigned int num_entity_dofs(unsigned int d) const

1367

{

1368

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

1369

}

1370

1371

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

1372

virtual void tabulate_dofs(unsigned int* dofs,

1373

const ufc::mesh& m,

1374

const ufc::cell& c) const

1375

{

1376

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

1377

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

1378

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

1379

unsigned int offset = m.num_entities[0];

1380

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

1381

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

1382

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

1383

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

1384

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

1385

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

1386

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

1387

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

1388

}

1389

1390

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

1391

virtual void tabulate_facet_dofs(unsigned int* dofs,

1392

unsigned int facet) const

1393

{

1394

switch ( facet )

1395

{

1396

case 0:

1397

dofs[0] = 1;

1398

dofs[1] = 2;

1399

dofs[2] = 3;

1400

dofs[3] = 4;

1401

break;

1402

case 1:

1403

dofs[0] = 0;

1404

dofs[1] = 2;

1405

dofs[2] = 5;

1406

dofs[3] = 6;

1407

break;

1408

case 2:

1409

dofs[0] = 0;

1410

dofs[1] = 1;

1411

dofs[2] = 7;

1412

dofs[3] = 8;

1413

break;

1414

}

1415

}

1416

1417

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

1418

virtual void tabulate_entity_dofs(unsigned int* dofs,

1419

unsigned int d, unsigned int i) const

1420

{

1421

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

1422

}

1423

1424

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

1425

virtual void tabulate_coordinates(double** coordinates,

1426

const ufc::cell& c) const

1427

{

1428

const double * const * x = c.coordinates;

1429

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

1430

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

1431

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

1432

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

1433

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

1434

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

1435

coordinates[3][0] = 0.666666666666667*x[1][0] + 0.333333333333333*x[2][0];

1436

coordinates[3][1] = 0.666666666666667*x[1][1] + 0.333333333333333*x[2][1];

1437

coordinates[4][0] = 0.333333333333333*x[1][0] + 0.666666666666667*x[2][0];

1438

coordinates[4][1] = 0.333333333333333*x[1][1] + 0.666666666666667*x[2][1];

1439

coordinates[5][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[2][0];

1440

coordinates[5][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[2][1];

1441

coordinates[6][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[2][0];

1442

coordinates[6][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[2][1];

1443

coordinates[7][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[1][0];

1444

coordinates[7][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[1][1];

1445

coordinates[8][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[1][0];

1446

coordinates[8][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[1][1];

1447

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

1448

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

1449

}

1450

1451

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

1452

virtual unsigned int num_sub_dof_maps() const

1453

{

1454

return 1;

1455

}

1456

1457

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

1458

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

1459

{

1460

return new UFC_PoissonP3BilinearForm_dof_map_1();

1461

}

1462

1463

};

1464

1465

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

1466

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

1467

/// the integral over a cell.

1468

1469

class UFC_PoissonP3BilinearForm_cell_integral_0: public ufc::cell_integral

1470

{

1471

public:

1472

1473

/// Constructor

1474

UFC_PoissonP3BilinearForm_cell_integral_0() : ufc::cell_integral()

1475

{

1476

// Do nothing

1477

}

1478

1479

/// Destructor

1480

virtual ~UFC_PoissonP3BilinearForm_cell_integral_0()

1481

{

1482

// Do nothing

1483

}

1484

1485

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

1486

virtual void tabulate_tensor(double* A,

1487

const double * const * w,

1488

const ufc::cell& c) const

1489

{

1490

// Extract vertex coordinates

1491

const double * const * x = c.coordinates;

1492

1493

// Compute Jacobian of affine map from reference cell

1494

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

1495

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

1496

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

1497

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

1498

1499

// Compute determinant of Jacobian

1500

double detJ = J_00*J_11 - J_01*J_10;

1501

1502

// Compute inverse of Jacobian

1503

const double Jinv_00 = J_11 / detJ;

1504

const double Jinv_01 = -J_01 / detJ;

1505

const double Jinv_10 = -J_10 / detJ;

1506

const double Jinv_11 = J_00 / detJ;

1507

1508

// Set scale factor

1509

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

1510

1511

// Number of operations to compute element tensor = 426

1512

// Compute geometry tensors

1513

// Number of operations to compute decalrations = 16

1514

const double G0_0_0 = det*(Jinv_00*Jinv_00 + Jinv_01*Jinv_01);

1515

const double G0_0_1 = det*(Jinv_00*Jinv_10 + Jinv_01*Jinv_11);

1516

const double G0_1_0 = det*(Jinv_10*Jinv_00 + Jinv_11*Jinv_01);

1517

const double G0_1_1 = det*(Jinv_10*Jinv_10 + Jinv_11*Jinv_11);

1518

1519

// Compute element tensor

1520

// Number of operations to compute tensor = 410

1521

A[0] = 0.425*G0_0_0 + 0.425*G0_0_1 + 0.425*G0_1_0 + 0.425*G0_1_1;

1522

A[1] = -0.0875*G0_0_0 - 0.0875*G0_1_0;

1523

A[2] = -0.0874999999999998*G0_0_1 - 0.0874999999999998*G0_1_1;

1524

A[3] = -0.0374999999999998*G0_0_0 - 0.0374999999999999*G0_0_1 - 0.0374999999999999*G0_1_0 - 0.0375000000000002*G0_1_1;

1525

A[4] = -0.0374999999999997*G0_0_0 - 0.0374999999999996*G0_0_1 - 0.0374999999999997*G0_1_0 - 0.0374999999999995*G0_1_1;

1526

A[5] = 0.0374999999999999*G0_0_0 - 0.674999999999999*G0_0_1 + 0.0374999999999997*G0_1_0 - 0.674999999999999*G0_1_1;

1527

A[6] = 0.0374999999999997*G0_0_0 + 0.3375*G0_0_1 + 0.0374999999999997*G0_1_0 + 0.3375*G0_1_1;

1528

A[7] = -0.674999999999999*G0_0_0 + 0.0374999999999999*G0_0_1 - 0.674999999999999*G0_1_0 + 0.0374999999999998*G0_1_1;

1529

A[8] = 0.3375*G0_0_0 + 0.0375000000000001*G0_0_1 + 0.3375*G0_1_0 + 0.0375000000000003*G0_1_1;

1530

A[9] = 0;

1531

A[10] = -0.0875*G0_0_0 - 0.0875*G0_0_1;

1532

A[11] = 0.425*G0_0_0;

1533

A[12] = 0.0874999999999999*G0_0_1;

1534

A[13] = 0.0375000000000002*G0_0_0 + 0.712499999999999*G0_0_1;

1535

A[14] = 0.0374999999999999*G0_0_0 - 0.3*G0_0_1;

1536

A[15] = -0.0374999999999997*G0_0_0;

1537

A[16] = -0.0375*G0_0_0;

1538

A[17] = 0.3375*G0_0_0 + 0.3*G0_0_1;

1539

A[18] = -0.674999999999999*G0_0_0 - 0.7125*G0_0_1;

1540

A[19] = 0;

1541

A[20] = -0.0874999999999997*G0_1_0 - 0.0874999999999998*G0_1_1;

1542

A[21] = 0.0874999999999999*G0_1_0;

1543

A[22] = 0.424999999999999*G0_1_1;

1544

A[23] = -0.3*G0_1_0 + 0.0375000000000001*G0_1_1;

1545

A[24] = 0.712499999999999*G0_1_0 + 0.0375000000000001*G0_1_1;

1546

A[25] = 0.299999999999999*G0_1_0 + 0.337499999999999*G0_1_1;

1547

A[26] = -0.712499999999999*G0_1_0 - 0.674999999999999*G0_1_1;

1548

A[27] = -0.0375*G0_1_1;

1549

A[28] = -0.0375*G0_1_1;

1550

A[29] = 0;

1551

A[30] = -0.0374999999999998*G0_0_0 - 0.0374999999999999*G0_0_1 - 0.0375*G0_1_0 - 0.0375000000000002*G0_1_1;

1552

A[31] = 0.0375000000000002*G0_0_0 + 0.712499999999999*G0_1_0;

1553

A[32] = -0.3*G0_0_1 + 0.0375000000000001*G0_1_1;

1554

A[33] = 1.6875*G0_0_0 + 0.84375*G0_0_1 + 0.84375*G0_1_0 + 1.6875*G0_1_1;

1555

A[34] = -0.3375*G0_0_0 + 0.843749999999998*G0_0_1 - 0.16875*G0_1_0 - 0.3375*G0_1_1;

1556

A[35] = 0.337500000000001*G0_0_0 + 0.16875*G0_0_1 + 0.168750000000001*G0_1_0;

1557

A[36] = 0.337499999999999*G0_0_0 + 0.168749999999999*G0_0_1 + 0.16875*G0_1_0;

1558

A[37] = 0.168750000000001*G0_0_1 + 0.16875*G0_1_0 + 0.3375*G0_1_1;

1559

A[38] = -0.84375*G0_0_1 - 0.843749999999999*G0_1_0 - 1.6875*G0_1_1;

1560

A[39] = -2.025*G0_0_0 - 1.0125*G0_0_1 - 1.0125*G0_1_0;

1561

A[40] = -0.0374999999999997*G0_0_0 - 0.0374999999999997*G0_0_1 - 0.0374999999999996*G0_1_0 - 0.0374999999999995*G0_1_1;

1562

A[41] = 0.0374999999999999*G0_0_0 - 0.3*G0_1_0;

1563

A[42] = 0.712499999999999*G0_0_1 + 0.0375000000000002*G0_1_1;

1564

A[43] = -0.3375*G0_0_0 - 0.16875*G0_0_1 + 0.843749999999998*G0_1_0 - 0.3375*G0_1_1;

1565

A[44] = 1.6875*G0_0_0 + 0.84375*G0_0_1 + 0.84375*G0_1_0 + 1.6875*G0_1_1;

1566

A[45] = 0.337499999999999*G0_0_0 + 0.168749999999999*G0_0_1 + 0.16875*G0_1_0;

1567

A[46] = -1.6875*G0_0_0 - 0.843749999999998*G0_0_1 - 0.84375*G0_1_0;

1568

A[47] = 0.16875*G0_0_1 + 0.16875*G0_1_0 + 0.3375*G0_1_1;

1569

A[48] = 0.16875*G0_0_1 + 0.168749999999999*G0_1_0 + 0.3375*G0_1_1;

1570

A[49] = -1.0125*G0_0_1 - 1.0125*G0_1_0 - 2.025*G0_1_1;

1571

A[50] = 0.0374999999999998*G0_0_0 + 0.0374999999999997*G0_0_1 - 0.674999999999999*G0_1_0 - 0.674999999999999*G0_1_1;

1572

A[51] = -0.0374999999999997*G0_0_0;

1573

A[52] = 0.299999999999999*G0_0_1 + 0.337499999999999*G0_1_1;

1574

A[53] = 0.337500000000001*G0_0_0 + 0.168750000000001*G0_0_1 + 0.16875*G0_1_0;

1575

A[54] = 0.337499999999999*G0_0_0 + 0.16875*G0_0_1 + 0.168749999999999*G0_1_0;

1576

A[55] = 1.6875*G0_0_0 + 0.843749999999999*G0_0_1 + 0.843749999999999*G0_1_0 + 1.6875*G0_1_1;

1577

A[56] = -0.337499999999999*G0_0_0 - 1.18125*G0_0_1 - 0.168749999999999*G0_1_0 - 1.35*G0_1_1;

1578

A[57] = 0.84375*G0_0_1 + 0.843749999999999*G0_1_0;

1579

A[58] = -0.168750000000001*G0_0_1 - 0.16875*G0_1_0;

1580

A[59] = -2.025*G0_0_0 - 1.0125*G0_0_1 - 1.0125*G0_1_0;

1581

A[60] = 0.0374999999999997*G0_0_0 + 0.0374999999999997*G0_0_1 + 0.3375*G0_1_0 + 0.3375*G0_1_1;

1582

A[61] = -0.0375*G0_0_0;

1583

A[62] = -0.712499999999999*G0_0_1 - 0.674999999999999*G0_1_1;

1584

A[63] = 0.337499999999999*G0_0_0 + 0.16875*G0_0_1 + 0.168749999999999*G0_1_0;

1585

A[64] = -1.6875*G0_0_0 - 0.84375*G0_0_1 - 0.843749999999998*G0_1_0;

1586

A[65] = -0.337499999999999*G0_0_0 - 0.168749999999999*G0_0_1 - 1.18125*G0_1_0 - 1.35*G0_1_1;

1587

A[66] = 1.6875*G0_0_0 + 0.843749999999998*G0_0_1 + 0.843749999999998*G0_1_0 + 1.6875*G0_1_1;

1588

A[67] = -0.16875*G0_0_1 - 0.16875*G0_1_0;

1589

A[68] = -0.168749999999999*G0_0_1 - 0.168749999999999*G0_1_0;

1590

A[69] = 1.0125*G0_0_1 + 1.0125*G0_1_0;

1591

A[70] = -0.674999999999999*G0_0_0 - 0.674999999999999*G0_0_1 + 0.0374999999999999*G0_1_0 + 0.0374999999999998*G0_1_1;

1592

A[71] = 0.3375*G0_0_0 + 0.3*G0_1_0;

1593

A[72] = -0.0375*G0_1_1;

1594

A[73] = 0.16875*G0_0_1 + 0.168750000000001*G0_1_0 + 0.3375*G0_1_1;

1595

A[74] = 0.16875*G0_0_1 + 0.16875*G0_1_0 + 0.3375*G0_1_1;

1596

A[75] = 0.843749999999999*G0_0_1 + 0.84375*G0_1_0;

1597

A[76] = -0.16875*G0_0_1 - 0.16875*G0_1_0;

1598

A[77] = 1.6875*G0_0_0 + 0.843749999999999*G0_0_1 + 0.84375*G0_1_0 + 1.6875*G0_1_1;

1599

A[78] = -1.35*G0_0_0 - 0.16875*G0_0_1 - 1.18125*G0_1_0 - 0.3375*G0_1_1;

1600

A[79] = -1.0125*G0_0_1 - 1.0125*G0_1_0 - 2.025*G0_1_1;

1601

A[80] = 0.3375*G0_0_0 + 0.3375*G0_0_1 + 0.0375000000000001*G0_1_0 + 0.0375000000000003*G0_1_1;

1602

A[81] = -0.674999999999999*G0_0_0 - 0.712499999999999*G0_1_0;

1603

A[82] = -0.0375*G0_1_1;

1604

A[83] = -0.843749999999999*G0_0_1 - 0.84375*G0_1_0 - 1.6875*G0_1_1;

1605

A[84] = 0.168749999999999*G0_0_1 + 0.16875*G0_1_0 + 0.3375*G0_1_1;

1606

A[85] = -0.16875*G0_0_1 - 0.168750000000001*G0_1_0;

1607

A[86] = -0.168749999999999*G0_0_1 - 0.168749999999999*G0_1_0;

1608

A[87] = -1.35*G0_0_0 - 1.18125*G0_0_1 - 0.16875*G0_1_0 - 0.3375*G0_1_1;

1609

A[88] = 1.6875*G0_0_0 + 0.84375*G0_0_1 + 0.84375*G0_1_0 + 1.6875*G0_1_1;

1610

A[89] = 1.0125*G0_0_1 + 1.0125*G0_1_0;

1611

A[90] = 0;

1612

A[91] = 0;

1613

A[92] = 0;

1614

A[93] = -2.025*G0_0_0 - 1.0125*G0_0_1 - 1.0125*G0_1_0;

1615

A[94] = -1.0125*G0_0_1 - 1.0125*G0_1_0 - 2.025*G0_1_1;

1616

A[95] = -2.025*G0_0_0 - 1.0125*G0_0_1 - 1.0125*G0_1_0;

1617

A[96] = 1.0125*G0_0_1 + 1.0125*G0_1_0;

1618

A[97] = -1.0125*G0_0_1 - 1.0125*G0_1_0 - 2.025*G0_1_1;

1619

A[98] = 1.0125*G0_0_1 + 1.0125*G0_1_0;

1620

A[99] = 4.05*G0_0_0 + 2.025*G0_0_1 + 2.025*G0_1_0 + 4.05*G0_1_1;

1621

}

1622

1623

};

1624

1625

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

1626

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

1627

/// mapping

1628

///

1629

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

1630

///

1631

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

1632

/// global tensor A is defined by

1633

///

1634

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

1635

///

1636

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

1637

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

1638

/// fixed functions (coefficients).

1639

1640

class UFC_PoissonP3BilinearForm: public ufc::form

1641

{

1642

public:

1643

1644

/// Constructor

1645

UFC_PoissonP3BilinearForm() : ufc::form()

1646

{

1647

// Do nothing

1648

}

1649

1650

/// Destructor

1651

virtual ~UFC_PoissonP3BilinearForm()

1652

{

1653

// Do nothing

1654

}

1655

1656

/// Return a string identifying the form

1657

virtual const char* signature() const

1658

{

1659

return "(dXa0[0, 1]/dxb0[0, 1])(dXa1[0, 1]/dxb0[0, 1]) | ((d/dXa0[0, 1])vi0[0, 1, 2, 3, 4, 5, 6, 7, 8, 9])*((d/dXa1[0, 1])vi1[0, 1, 2, 3, 4, 5, 6, 7, 8, 9])*dX(0)";

1660

}

1661

1662

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

1663

virtual unsigned int rank() const

1664

{

1665

return 2;

1666

}

1667

1668

/// Return the number of coefficients (n)

1669

virtual unsigned int num_coefficients() const

1670

{

1671

return 0;

1672

}

1673

1674

/// Return the number of cell integrals

1675

virtual unsigned int num_cell_integrals() const

1676

{

1677

return 1;

1678

}

1679

1680

/// Return the number of exterior facet integrals

1681

virtual unsigned int num_exterior_facet_integrals() const

1682

{

1683

return 0;

1684

}

1685

1686

/// Return the number of interior facet integrals

1687

virtual unsigned int num_interior_facet_integrals() const

1688

{

1689

return 0;

1690

}

1691

1692

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

1693

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

1694

{

1695

switch ( i )

1696

{

1697

case 0:

1698

return new UFC_PoissonP3BilinearForm_finite_element_0();

1699

break;

1700

case 1:

1701

return new UFC_PoissonP3BilinearForm_finite_element_1();

1702

break;

1703

}

1704

return 0;

1705

}

1706

1707

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

1708

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

1709

{

1710

switch ( i )

1711

{

1712

case 0:

1713

return new UFC_PoissonP3BilinearForm_dof_map_0();

1714

break;

1715

case 1:

1716

return new UFC_PoissonP3BilinearForm_dof_map_1();

1717

break;

1718

}

1719

return 0;

1720

}

1721

1722

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

1723

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

1724

{

1725

return new UFC_PoissonP3BilinearForm_cell_integral_0();

1726

}

1727

1728

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

1729

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

1730

{

1731

return 0;

1732

}

1733

1734

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

1735

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

1736

{

1737

return 0;

1738

}

1739

1740

};

1741

1742

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

1743

1744

class UFC_PoissonP3LinearForm_finite_element_0: public ufc::finite_element

1745

{

1746

public:

1747

1748

/// Constructor

1749

UFC_PoissonP3LinearForm_finite_element_0() : ufc::finite_element()

1750

{

1751

// Do nothing

1752

}

1753

1754

/// Destructor

1755

virtual ~UFC_PoissonP3LinearForm_finite_element_0()

1756

{

1757

// Do nothing

1758

}

1759

1760

/// Return a string identifying the finite element

1761

virtual const char* signature() const

1762

{

1763

return "FiniteElement('Lagrange', 'triangle', 3)";

1764

}

1765

1766

/// Return the cell shape

1767

virtual ufc::shape cell_shape() const

1768

{

1769

return ufc::triangle;

1770

}

1771

1772

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

1773

virtual unsigned int space_dimension() const

1774

{

1775

return 10;

1776

}

1777

1778

/// Return the rank of the value space

1779

virtual unsigned int value_rank() const

1780

{

1781

return 0;

1782

}

1783

1784

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

1785

virtual unsigned int value_dimension(unsigned int i) const

1786

{

1787

return 1;

1788

}

1789

1790

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

1791

virtual void evaluate_basis(unsigned int i,

1792

double* values,

1793

const double* coordinates,

1794

const ufc::cell& c) const

1795

{

1796

// Extract vertex coordinates

1797

const double * const * element_coordinates = c.coordinates;

1798

1799

// Compute Jacobian of affine map from reference cell

1800

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

1801

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

1802

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

1803

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

1804

1805

// Compute determinant of Jacobian

1806

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

1807

1808

// Compute inverse of Jacobian

1809

1810

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

1811

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

1812

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

1813

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

1814

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

1815

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

1816

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

1817

1818

// Map coordinates to the reference square

1819

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

1820

x = -1.0;

1821

else

1822

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

1823

y = 2.0*y - 1.0;

1824

1825

// Reset values

1826

*values = 0;

1827

1828

// Map degree of freedom to element degree of freedom

1829

const unsigned int dof = i;

1830

1831

// Generate scalings

1832

const double scalings_y_0 = 1;

1833

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

1834

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

1835

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

1836

1837

// Compute psitilde_a

1838

const double psitilde_a_0 = 1;

1839

const double psitilde_a_1 = x;

1840

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

1841

const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;

1842

1843

// Compute psitilde_bs

1844

const double psitilde_bs_0_0 = 1;

1845

const double psitilde_bs_0_1 = 1.5*y + 0.5;

1846

const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;

1847

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;

1848

const double psitilde_bs_1_0 = 1;

1849

const double psitilde_bs_1_1 = 2.5*y + 1.5;

1850

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;

1851

const double psitilde_bs_2_0 = 1;

1852

const double psitilde_bs_2_1 = 3.5*y + 2.5;

1853

const double psitilde_bs_3_0 = 1;

1854

1855

// Compute basisvalues

1856

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

1857

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

1858

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

1859

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

1860

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

1861

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

1862

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

1863

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

1864

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

1865

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

1866

1867

// Table(s) of coefficients

1868

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

1869

1870

1871

{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},

1872

{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},

1873

{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},

1874

{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},

1875

{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},

1876

{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},

1877

{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},

1878

{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};

1879

1880

// Extract relevant coefficients

1881

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

1882

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

1883

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

1884

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

1885

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

1886

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

1887

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

1888

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

1889

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

1890

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

1891

1892

// Compute value(s)

1893

1894

}

1895

1896

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

1897

virtual void evaluate_basis_all(double* values,

1898

const double* coordinates,

1899

const ufc::cell& c) const

1900

{

1901

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

1902

}

1903

1904

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

1905

virtual void evaluate_basis_derivatives(unsigned int i,

1906

unsigned int n,

1907

double* values,

1908

const double* coordinates,

1909

const ufc::cell& c) const

1910

{

1911

// Extract vertex coordinates

1912

const double * const * element_coordinates = c.coordinates;

1913

1914

// Compute Jacobian of affine map from reference cell

1915

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

1916

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

1917

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

1918

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

1919

1920

// Compute determinant of Jacobian

1921

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

1922

1923

// Compute inverse of Jacobian

1924

1925

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

1926

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

1927

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

1928

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

1929

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

1930

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

1931

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

1932

1933

// Map coordinates to the reference square

1934

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

1935

x = -1.0;

1936

else

1937

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

1938

y = 2.0*y - 1.0;

1939

1940

// Compute number of derivatives

1941

unsigned int num_derivatives = 1;

1942

1943

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

1944

num_derivatives *= 2;

1945

1946

1947

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

1948

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

1949

1950

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

1951

{

1952

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

1953

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

1954

combinations[j][k] = 0;

1955

}

1956

1957

// Generate combinations of derivatives

1958

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

1959

{

1960

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

1961

{

1962

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

1963

{

1964

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

1965

combinations[row][col] = 0;

1966

else

1967

{

1968

combinations[row][col] += 1;

1969

break;

1970

}

1971

}

1972

}

1973

}

1974

1975

// Compute inverse of Jacobian

1976

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

1977

1978

// Declare transformation matrix

1979

// Declare pointer to two dimensional array and initialise

1980

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

1981

1982

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

1983

{

1984

transform[j] = new double [num_derivatives];

1985

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

1986

transform[j][k] = 1;

1987

}

1988

1989

// Construct transformation matrix

1990

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

1991

{

1992

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

1993

{

1994

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

1995

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

1996

}

1997

}

1998

1999

// Reset values

2000

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

2001

values[j] = 0;

2002

2003

// Map degree of freedom to element degree of freedom

2004

const unsigned int dof = i;

2005

2006

// Generate scalings

2007

const double scalings_y_0 = 1;

2008

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

2009

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

2010

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

2011

2012

// Compute psitilde_a

2013

const double psitilde_a_0 = 1;

2014

const double psitilde_a_1 = x;

2015

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

2016

const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;

2017

2018

// Compute psitilde_bs

2019

const double psitilde_bs_0_0 = 1;

2020

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2021

const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;

2022

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;

2023

const double psitilde_bs_1_0 = 1;

2024

const double psitilde_bs_1_1 = 2.5*y + 1.5;

2025

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;

2026

const double psitilde_bs_2_0 = 1;

2027

const double psitilde_bs_2_1 = 3.5*y + 2.5;

2028

const double psitilde_bs_3_0 = 1;

2029

2030

// Compute basisvalues

2031

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

2032

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

2033

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2034

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

2035

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

2036

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

2037

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

2038

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

2039

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

2040

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

2041

2042

// Table(s) of coefficients

2043

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

2044

2045

2046

{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},

2047

{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},

2048

{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},

2049

{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},

2050

{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},

2051

{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},

2052

{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},

2053

{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};

2054

2055

// Interesting (new) part

2056

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

2057

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

2058

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

2059

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

2060

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

2061

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

2062

{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0},

2063

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

2064

{5.29150262212918, 0, -2.99332590941915, 13.6626010212795, 0, 0.611010092660779, 0, 0, 0, 0},

2065

{0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0},

2066

{3.46410161513775, 0, 7.83836717690617, 0, 0, 8.4, 0, 0, 0, 0},

2067

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

2068

2069

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

2070

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

2071

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

2072

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

2073

{2.58198889747161, 4.74341649025257, -0.912870929175276, 0, 0, 0, 0, 0, 0, 0},

2074

{2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0},

2075

{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0},

2076

{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330391, 0, 0, 0, 0},

2077

{2.23606797749979, 2.19089023002067, 2.5298221281347, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0, 0, 0, 0},

2078

{1.73205080756888, -5.09116882454314, 3.91918358845309, 0, 9.69948452238571, 4.2, 0, 0, 0, 0},

2079

{5, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0}};

2080

2081

// Compute reference derivatives

2082

// Declare pointer to array of derivatives on FIAT element

2083

double *derivatives = new double [num_derivatives];

2084

2085

// Declare coefficients

2086

double coeff0_0 = 0;

2087

double coeff0_1 = 0;

2088

double coeff0_2 = 0;

2089

double coeff0_3 = 0;

2090

double coeff0_4 = 0;

2091

double coeff0_5 = 0;

2092

double coeff0_6 = 0;

2093

double coeff0_7 = 0;

2094

double coeff0_8 = 0;

2095

double coeff0_9 = 0;

2096

2097

// Declare new coefficients

2098

double new_coeff0_0 = 0;

2099

double new_coeff0_1 = 0;

2100

double new_coeff0_2 = 0;

2101

double new_coeff0_3 = 0;

2102

double new_coeff0_4 = 0;

2103

double new_coeff0_5 = 0;

2104

double new_coeff0_6 = 0;

2105

double new_coeff0_7 = 0;

2106

double new_coeff0_8 = 0;

2107

double new_coeff0_9 = 0;

2108

2109

// Loop possible derivatives

2110

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

2111

{

2112

// Get values from coefficients array

2113

new_coeff0_0 = coefficients0[dof][0];

2114

new_coeff0_1 = coefficients0[dof][1];

2115

new_coeff0_2 = coefficients0[dof][2];

2116

new_coeff0_3 = coefficients0[dof][3];

2117

new_coeff0_4 = coefficients0[dof][4];

2118

new_coeff0_5 = coefficients0[dof][5];

2119

new_coeff0_6 = coefficients0[dof][6];

2120

new_coeff0_7 = coefficients0[dof][7];

2121

new_coeff0_8 = coefficients0[dof][8];

2122

new_coeff0_9 = coefficients0[dof][9];

2123

2124

// Loop derivative order

2125

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

2126

{

2127

// Update old coefficients

2128

coeff0_0 = new_coeff0_0;

2129

coeff0_1 = new_coeff0_1;

2130

coeff0_2 = new_coeff0_2;

2131

coeff0_3 = new_coeff0_3;

2132

coeff0_4 = new_coeff0_4;

2133

coeff0_5 = new_coeff0_5;

2134

coeff0_6 = new_coeff0_6;

2135

coeff0_7 = new_coeff0_7;

2136

coeff0_8 = new_coeff0_8;

2137

coeff0_9 = new_coeff0_9;

2138

2139

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

2140

{

2141

2142

2143

2144

2145

2146

2147

2148

2149

2150

2151

}

2152

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

2153

{

2154

2155

2156

2157

2158

2159

2160

2161

2162

2163

2164

}

2165

2166

}

2167

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

2168

2169

}

2170

2171

// Transform derivatives back to physical element

2172

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

2173

{

2174

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

2175

{

2176

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

2177

}

2178

}

2179

// Delete pointer to array of derivatives on FIAT element

2180

delete [] derivatives;

2181

2182

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

2183

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

2184

{

2185

delete [] combinations[row];

2186

delete [] transform[row];

2187

}

2188

2189

delete [] combinations;

2190

delete [] transform;

2191

}

2192

2193

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

2194

virtual void evaluate_basis_derivatives_all(unsigned int n,

2195

double* values,

2196

const double* coordinates,

2197

const ufc::cell& c) const

2198

{

2199

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

2200

}

2201

2202

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

2203

virtual double evaluate_dof(unsigned int i,

2204

const ufc::function& f,

2205

const ufc::cell& c) const

2206

{

2207

// The reference points, direction and weights:

2208

2209

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

2210

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

2211

2212

const double * const * x = c.coordinates;

2213

double result = 0.0;

2214

// Iterate over the points:

2215

// Evaluate basis functions for affine mapping

2216

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

2217

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

2218

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

2219

2220

// Compute affine mapping y = F(X)

2221

double y[2];

2222

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

2223

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

2224

2225

// Evaluate function at physical points

2226

double values[1];

2227

f.evaluate(values, y, c);

2228

2229

// Map function values using appropriate mapping

2230

// Affine map: Do nothing

2231

2232

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

2233

2234

// Take directional components

2235

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

2236

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

2237

// Multiply by weights

2238

result *= W[i][0];

2239

2240

return result;

2241

}

2242

2243

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

2244

virtual void evaluate_dofs(double* values,

2245

const ufc::function& f,

2246

const ufc::cell& c) const

2247

{

2248

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

2249

}

2250

2251

/// Interpolate vertex values from dof values

2252

virtual void interpolate_vertex_values(double* vertex_values,

2253

const double* dof_values,

2254

const ufc::cell& c) const

2255

{

2256

// Evaluate at vertices and use affine mapping

2257

vertex_values[0] = dof_values[0];

2258

vertex_values[1] = dof_values[1];

2259

vertex_values[2] = dof_values[2];

2260

}

2261

2262

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

2263

virtual unsigned int num_sub_elements() const

2264

{

2265

return 1;

2266

}

2267

2268

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

2269

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

2270

{

2271

return new UFC_PoissonP3LinearForm_finite_element_0();

2272

}

2273

2274

};

2275

2276

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

2277

2278

class UFC_PoissonP3LinearForm_finite_element_1: public ufc::finite_element

2279

{

2280

public:

2281

2282

/// Constructor

2283

UFC_PoissonP3LinearForm_finite_element_1() : ufc::finite_element()

2284

{

2285

// Do nothing

2286

}

2287

2288

/// Destructor

2289

virtual ~UFC_PoissonP3LinearForm_finite_element_1()

2290

{

2291

// Do nothing

2292

}

2293

2294

/// Return a string identifying the finite element

2295

virtual const char* signature() const

2296

{

2297

return "FiniteElement('Lagrange', 'triangle', 3)";

2298

}

2299

2300

/// Return the cell shape

2301

virtual ufc::shape cell_shape() const

2302

{

2303

return ufc::triangle;

2304

}

2305

2306

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

2307

virtual unsigned int space_dimension() const

2308

{

2309

return 10;

2310

}

2311

2312

/// Return the rank of the value space

2313

virtual unsigned int value_rank() const

2314

{

2315

return 0;

2316

}

2317

2318

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

2319

virtual unsigned int value_dimension(unsigned int i) const

2320

{

2321

return 1;

2322

}

2323

2324

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

2325

virtual void evaluate_basis(unsigned int i,

2326

double* values,

2327

const double* coordinates,

2328

const ufc::cell& c) const

2329

{

2330

// Extract vertex coordinates

2331

const double * const * element_coordinates = c.coordinates;

2332

2333

// Compute Jacobian of affine map from reference cell

2334

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

2335

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

2336

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

2337

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

2338

2339

// Compute determinant of Jacobian

2340

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

2341

2342

// Compute inverse of Jacobian

2343

2344

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

2345

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

2346

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

2347

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

2348

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

2349

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

2350

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

2351

2352

// Map coordinates to the reference square

2353

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

2354

x = -1.0;

2355

else

2356

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

2357

y = 2.0*y - 1.0;

2358

2359

// Reset values

2360

*values = 0;

2361

2362

// Map degree of freedom to element degree of freedom

2363

const unsigned int dof = i;

2364

2365

// Generate scalings

2366

const double scalings_y_0 = 1;

2367

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

2368

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

2369

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

2370

2371

// Compute psitilde_a

2372

const double psitilde_a_0 = 1;

2373

const double psitilde_a_1 = x;

2374

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

2375

const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;

2376

2377

// Compute psitilde_bs

2378

const double psitilde_bs_0_0 = 1;

2379

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2380

const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;

2381

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;

2382

const double psitilde_bs_1_0 = 1;

2383

const double psitilde_bs_1_1 = 2.5*y + 1.5;

2384

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;

2385

const double psitilde_bs_2_0 = 1;

2386

const double psitilde_bs_2_1 = 3.5*y + 2.5;

2387

const double psitilde_bs_3_0 = 1;

2388

2389

// Compute basisvalues

2390

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

2391

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

2392

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2393

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

2394

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

2395

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

2396

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

2397

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

2398

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

2399

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

2400

2401

// Table(s) of coefficients

2402

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

2403

2404

2405

{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},

2406

{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},

2407

{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},

2408

{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},

2409

{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},

2410

{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},

2411

{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},

2412

{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};

2413

2414

// Extract relevant coefficients

2415

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

2416

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

2417

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

2418

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

2419

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

2420

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

2421

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

2422

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

2423

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

2424

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

2425

2426

// Compute value(s)

2427

2428

}

2429

2430

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

2431

virtual void evaluate_basis_all(double* values,

2432

const double* coordinates,

2433

const ufc::cell& c) const

2434

{

2435

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

2436

}

2437

2438

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

2439

virtual void evaluate_basis_derivatives(unsigned int i,

2440

unsigned int n,

2441

double* values,

2442

const double* coordinates,

2443

const ufc::cell& c) const

2444

{

2445

// Extract vertex coordinates

2446

const double * const * element_coordinates = c.coordinates;

2447

2448

// Compute Jacobian of affine map from reference cell

2449

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

2450

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

2451

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

2452

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

2453

2454

// Compute determinant of Jacobian

2455

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

2456

2457

// Compute inverse of Jacobian

2458

2459

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

2460

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

2461

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

2462

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

2463

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

2464

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

2465

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

2466

2467

// Map coordinates to the reference square

2468

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

2469

x = -1.0;

2470

else

2471

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

2472

y = 2.0*y - 1.0;

2473

2474

// Compute number of derivatives

2475

unsigned int num_derivatives = 1;

2476

2477

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

2478

num_derivatives *= 2;

2479

2480

2481

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

2482

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

2483

2484

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

2485

{

2486

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

2487

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

2488

combinations[j][k] = 0;

2489

}

2490

2491

// Generate combinations of derivatives

2492

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

2493

{

2494

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

2495

{

2496

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

2497

{

2498

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

2499

combinations[row][col] = 0;

2500

else

2501

{

2502

combinations[row][col] += 1;

2503

break;

2504

}

2505

}

2506

}

2507

}

2508

2509

// Compute inverse of Jacobian

2510

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

2511

2512

// Declare transformation matrix

2513

// Declare pointer to two dimensional array and initialise

2514

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

2515

2516

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

2517

{

2518

transform[j] = new double [num_derivatives];

2519

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

2520

transform[j][k] = 1;

2521

}

2522

2523

// Construct transformation matrix

2524

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

2525

{

2526

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

2527

{

2528

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

2529

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

2530

}

2531

}

2532

2533

// Reset values

2534

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

2535

values[j] = 0;

2536

2537

// Map degree of freedom to element degree of freedom

2538

const unsigned int dof = i;

2539

2540

// Generate scalings

2541

const double scalings_y_0 = 1;

2542

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

2543

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

2544

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

2545

2546

// Compute psitilde_a

2547

const double psitilde_a_0 = 1;

2548

const double psitilde_a_1 = x;

2549

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

2550

const double psitilde_a_3 = 1.66666666666667*x*psitilde_a_2 - 0.666666666666667*psitilde_a_1;

2551

2552

// Compute psitilde_bs

2553

const double psitilde_bs_0_0 = 1;

2554

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2555

const double psitilde_bs_0_2 = 0.111111111111111*psitilde_bs_0_1 + 1.66666666666667*y*psitilde_bs_0_1 - 0.555555555555556*psitilde_bs_0_0;

2556

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;

2557

const double psitilde_bs_1_0 = 1;

2558

const double psitilde_bs_1_1 = 2.5*y + 1.5;

2559

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;

2560

const double psitilde_bs_2_0 = 1;

2561

const double psitilde_bs_2_1 = 3.5*y + 2.5;

2562

const double psitilde_bs_3_0 = 1;

2563

2564

// Compute basisvalues

2565

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

2566

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

2567

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2568

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

2569

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

2570

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

2571

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

2572

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

2573

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

2574

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

2575

2576

// Table(s) of coefficients

2577

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

2578

2579

2580

{0.0471404520791031, 0, 0.0333333333333334, 0, 0, 0.104978131833565, 0, 0, 0, 0.0909137290096989},

2581

{0.106066017177982, 0.259807621135332, -0.15, 0.117369119465393, 0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, -0.131222664791956, 0.090913729009699},

2582

{0.106066017177982, 0, 0.3, 0, 0.151522881682832, 0.0262445329583912, 0, 0, 0.131222664791956, -0.136370593514548},

2583

{0.106066017177982, -0.259807621135332, -0.15, 0.117369119465393, -0.0606091526731326, -0.0787335988751736, 0, 0.101644639076841, 0.131222664791956, 0.090913729009699},

2584

{0.106066017177982, 0, 0.3, 0, -0.151522881682832, 0.0262445329583912, 0, 0, -0.131222664791956, -0.136370593514548},

2585

{0.106066017177982, -0.259807621135332, -0.15, -0.0782460796435951, 0.090913729009699, 0.0962299541807677, 0.180401338290886, 0.0508223195384204, -0.0131222664791956, -0.0227284322524248},

2586

{0.106066017177982, 0.259807621135332, -0.15, -0.0782460796435952, -0.090913729009699, 0.0962299541807677, -0.180401338290886, 0.0508223195384203, 0.0131222664791956, -0.0227284322524247},

2587

{0.636396103067893, 0, 0, -0.234738238930785, 0, -0.262445329583912, 0, -0.203289278153681, 0, 0.090913729009699}};

2588

2589

// Interesting (new) part

2590

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

2591

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

2592

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

2593

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

2594

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

2595

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

2596

{4, 0, 7.07106781186548, 0, 0, 0, 0, 0, 0, 0},

2597

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

2598

{5.29150262212918, 0, -2.99332590941915, 13.6626010212795, 0, 0.611010092660779, 0, 0, 0, 0},

2599

{0, 4.38178046004133, 0, 0, 12.5219806739988, 0, 0, 0, 0, 0},

2600

{3.46410161513775, 0, 7.83836717690617, 0, 0, 8.4, 0, 0, 0, 0},

2601

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

2602

2603

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

2604

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

2605

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

2606

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

2607

{2.58198889747161, 4.74341649025257, -0.912870929175276, 0, 0, 0, 0, 0, 0, 0},

2608

{2, 6.12372435695795, 3.53553390593274, 0, 0, 0, 0, 0, 0, 0},

2609

{-2.3094010767585, 0, 8.16496580927726, 0, 0, 0, 0, 0, 0, 0},

2610

{2.64575131106459, 5.18459255872629, -1.49666295470958, 6.83130051063973, -1.05830052442584, 0.305505046330391, 0, 0, 0, 0},

2611

{2.23606797749979, 2.19089023002067, 2.5298221281347, 8.08290376865476, 6.26099033699941, -1.80739222823013, 0, 0, 0, 0},

2612

{1.73205080756888, -5.09116882454314, 3.91918358845309, 0, 9.69948452238571, 4.2, 0, 0, 0, 0},

2613

{5, 0, -2.82842712474619, 0, 0, 12.1243556529821, 0, 0, 0, 0}};

2614

2615

// Compute reference derivatives

2616

// Declare pointer to array of derivatives on FIAT element

2617

double *derivatives = new double [num_derivatives];

2618

2619

// Declare coefficients

2620

double coeff0_0 = 0;

2621

double coeff0_1 = 0;

2622

double coeff0_2 = 0;

2623

double coeff0_3 = 0;

2624

double coeff0_4 = 0;

2625

double coeff0_5 = 0;

2626

double coeff0_6 = 0;

2627

double coeff0_7 = 0;

2628

double coeff0_8 = 0;

2629

double coeff0_9 = 0;

2630

2631

// Declare new coefficients

2632

double new_coeff0_0 = 0;

2633

double new_coeff0_1 = 0;

2634

double new_coeff0_2 = 0;

2635

double new_coeff0_3 = 0;

2636

double new_coeff0_4 = 0;

2637

double new_coeff0_5 = 0;

2638

double new_coeff0_6 = 0;

2639

double new_coeff0_7 = 0;

2640

double new_coeff0_8 = 0;

2641

double new_coeff0_9 = 0;

2642

2643

// Loop possible derivatives

2644

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

2645

{

2646

// Get values from coefficients array

2647

new_coeff0_0 = coefficients0[dof][0];

2648

new_coeff0_1 = coefficients0[dof][1];

2649

new_coeff0_2 = coefficients0[dof][2];

2650

new_coeff0_3 = coefficients0[dof][3];

2651

new_coeff0_4 = coefficients0[dof][4];

2652

new_coeff0_5 = coefficients0[dof][5];

2653

new_coeff0_6 = coefficients0[dof][6];

2654

new_coeff0_7 = coefficients0[dof][7];

2655

new_coeff0_8 = coefficients0[dof][8];

2656

new_coeff0_9 = coefficients0[dof][9];

2657

2658

// Loop derivative order

2659

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

2660

{

2661

// Update old coefficients

2662

coeff0_0 = new_coeff0_0;

2663

coeff0_1 = new_coeff0_1;

2664

coeff0_2 = new_coeff0_2;

2665

coeff0_3 = new_coeff0_3;

2666

coeff0_4 = new_coeff0_4;

2667

coeff0_5 = new_coeff0_5;

2668

coeff0_6 = new_coeff0_6;

2669

coeff0_7 = new_coeff0_7;

2670

coeff0_8 = new_coeff0_8;

2671

coeff0_9 = new_coeff0_9;

2672

2673

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

2674

{

2675

2676

2677

2678

2679

2680

2681

2682

2683

2684

2685

}

2686

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

2687

{

2688

2689

2690

2691

2692

2693

2694

2695

2696

2697

2698

}

2699

2700

}

2701

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

2702

2703

}

2704

2705

// Transform derivatives back to physical element

2706

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

2707

{

2708

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

2709

{

2710

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

2711

}

2712

}

2713

// Delete pointer to array of derivatives on FIAT element

2714

delete [] derivatives;

2715

2716

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

2717

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

2718

{

2719

delete [] combinations[row];

2720

delete [] transform[row];

2721

}

2722

2723

delete [] combinations;

2724

delete [] transform;

2725

}

2726

2727

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

2728

virtual void evaluate_basis_derivatives_all(unsigned int n,

2729

double* values,

2730

const double* coordinates,

2731

const ufc::cell& c) const

2732

{

2733

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

2734

}

2735

2736

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

2737

virtual double evaluate_dof(unsigned int i,

2738

const ufc::function& f,

2739

const ufc::cell& c) const

2740

{

2741

// The reference points, direction and weights:

2742

2743

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

2744

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

2745

2746

const double * const * x = c.coordinates;

2747

double result = 0.0;

2748

// Iterate over the points:

2749

// Evaluate basis functions for affine mapping

2750

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

2751

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

2752

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

2753

2754

// Compute affine mapping y = F(X)

2755

double y[2];

2756

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

2757

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

2758

2759

// Evaluate function at physical points

2760

double values[1];

2761

f.evaluate(values, y, c);

2762

2763

// Map function values using appropriate mapping

2764

// Affine map: Do nothing

2765

2766

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

2767

2768

// Take directional components

2769

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

2770

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

2771

// Multiply by weights

2772

result *= W[i][0];

2773

2774

return result;

2775

}

2776

2777

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

2778

virtual void evaluate_dofs(double* values,

2779

const ufc::function& f,

2780

const ufc::cell& c) const

2781

{

2782

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

2783

}

2784

2785

/// Interpolate vertex values from dof values

2786

virtual void interpolate_vertex_values(double* vertex_values,

2787

const double* dof_values,

2788

const ufc::cell& c) const

2789

{

2790

// Evaluate at vertices and use affine mapping

2791

vertex_values[0] = dof_values[0];

2792

vertex_values[1] = dof_values[1];

2793

vertex_values[2] = dof_values[2];

2794

}

2795

2796

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

2797

virtual unsigned int num_sub_elements() const

2798

{

2799

return 1;

2800

}

2801

2802

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

2803

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

2804

{

2805

return new UFC_PoissonP3LinearForm_finite_element_1();

2806

}

2807

2808

};

2809

2810

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

2811

/// degrees of freedom (dofs).

2812

2813

class UFC_PoissonP3LinearForm_dof_map_0: public ufc::dof_map

2814

{

2815

private:

2816

2817

unsigned int __global_dimension;

2818

2819

public:

2820

2821

/// Constructor

2822

UFC_PoissonP3LinearForm_dof_map_0() : ufc::dof_map()

2823

{

2824

__global_dimension = 0;

2825

}

2826

2827

/// Destructor

2828

virtual ~UFC_PoissonP3LinearForm_dof_map_0()

2829

{

2830

// Do nothing

2831

}

2832

2833

/// Return a string identifying the dof map

2834

virtual const char* signature() const

2835

{

2836

return "FFC dof map for FiniteElement('Lagrange', 'triangle', 3)";

2837

}

2838

2839

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

2840

virtual bool needs_mesh_entities(unsigned int d) const

2841

{

2842

switch ( d )

2843

{

2844

case 0:

2845

return true;

2846

break;

2847

case 1:

2848

return true;

2849

break;

2850

case 2:

2851

return true;

2852

break;

2853

}

2854

return false;

2855

}

2856

2857

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

2858

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

2859

{

2860

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

2861

return false;

2862

}

2863

2864

/// Initialize dof map for given cell

2865

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

2866

const ufc::cell& c)

2867

{

2868

// Do nothing

2869

}

2870

2871

/// Finish initialization of dof map for cells

2872

virtual void init_cell_finalize()

2873

{

2874

// Do nothing

2875

}

2876

2877

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

2878

virtual unsigned int global_dimension() const

2879

{

2880

return __global_dimension;

2881

}

2882

2883

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

2884

virtual unsigned int local_dimension() const

2885

{

2886

return 10;

2887

}

2888

2889

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

2890

virtual unsigned int geometric_dimension() const

2891

{

2892

return 2;

2893

}

2894

2895

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

2896

virtual unsigned int num_facet_dofs() const

2897

{

2898

return 4;

2899

}

2900

2901

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

2902

virtual unsigned int num_entity_dofs(unsigned int d) const

2903

{

2904

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

2905

}

2906

2907

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

2908

virtual void tabulate_dofs(unsigned int* dofs,

2909

const ufc::mesh& m,

2910

const ufc::cell& c) const

2911

{

2912

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

2913

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

2914

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

2915

unsigned int offset = m.num_entities[0];

2916

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

2917

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

2918

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

2919

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

2920

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

2921

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

2922

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

2923

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

2924

}

2925

2926

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

2927

virtual void tabulate_facet_dofs(unsigned int* dofs,

2928

unsigned int facet) const

2929

{

2930

switch ( facet )

2931

{

2932

case 0:

2933

dofs[0] = 1;

2934

dofs[1] = 2;

2935

dofs[2] = 3;

2936

dofs[3] = 4;

2937

break;

2938

case 1:

2939

dofs[0] = 0;

2940

dofs[1] = 2;

2941

dofs[2] = 5;

2942

dofs[3] = 6;

2943

break;

2944

case 2:

2945

dofs[0] = 0;

2946

dofs[1] = 1;

2947

dofs[2] = 7;

2948

dofs[3] = 8;

2949

break;

2950

}

2951

}

2952

2953

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

2954

virtual void tabulate_entity_dofs(unsigned int* dofs,

2955

unsigned int d, unsigned int i) const

2956

{

2957

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

2958

}

2959

2960

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

2961

virtual void tabulate_coordinates(double** coordinates,

2962

const ufc::cell& c) const

2963

{

2964

const double * const * x = c.coordinates;

2965

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

2966

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

2967

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

2968

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

2969

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

2970

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

2971

coordinates[3][0] = 0.666666666666667*x[1][0] + 0.333333333333333*x[2][0];

2972

coordinates[3][1] = 0.666666666666667*x[1][1] + 0.333333333333333*x[2][1];

2973

coordinates[4][0] = 0.333333333333333*x[1][0] + 0.666666666666667*x[2][0];

2974

coordinates[4][1] = 0.333333333333333*x[1][1] + 0.666666666666667*x[2][1];

2975

coordinates[5][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[2][0];

2976

coordinates[5][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[2][1];

2977

coordinates[6][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[2][0];

2978

coordinates[6][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[2][1];

2979

coordinates[7][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[1][0];

2980

coordinates[7][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[1][1];

2981

coordinates[8][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[1][0];

2982

coordinates[8][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[1][1];

2983

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

2984

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

2985

}

2986

2987

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

2988

virtual unsigned int num_sub_dof_maps() const

2989

{

2990

return 1;

2991

}

2992

2993

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

2994

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

2995

{

2996

return new UFC_PoissonP3LinearForm_dof_map_0();

2997

}

2998

2999

};

3000

3001

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

3002

/// degrees of freedom (dofs).

3003

3004

class UFC_PoissonP3LinearForm_dof_map_1: public ufc::dof_map

3005

{

3006

private:

3007

3008

unsigned int __global_dimension;

3009

3010

public:

3011

3012

/// Constructor

3013

UFC_PoissonP3LinearForm_dof_map_1() : ufc::dof_map()

3014

{

3015

__global_dimension = 0;

3016

}

3017

3018

/// Destructor

3019

virtual ~UFC_PoissonP3LinearForm_dof_map_1()

3020

{

3021

// Do nothing

3022

}

3023

3024

/// Return a string identifying the dof map

3025

virtual const char* signature() const

3026

{

3027

return "FFC dof map for FiniteElement('Lagrange', 'triangle', 3)";

3028

}

3029

3030

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

3031

virtual bool needs_mesh_entities(unsigned int d) const

3032

{

3033

switch ( d )

3034

{

3035

case 0:

3036

return true;

3037

break;

3038

case 1:

3039

return true;

3040

break;

3041

case 2:

3042

return true;

3043

break;

3044

}

3045

return false;

3046

}

3047

3048

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

3049

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

3050

{

3051

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

3052

return false;

3053

}

3054

3055

/// Initialize dof map for given cell

3056

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

3057

const ufc::cell& c)

3058

{

3059

// Do nothing

3060

}

3061

3062

/// Finish initialization of dof map for cells

3063

virtual void init_cell_finalize()

3064

{

3065

// Do nothing

3066

}

3067

3068

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

3069

virtual unsigned int global_dimension() const

3070

{

3071

return __global_dimension;

3072

}

3073

3074

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

3075

virtual unsigned int local_dimension() const

3076

{

3077

return 10;

3078

}

3079

3080

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

3081

virtual unsigned int geometric_dimension() const

3082

{

3083

return 2;

3084

}

3085

3086

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

3087

virtual unsigned int num_facet_dofs() const

3088

{

3089

return 4;

3090

}

3091

3092

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

3093

virtual unsigned int num_entity_dofs(unsigned int d) const

3094

{

3095

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

3096

}

3097

3098

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

3099

virtual void tabulate_dofs(unsigned int* dofs,

3100

const ufc::mesh& m,

3101

const ufc::cell& c) const

3102

{

3103

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

3104

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

3105

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

3106

unsigned int offset = m.num_entities[0];

3107

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

3108

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

3109

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

3110

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

3111

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

3112

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

3113

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

3114

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

3115

}

3116

3117

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

3118

virtual void tabulate_facet_dofs(unsigned int* dofs,

3119

unsigned int facet) const

3120

{

3121

switch ( facet )

3122

{

3123

case 0:

3124

dofs[0] = 1;

3125

dofs[1] = 2;

3126

dofs[2] = 3;

3127

dofs[3] = 4;

3128

break;

3129

case 1:

3130

dofs[0] = 0;

3131

dofs[1] = 2;

3132

dofs[2] = 5;

3133

dofs[3] = 6;

3134

break;

3135

case 2:

3136

dofs[0] = 0;

3137

dofs[1] = 1;

3138

dofs[2] = 7;

3139

dofs[3] = 8;

3140

break;

3141

}

3142

}

3143

3144

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

3145

virtual void tabulate_entity_dofs(unsigned int* dofs,

3146

unsigned int d, unsigned int i) const

3147

{

3148

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

3149

}

3150

3151

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

3152

virtual void tabulate_coordinates(double** coordinates,

3153

const ufc::cell& c) const

3154

{

3155

const double * const * x = c.coordinates;

3156

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

3157

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

3158

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

3159

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

3160

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

3161

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

3162

coordinates[3][0] = 0.666666666666667*x[1][0] + 0.333333333333333*x[2][0];

3163

coordinates[3][1] = 0.666666666666667*x[1][1] + 0.333333333333333*x[2][1];

3164

coordinates[4][0] = 0.333333333333333*x[1][0] + 0.666666666666667*x[2][0];

3165

coordinates[4][1] = 0.333333333333333*x[1][1] + 0.666666666666667*x[2][1];

3166

coordinates[5][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[2][0];

3167

coordinates[5][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[2][1];

3168

coordinates[6][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[2][0];

3169

coordinates[6][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[2][1];

3170

coordinates[7][0] = 0.666666666666667*x[0][0] + 0.333333333333333*x[1][0];

3171

coordinates[7][1] = 0.666666666666667*x[0][1] + 0.333333333333333*x[1][1];

3172

coordinates[8][0] = 0.333333333333333*x[0][0] + 0.666666666666667*x[1][0];

3173

coordinates[8][1] = 0.333333333333333*x[0][1] + 0.666666666666667*x[1][1];

3174

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

3175

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

3176

}

3177

3178

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

3179

virtual unsigned int num_sub_dof_maps() const

3180

{

3181

return 1;

3182

}

3183

3184

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

3185

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

3186

{

3187

return new UFC_PoissonP3LinearForm_dof_map_1();

3188

}

3189

3190

};

3191

3192

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

3193

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

3194

/// the integral over a cell.

3195

3196

class UFC_PoissonP3LinearForm_cell_integral_0: public ufc::cell_integral

3197

{

3198

public:

3199

3200

/// Constructor

3201

UFC_PoissonP3LinearForm_cell_integral_0() : ufc::cell_integral()

3202

{

3203

// Do nothing

3204

}

3205

3206

/// Destructor

3207

virtual ~UFC_PoissonP3LinearForm_cell_integral_0()

3208

{

3209

// Do nothing

3210

}

3211

3212

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

3213

virtual void tabulate_tensor(double* A,

3214

const double * const * w,

3215

const ufc::cell& c) const

3216

{

3217

// Extract vertex coordinates

3218

const double * const * x = c.coordinates;

3219

3220

// Compute Jacobian of affine map from reference cell

3221

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

3222

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

3223

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

3224

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

3225

3226

// Compute determinant of Jacobian

3227

double detJ = J_00*J_11 - J_01*J_10;

3228

3229

// Compute inverse of Jacobian

3230

3231

// Set scale factor

3232

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

3233

3234

// Number of operations to compute element tensor = 176

3235

// Compute coefficients

3236

const double c0_0_0_0 = w[0][0];

3237

const double c0_0_0_1 = w[0][1];

3238

const double c0_0_0_2 = w[0][2];

3239

const double c0_0_0_3 = w[0][3];

3240

const double c0_0_0_4 = w[0][4];

3241

const double c0_0_0_5 = w[0][5];

3242

const double c0_0_0_6 = w[0][6];

3243

const double c0_0_0_7 = w[0][7];

3244

const double c0_0_0_8 = w[0][8];

3245

const double c0_0_0_9 = w[0][9];

3246

3247

// Compute geometry tensors

3248

// Number of operations to compute decalrations = 10

3249

const double G0_0 = det*c0_0_0_0;

3250

const double G0_1 = det*c0_0_0_1;

3251

const double G0_2 = det*c0_0_0_2;

3252

const double G0_3 = det*c0_0_0_3;

3253

const double G0_4 = det*c0_0_0_4;

3254

const double G0_5 = det*c0_0_0_5;

3255

const double G0_6 = det*c0_0_0_6;

3256

const double G0_7 = det*c0_0_0_7;

3257

const double G0_8 = det*c0_0_0_8;

3258

const double G0_9 = det*c0_0_0_9;

3259

3260

// Compute element tensor

3261

// Number of operations to compute tensor = 166

3262

A[0] = 0.0056547619047619*G0_0 + 0.000818452380952379*G0_1 + 0.000818452380952379*G0_2 + 0.00200892857142857*G0_3 + 0.00200892857142857*G0_4 + 0.0013392857142857*G0_5 + 0.0013392857142857*G0_7 + 0.00267857142857142*G0_9;

3263

A[1] = 0.000818452380952379*G0_0 + 0.0056547619047619*G0_1 + 0.00081845238095238*G0_2 + 0.0013392857142857*G0_3 + 0.00200892857142857*G0_5 + 0.00200892857142858*G0_6 + 0.0013392857142857*G0_8 + 0.00267857142857143*G0_9;

3264

A[2] = 0.000818452380952379*G0_0 + 0.00081845238095238*G0_1 + 0.0056547619047619*G0_2 + 0.00133928571428572*G0_4 + 0.00133928571428571*G0_6 + 0.00200892857142857*G0_7 + 0.00200892857142857*G0_8 + 0.00267857142857142*G0_9;

3265

A[3] = 0.00200892857142857*G0_0 + 0.0013392857142857*G0_1 + 0.0401785714285714*G0_3 - 0.0140625*G0_4 - 0.00401785714285713*G0_5 - 0.0100446428571429*G0_6 - 0.0100446428571428*G0_7 + 0.0200892857142857*G0_8 + 0.0120535714285714*G0_9;

3266

A[4] = 0.00200892857142857*G0_0 + 0.00133928571428572*G0_2 - 0.0140625*G0_3 + 0.0401785714285714*G0_4 - 0.0100446428571429*G0_5 + 0.0200892857142857*G0_6 - 0.00401785714285714*G0_7 - 0.0100446428571428*G0_8 + 0.0120535714285714*G0_9;

3267

A[5] = 0.0013392857142857*G0_0 + 0.00200892857142857*G0_1 - 0.00401785714285713*G0_3 - 0.0100446428571429*G0_4 + 0.0401785714285714*G0_5 - 0.0140625*G0_6 + 0.0200892857142857*G0_7 - 0.0100446428571428*G0_8 + 0.0120535714285714*G0_9;

3268

A[6] = 0.00200892857142858*G0_1 + 0.00133928571428571*G0_2 - 0.0100446428571429*G0_3 + 0.0200892857142857*G0_4 - 0.0140625*G0_5 + 0.0401785714285714*G0_6 - 0.0100446428571429*G0_7 - 0.00401785714285714*G0_8 + 0.0120535714285714*G0_9;

3269

A[7] = 0.0013392857142857*G0_0 + 0.00200892857142857*G0_2 - 0.0100446428571428*G0_3 - 0.00401785714285714*G0_4 + 0.0200892857142857*G0_5 - 0.0100446428571429*G0_6 + 0.0401785714285714*G0_7 - 0.0140625*G0_8 + 0.0120535714285714*G0_9;

3270

A[8] = 0.0013392857142857*G0_1 + 0.00200892857142857*G0_2 + 0.0200892857142857*G0_3 - 0.0100446428571428*G0_4 - 0.0100446428571428*G0_5 - 0.00401785714285714*G0_6 - 0.0140625*G0_7 + 0.0401785714285714*G0_8 + 0.0120535714285714*G0_9;

3271

A[9] = 0.00267857142857142*G0_0 + 0.00267857142857143*G0_1 + 0.00267857142857142*G0_2 + 0.0120535714285714*G0_3 + 0.0120535714285714*G0_4 + 0.0120535714285714*G0_5 + 0.0120535714285714*G0_6 + 0.0120535714285714*G0_7 + 0.0120535714285714*G0_8 + 0.144642857142857*G0_9;

3272

}

3273

3274

};

3275

3276

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

3277

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

3278

/// mapping

3279

///

3280

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

3281

///

3282

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

3283

/// global tensor A is defined by

3284

///

3285

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

3286

///

3287

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

3288

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

3289

/// fixed functions (coefficients).

3290

3291

class UFC_PoissonP3LinearForm: public ufc::form

3292

{

3293

public:

3294

3295

/// Constructor

3296

UFC_PoissonP3LinearForm() : ufc::form()

3297

{

3298

// Do nothing

3299

}

3300

3301

/// Destructor

3302

virtual ~UFC_PoissonP3LinearForm()

3303

{

3304

// Do nothing

3305

}

3306

3307

/// Return a string identifying the form

3308

virtual const char* signature() const

3309

{

3310

return "w0_a0[0, 1, 2, 3, 4, 5, 6, 7, 8, 9] | vi0[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]*va0[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]*dX(0)";

3311

}

3312

3313

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

3314

virtual unsigned int rank() const

3315

{

3316

return 1;

3317

}

3318

3319

/// Return the number of coefficients (n)

3320

virtual unsigned int num_coefficients() const

3321

{

3322

return 1;

3323

}

3324

3325

/// Return the number of cell integrals

3326

virtual unsigned int num_cell_integrals() const

3327

{

3328

return 1;

3329

}

3330

3331

/// Return the number of exterior facet integrals

3332

virtual unsigned int num_exterior_facet_integrals() const

3333

{

3334

return 0;

3335

}

3336

3337

/// Return the number of interior facet integrals

3338

virtual unsigned int num_interior_facet_integrals() const

3339

{

3340

return 0;

3341

}

3342

3343

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

3344

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

3345

{

3346

switch ( i )

3347

{

3348

case 0:

3349

return new UFC_PoissonP3LinearForm_finite_element_0();

3350

break;

3351

case 1:

3352

return new UFC_PoissonP3LinearForm_finite_element_1();

3353

break;

3354

}

3355

return 0;

3356

}

3357

3358

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

3359

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

3360

{

3361

switch ( i )

3362

{

3363

case 0:

3364

return new UFC_PoissonP3LinearForm_dof_map_0();

3365

break;

3366

case 1:

3367

return new UFC_PoissonP3LinearForm_dof_map_1();

3368

break;

3369

}

3370

return 0;

3371

}

3372

3373

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

3374

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

3375

{

3376

return new UFC_PoissonP3LinearForm_cell_integral_0();

3377

}

3378

3379

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

3380

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

3381

{

3382

return 0;

3383

}

3384

3385

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

3386

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

3387

{

3388

return 0;

3389

}

3390

3391

};

3392

3393

// DOLFIN wrappers

3394

3395

#include <dolfin/fem/Form.h>

3396

#include <dolfin/fem/FiniteElement.h>

3397

#include <dolfin/fem/DofMap.h>

3398

#include <dolfin/function/Coefficient.h>

3399

#include <dolfin/function/Function.h>

3400

#include <dolfin/function/FunctionSpace.h>

3401

3402

class PoissonP3BilinearFormFunctionSpace0 : public dolfin::FunctionSpace

3403

{

3404

public:

3405

3406

PoissonP3BilinearFormFunctionSpace0(const dolfin::Mesh& mesh)

3407

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3408

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_PoissonP3LinearForm_finite_element_1()))),

3409

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_PoissonP3LinearForm_dof_map_1()), mesh)))

3410

{

3411

// Do nothing

3412

}

3413

3414

};

3415

3416

class PoissonP3BilinearFormFunctionSpace1 : public dolfin::FunctionSpace

3417

{

3418

public:

3419

3420

PoissonP3BilinearFormFunctionSpace1(const dolfin::Mesh& mesh)

3421

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3422

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_PoissonP3LinearForm_finite_element_1()))),

3423

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_PoissonP3LinearForm_dof_map_1()), mesh)))

3424

{

3425

// Do nothing

3426

}

3427

3428

};

3429

3430

class PoissonP3LinearFormFunctionSpace0 : public dolfin::FunctionSpace

3431

{

3432

public:

3433

3434

PoissonP3LinearFormFunctionSpace0(const dolfin::Mesh& mesh)

3435

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3436

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_PoissonP3LinearForm_finite_element_1()))),

3437

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_PoissonP3LinearForm_dof_map_1()), mesh)))

3438

{

3439

// Do nothing

3440

}

3441

3442

};

3443

3444

class PoissonP3LinearFormCoefficientSpace0 : public dolfin::FunctionSpace

3445

{

3446

public:

3447

3448

PoissonP3LinearFormCoefficientSpace0(const dolfin::Mesh& mesh)

3449

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3450

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_PoissonP3LinearForm_finite_element_1()))),

3451

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_PoissonP3LinearForm_dof_map_1()), mesh)))

3452

{

3453

// Do nothing

3454

}

3455

3456

};

3457

3458

class PoissonP3TestSpace : public dolfin::FunctionSpace

3459

{

3460

public:

3461

3462

PoissonP3TestSpace(const dolfin::Mesh& mesh)

3463

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3464

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_PoissonP3LinearForm_finite_element_1()))),

3465

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_PoissonP3LinearForm_dof_map_1()), mesh)))

3466

{

3467

// Do nothing

3468

}

3469

3470

};

3471

3472

class PoissonP3TrialSpace : public dolfin::FunctionSpace

3473

{

3474

public:

3475

3476

PoissonP3TrialSpace(const dolfin::Mesh& mesh)

3477

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3478

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_PoissonP3LinearForm_finite_element_1()))),

3479

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_PoissonP3LinearForm_dof_map_1()), mesh)))

3480

{

3481

// Do nothing

3482

}

3483

3484

};

3485

3486

class PoissonP3CoefficientSpace : public dolfin::FunctionSpace

3487

{

3488

public:

3489

3490

PoissonP3CoefficientSpace(const dolfin::Mesh& mesh)

3491

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3492

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_PoissonP3LinearForm_finite_element_1()))),

3493

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_PoissonP3LinearForm_dof_map_1()), mesh)))

3494

{

3495

// Do nothing

3496

}

3497

3498

};

3499

3500

class PoissonP3FunctionSpace : public dolfin::FunctionSpace

3501

{

3502

public:

3503

3504

PoissonP3FunctionSpace(const dolfin::Mesh& mesh)

3505

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3506

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_PoissonP3LinearForm_finite_element_1()))),

3507

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_PoissonP3LinearForm_dof_map_1()), mesh)))

3508

{

3509

// Do nothing

3510

}

3511

3512

};

3513

3514

class PoissonP3BilinearForm : public dolfin::Form

3515

{

3516

public:

3517

3518

// Create form on given function space(s)

3519

PoissonP3BilinearForm(const dolfin::FunctionSpace& V0, const dolfin::FunctionSpace& V1) : dolfin::Form()

3520

{

3521

boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());

3522

_function_spaces.push_back(_V0);

3523

boost::shared_ptr<const dolfin::FunctionSpace> _V1(&V1, dolfin::NoDeleter<const dolfin::FunctionSpace>());

3524

_function_spaces.push_back(_V1);

3525

3526

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_PoissonP3BilinearForm());

3527

}

3528

3529

// Create form on given function space(s) (shared data)

3530

PoissonP3BilinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0, boost::shared_ptr<const dolfin::FunctionSpace> V1) : dolfin::Form()

3531

{

3532

_function_spaces.push_back(V0);

3533

_function_spaces.push_back(V1);

3534

3535

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_PoissonP3BilinearForm());

3536

}

3537

3538

// Destructor

3539

~PoissonP3BilinearForm() {}

3540

3541

};

3542

3543

class PoissonP3LinearFormCoefficient0 : public dolfin::Coefficient

3544

{

3545

public:

3546

3547

// Constructor

3548

PoissonP3LinearFormCoefficient0(dolfin::Form& form) : dolfin::Coefficient(form) {}

3549

3550

// Destructor

3551

~PoissonP3LinearFormCoefficient0() {}

3552

3553

// Attach function to coefficient

3554

const PoissonP3LinearFormCoefficient0& operator= (dolfin::Function& v)

3555

{

3556

attach(v);

3557

return *this;

3558

}

3559

3560

/// Create function space for coefficient

3561

const dolfin::FunctionSpace* create_function_space() const

3562

{

3563

return new PoissonP3LinearFormCoefficientSpace0(form.mesh());

3564

}

3565

3566

/// Return coefficient number

3567

dolfin::uint number() const

3568

{

3569

return 0;

3570

}

3571

3572

/// Return coefficient name

3573

virtual std::string name() const

3574

{

3575

return "f";

3576

}

3577

3578

};

3579

class PoissonP3LinearForm : public dolfin::Form

3580

{

3581

public:

3582

3583

// Create form on given function space(s)

3584

PoissonP3LinearForm(const dolfin::FunctionSpace& V0) : dolfin::Form(), f(*this)

3585

{

3586

boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());

3587

_function_spaces.push_back(_V0);

3588

3589

_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));

3590

3591

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_PoissonP3LinearForm());

3592

}

3593

3594

// Create form on given function space(s) (shared data)

3595

PoissonP3LinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0) : dolfin::Form(), f(*this)

3596

{

3597

_function_spaces.push_back(V0);

3598

3599

_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));

3600

3601

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_PoissonP3LinearForm());

3602

}

3603

3604

// Create form on given function space(s) with given coefficient(s)

3605

PoissonP3LinearForm(const dolfin::FunctionSpace& V0, dolfin::Function& w0) : dolfin::Form(), f(*this)

3606

{

3607

boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());

3608

_function_spaces.push_back(_V0);

3609

3610

_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));

3611

3612

this->f = w0;

3613

3614

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_PoissonP3LinearForm());

3615

}

3616

3617

// Create form on given function space(s) with given coefficient(s) (shared data)

3618

PoissonP3LinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0, dolfin::Function& w0) : dolfin::Form(), f(*this)

3619

{

3620

_function_spaces.push_back(V0);

3621

3622

_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));

3623

3624

this->f = w0;

3625

3626

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_PoissonP3LinearForm());

3627

}

3628

3629

// Destructor

3630

~PoissonP3LinearForm() {}

3631

3632

// Coefficients

3633

PoissonP3LinearFormCoefficient0 f;

3634

3635

};

3636

3637

#endif

Older »