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

1

// This code conforms with the UFC specification version 1.0

2

// and was automatically generated by FFC version 0.7.0.

3

4

#ifndef __POISSONDG_H

5

#define __POISSONDG_H

6

7

#include <cmath>

8

#include <stdexcept>

9

#include <ufc.h>

10

11

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

12

13

class poissondg_0_finite_element_0: public ufc::finite_element

14

{

15

public:

16

17

/// Constructor

18

poissondg_0_finite_element_0() : ufc::finite_element()

19

{

20

// Do nothing

21

}

22

23

/// Destructor

24

virtual ~poissondg_0_finite_element_0()

25

{

26

// Do nothing

27

}

28

29

/// Return a string identifying the finite element

30

virtual const char* signature() const

31

{

32

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

33

}

34

35

/// Return the cell shape

36

virtual ufc::shape cell_shape() const

37

{

38

return ufc::triangle;

39

}

40

41

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

42

virtual unsigned int space_dimension() const

43

{

44

return 3;

45

}

46

47

/// Return the rank of the value space

48

virtual unsigned int value_rank() const

49

{

50

return 0;

51

}

52

53

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

54

virtual unsigned int value_dimension(unsigned int i) const

55

{

56

return 1;

57

}

58

59

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

60

virtual void evaluate_basis(unsigned int i,

61

double* values,

62

const double* coordinates,

63

const ufc::cell& c) const

64

{

65

// Extract vertex coordinates

66

const double * const * element_coordinates = c.coordinates;

67

68

// Compute Jacobian of affine map from reference cell

69

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

70

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

71

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

72

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

73

74

// Compute determinant of Jacobian

75

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

76

77

// Compute inverse of Jacobian

78

79

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

80

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

81

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

82

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

83

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

84

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

85

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

86

87

// Map coordinates to the reference square

88

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

89

x = -1.0;

90

else

91

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

92

y = 2.0*y - 1.0;

93

94

// Reset values

95

*values = 0;

96

97

// Map degree of freedom to element degree of freedom

98

const unsigned int dof = i;

99

100

// Generate scalings

101

const double scalings_y_0 = 1;

102

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

103

104

// Compute psitilde_a

105

const double psitilde_a_0 = 1;

106

const double psitilde_a_1 = x;

107

108

// Compute psitilde_bs

109

const double psitilde_bs_0_0 = 1;

110

const double psitilde_bs_0_1 = 1.5*y + 0.5;

111

const double psitilde_bs_1_0 = 1;

112

113

// Compute basisvalues

114

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

115

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

116

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

117

118

// Table(s) of coefficients

119

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

120

{{0.471404521, -0.288675135, -0.166666667},

121

{0.471404521, 0.288675135, -0.166666667},

122

{0.471404521, 0, 0.333333333}};

123

124

// Extract relevant coefficients

125

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

126

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

127

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

128

129

// Compute value(s)

130

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;

131

}

132

133

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

134

virtual void evaluate_basis_all(double* values,

135

const double* coordinates,

136

const ufc::cell& c) const

137

{

138

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

139

}

140

141

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

142

virtual void evaluate_basis_derivatives(unsigned int i,

143

unsigned int n,

144

double* values,

145

const double* coordinates,

146

const ufc::cell& c) const

147

{

148

// Extract vertex coordinates

149

const double * const * element_coordinates = c.coordinates;

150

151

// Compute Jacobian of affine map from reference cell

152

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

153

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

154

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

155

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

156

157

// Compute determinant of Jacobian

158

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

159

160

// Compute inverse of Jacobian

161

162

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

163

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

164

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

165

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

166

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

167

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

168

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

169

170

// Map coordinates to the reference square

171

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

172

x = -1.0;

173

else

174

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

175

y = 2.0*y - 1.0;

176

177

// Compute number of derivatives

178

unsigned int num_derivatives = 1;

179

180

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

181

num_derivatives *= 2;

182

183

184

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

185

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

186

187

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

188

{

189

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

190

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

191

combinations[j][k] = 0;

192

}

193

194

// Generate combinations of derivatives

195

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

196

{

197

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

198

{

199

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

200

{

201

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

202

combinations[row][col] = 0;

203

else

204

{

205

combinations[row][col] += 1;

206

break;

207

}

208

}

209

}

210

}

211

212

// Compute inverse of Jacobian

213

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

214

215

// Declare transformation matrix

216

// Declare pointer to two dimensional array and initialise

217

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

218

219

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

220

{

221

transform[j] = new double [num_derivatives];

222

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

223

transform[j][k] = 1;

224

}

225

226

// Construct transformation matrix

227

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

228

{

229

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

230

{

231

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

232

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

233

}

234

}

235

236

// Reset values

237

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

238

values[j] = 0;

239

240

// Map degree of freedom to element degree of freedom

241

const unsigned int dof = i;

242

243

// Generate scalings

244

const double scalings_y_0 = 1;

245

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

246

247

// Compute psitilde_a

248

const double psitilde_a_0 = 1;

249

const double psitilde_a_1 = x;

250

251

// Compute psitilde_bs

252

const double psitilde_bs_0_0 = 1;

253

const double psitilde_bs_0_1 = 1.5*y + 0.5;

254

const double psitilde_bs_1_0 = 1;

255

256

// Compute basisvalues

257

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

258

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

259

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

260

261

// Table(s) of coefficients

262

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

263

{{0.471404521, -0.288675135, -0.166666667},

264

{0.471404521, 0.288675135, -0.166666667},

265

{0.471404521, 0, 0.333333333}};

266

267

// Interesting (new) part

268

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

269

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

270

{{0, 0, 0},

271

{4.89897949, 0, 0},

272

{0, 0, 0}};

273

274

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

275

{{0, 0, 0},

276

{2.44948974, 0, 0},

277

{4.24264069, 0, 0}};

278

279

// Compute reference derivatives

280

// Declare pointer to array of derivatives on FIAT element

281

double *derivatives = new double [num_derivatives];

282

283

// Declare coefficients

284

double coeff0_0 = 0;

285

double coeff0_1 = 0;

286

double coeff0_2 = 0;

287

288

// Declare new coefficients

289

double new_coeff0_0 = 0;

290

double new_coeff0_1 = 0;

291

double new_coeff0_2 = 0;

292

293

// Loop possible derivatives

294

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

295

{

296

// Get values from coefficients array

297

new_coeff0_0 = coefficients0[dof][0];

298

new_coeff0_1 = coefficients0[dof][1];

299

new_coeff0_2 = coefficients0[dof][2];

300

301

// Loop derivative order

302

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

303

{

304

// Update old coefficients

305

coeff0_0 = new_coeff0_0;

306

coeff0_1 = new_coeff0_1;

307

coeff0_2 = new_coeff0_2;

308

309

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

310

{

311

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];

312

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];

313

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];

314

}

315

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

316

{

317

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];

318

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];

319

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];

320

}

321

322

}

323

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

324

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;

325

}

326

327

// Transform derivatives back to physical element

328

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

329

{

330

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

331

{

332

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

333

}

334

}

335

// Delete pointer to array of derivatives on FIAT element

336

delete [] derivatives;

337

338

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

339

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

340

{

341

delete [] combinations[row];

342

delete [] transform[row];

343

}

344

345

delete [] combinations;

346

delete [] transform;

347

}

348

349

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

350

virtual void evaluate_basis_derivatives_all(unsigned int n,

351

double* values,

352

const double* coordinates,

353

const ufc::cell& c) const

354

{

355

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

356

}

357

358

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

359

virtual double evaluate_dof(unsigned int i,

360

const ufc::function& f,

361

const ufc::cell& c) const

362

{

363

// The reference points, direction and weights:

364

static const double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};

365

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

366

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

367

368

const double * const * x = c.coordinates;

369

double result = 0.0;

370

// Iterate over the points:

371

// Evaluate basis functions for affine mapping

372

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

373

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

374

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

375

376

// Compute affine mapping y = F(X)

377

double y[2];

378

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

379

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

380

381

// Evaluate function at physical points

382

double values[1];

383

f.evaluate(values, y, c);

384

385

// Map function values using appropriate mapping

386

// Affine map: Do nothing

387

388

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

389

390

// Take directional components

391

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

392

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

393

// Multiply by weights

394

result *= W[i][0];

395

396

return result;

397

}

398

399

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

400

virtual void evaluate_dofs(double* values,

401

const ufc::function& f,

402

const ufc::cell& c) const

403

{

404

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

405

}

406

407

/// Interpolate vertex values from dof values

408

virtual void interpolate_vertex_values(double* vertex_values,

409

const double* dof_values,

410

const ufc::cell& c) const

411

{

412

// Evaluate at vertices and use affine mapping

413

vertex_values[0] = dof_values[0];

414

vertex_values[1] = dof_values[1];

415

vertex_values[2] = dof_values[2];

416

}

417

418

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

419

virtual unsigned int num_sub_elements() const

420

{

421

return 1;

422

}

423

424

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

425

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

426

{

427

return new poissondg_0_finite_element_0();

428

}

429

430

};

431

432

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

433

434

class poissondg_0_finite_element_1: public ufc::finite_element

435

{

436

public:

437

438

/// Constructor

439

poissondg_0_finite_element_1() : ufc::finite_element()

440

{

441

// Do nothing

442

}

443

444

/// Destructor

445

virtual ~poissondg_0_finite_element_1()

446

{

447

// Do nothing

448

}

449

450

/// Return a string identifying the finite element

451

virtual const char* signature() const

452

{

453

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

454

}

455

456

/// Return the cell shape

457

virtual ufc::shape cell_shape() const

458

{

459

return ufc::triangle;

460

}

461

462

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

463

virtual unsigned int space_dimension() const

464

{

465

return 3;

466

}

467

468

/// Return the rank of the value space

469

virtual unsigned int value_rank() const

470

{

471

return 0;

472

}

473

474

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

475

virtual unsigned int value_dimension(unsigned int i) const

476

{

477

return 1;

478

}

479

480

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

481

virtual void evaluate_basis(unsigned int i,

482

double* values,

483

const double* coordinates,

484

const ufc::cell& c) const

485

{

486

// Extract vertex coordinates

487

const double * const * element_coordinates = c.coordinates;

488

489

// Compute Jacobian of affine map from reference cell

490

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

491

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

492

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

493

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

494

495

// Compute determinant of Jacobian

496

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

497

498

// Compute inverse of Jacobian

499

500

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

501

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

502

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

503

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

504

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

505

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

506

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

507

508

// Map coordinates to the reference square

509

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

510

x = -1.0;

511

else

512

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

513

y = 2.0*y - 1.0;

514

515

// Reset values

516

*values = 0;

517

518

// Map degree of freedom to element degree of freedom

519

const unsigned int dof = i;

520

521

// Generate scalings

522

const double scalings_y_0 = 1;

523

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

524

525

// Compute psitilde_a

526

const double psitilde_a_0 = 1;

527

const double psitilde_a_1 = x;

528

529

// Compute psitilde_bs

530

const double psitilde_bs_0_0 = 1;

531

const double psitilde_bs_0_1 = 1.5*y + 0.5;

532

const double psitilde_bs_1_0 = 1;

533

534

// Compute basisvalues

535

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

536

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

537

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

538

539

// Table(s) of coefficients

540

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

541

{{0.471404521, -0.288675135, -0.166666667},

542

{0.471404521, 0.288675135, -0.166666667},

543

{0.471404521, 0, 0.333333333}};

544

545

// Extract relevant coefficients

546

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

547

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

548

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

549

550

// Compute value(s)

551

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;

552

}

553

554

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

555

virtual void evaluate_basis_all(double* values,

556

const double* coordinates,

557

const ufc::cell& c) const

558

{

559

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

560

}

561

562

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

563

virtual void evaluate_basis_derivatives(unsigned int i,

564

unsigned int n,

565

double* values,

566

const double* coordinates,

567

const ufc::cell& c) const

568

{

569

// Extract vertex coordinates

570

const double * const * element_coordinates = c.coordinates;

571

572

// Compute Jacobian of affine map from reference cell

573

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

574

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

575

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

576

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

577

578

// Compute determinant of Jacobian

579

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

580

581

// Compute inverse of Jacobian

582

583

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

584

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

585

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

586

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

587

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

588

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

589

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

590

591

// Map coordinates to the reference square

592

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

593

x = -1.0;

594

else

595

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

596

y = 2.0*y - 1.0;

597

598

// Compute number of derivatives

599

unsigned int num_derivatives = 1;

600

601

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

602

num_derivatives *= 2;

603

604

605

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

606

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

607

608

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

609

{

610

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

611

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

612

combinations[j][k] = 0;

613

}

614

615

// Generate combinations of derivatives

616

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

617

{

618

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

619

{

620

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

621

{

622

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

623

combinations[row][col] = 0;

624

else

625

{

626

combinations[row][col] += 1;

627

break;

628

}

629

}

630

}

631

}

632

633

// Compute inverse of Jacobian

634

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

635

636

// Declare transformation matrix

637

// Declare pointer to two dimensional array and initialise

638

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

639

640

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

641

{

642

transform[j] = new double [num_derivatives];

643

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

644

transform[j][k] = 1;

645

}

646

647

// Construct transformation matrix

648

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

649

{

650

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

651

{

652

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

653

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

654

}

655

}

656

657

// Reset values

658

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

659

values[j] = 0;

660

661

// Map degree of freedom to element degree of freedom

662

const unsigned int dof = i;

663

664

// Generate scalings

665

const double scalings_y_0 = 1;

666

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

667

668

// Compute psitilde_a

669

const double psitilde_a_0 = 1;

670

const double psitilde_a_1 = x;

671

672

// Compute psitilde_bs

673

const double psitilde_bs_0_0 = 1;

674

const double psitilde_bs_0_1 = 1.5*y + 0.5;

675

const double psitilde_bs_1_0 = 1;

676

677

// Compute basisvalues

678

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

679

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

680

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

681

682

// Table(s) of coefficients

683

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

684

{{0.471404521, -0.288675135, -0.166666667},

685

{0.471404521, 0.288675135, -0.166666667},

686

{0.471404521, 0, 0.333333333}};

687

688

// Interesting (new) part

689

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

690

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

691

{{0, 0, 0},

692

{4.89897949, 0, 0},

693

{0, 0, 0}};

694

695

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

696

{{0, 0, 0},

697

{2.44948974, 0, 0},

698

{4.24264069, 0, 0}};

699

700

// Compute reference derivatives

701

// Declare pointer to array of derivatives on FIAT element

702

double *derivatives = new double [num_derivatives];

703

704

// Declare coefficients

705

double coeff0_0 = 0;

706

double coeff0_1 = 0;

707

double coeff0_2 = 0;

708

709

// Declare new coefficients

710

double new_coeff0_0 = 0;

711

double new_coeff0_1 = 0;

712

double new_coeff0_2 = 0;

713

714

// Loop possible derivatives

715

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

716

{

717

// Get values from coefficients array

718

new_coeff0_0 = coefficients0[dof][0];

719

new_coeff0_1 = coefficients0[dof][1];

720

new_coeff0_2 = coefficients0[dof][2];

721

722

// Loop derivative order

723

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

724

{

725

// Update old coefficients

726

coeff0_0 = new_coeff0_0;

727

coeff0_1 = new_coeff0_1;

728

coeff0_2 = new_coeff0_2;

729

730

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

731

{

732

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];

733

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];

734

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];

735

}

736

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

737

{

738

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];

739

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];

740

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];

741

}

742

743

}

744

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

745

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;

746

}

747

748

// Transform derivatives back to physical element

749

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

750

{

751

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

752

{

753

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

754

}

755

}

756

// Delete pointer to array of derivatives on FIAT element

757

delete [] derivatives;

758

759

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

760

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

761

{

762

delete [] combinations[row];

763

delete [] transform[row];

764

}

765

766

delete [] combinations;

767

delete [] transform;

768

}

769

770

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

771

virtual void evaluate_basis_derivatives_all(unsigned int n,

772

double* values,

773

const double* coordinates,

774

const ufc::cell& c) const

775

{

776

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

777

}

778

779

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

780

virtual double evaluate_dof(unsigned int i,

781

const ufc::function& f,

782

const ufc::cell& c) const

783

{

784

// The reference points, direction and weights:

785

static const double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};

786

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

787

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

788

789

const double * const * x = c.coordinates;

790

double result = 0.0;

791

// Iterate over the points:

792

// Evaluate basis functions for affine mapping

793

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

794

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

795

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

796

797

// Compute affine mapping y = F(X)

798

double y[2];

799

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

800

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

801

802

// Evaluate function at physical points

803

double values[1];

804

f.evaluate(values, y, c);

805

806

// Map function values using appropriate mapping

807

// Affine map: Do nothing

808

809

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

810

811

// Take directional components

812

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

813

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

814

// Multiply by weights

815

result *= W[i][0];

816

817

return result;

818

}

819

820

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

821

virtual void evaluate_dofs(double* values,

822

const ufc::function& f,

823

const ufc::cell& c) const

824

{

825

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

826

}

827

828

/// Interpolate vertex values from dof values

829

virtual void interpolate_vertex_values(double* vertex_values,

830

const double* dof_values,

831

const ufc::cell& c) const

832

{

833

// Evaluate at vertices and use affine mapping

834

vertex_values[0] = dof_values[0];

835

vertex_values[1] = dof_values[1];

836

vertex_values[2] = dof_values[2];

837

}

838

839

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

840

virtual unsigned int num_sub_elements() const

841

{

842

return 1;

843

}

844

845

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

846

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

847

{

848

return new poissondg_0_finite_element_1();

849

}

850

851

};

852

853

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

854

855

class poissondg_0_finite_element_2: public ufc::finite_element

856

{

857

public:

858

859

/// Constructor

860

poissondg_0_finite_element_2() : ufc::finite_element()

861

{

862

// Do nothing

863

}

864

865

/// Destructor

866

virtual ~poissondg_0_finite_element_2()

867

{

868

// Do nothing

869

}

870

871

/// Return a string identifying the finite element

872

virtual const char* signature() const

873

{

874

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

875

}

876

877

/// Return the cell shape

878

virtual ufc::shape cell_shape() const

879

{

880

return ufc::triangle;

881

}

882

883

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

884

virtual unsigned int space_dimension() const

885

{

886

return 1;

887

}

888

889

/// Return the rank of the value space

890

virtual unsigned int value_rank() const

891

{

892

return 0;

893

}

894

895

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

896

virtual unsigned int value_dimension(unsigned int i) const

897

{

898

return 1;

899

}

900

901

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

902

virtual void evaluate_basis(unsigned int i,

903

double* values,

904

const double* coordinates,

905

const ufc::cell& c) const

906

{

907

// Extract vertex coordinates

908

const double * const * element_coordinates = c.coordinates;

909

910

// Compute Jacobian of affine map from reference cell

911

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

912

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

913

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

914

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

915

916

// Compute determinant of Jacobian

917

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

918

919

// Compute inverse of Jacobian

920

921

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

922

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

923

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

924

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

925

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

926

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

927

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

928

929

// Map coordinates to the reference square

930

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

931

x = -1.0;

932

else

933

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

934

y = 2.0*y - 1.0;

935

936

// Reset values

937

*values = 0;

938

939

// Map degree of freedom to element degree of freedom

940

const unsigned int dof = i;

941

942

// Generate scalings

943

const double scalings_y_0 = 1;

944

945

// Compute psitilde_a

946

const double psitilde_a_0 = 1;

947

948

// Compute psitilde_bs

949

const double psitilde_bs_0_0 = 1;

950

951

// Compute basisvalues

952

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

953

954

// Table(s) of coefficients

955

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

956

{{1.41421356}};

957

958

// Extract relevant coefficients

959

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

960

961

// Compute value(s)

962

*values = coeff0_0*basisvalue0;

963

}

964

965

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

966

virtual void evaluate_basis_all(double* values,

967

const double* coordinates,

968

const ufc::cell& c) const

969

{

970

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

971

}

972

973

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

974

virtual void evaluate_basis_derivatives(unsigned int i,

975

unsigned int n,

976

double* values,

977

const double* coordinates,

978

const ufc::cell& c) const

979

{

980

// Extract vertex coordinates

981

const double * const * element_coordinates = c.coordinates;

982

983

// Compute Jacobian of affine map from reference cell

984

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

985

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

986

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

987

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

988

989

// Compute determinant of Jacobian

990

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

991

992

// Compute inverse of Jacobian

993

994

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

995

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

996

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

997

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

998

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

999

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

1000

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

1001

1002

// Map coordinates to the reference square

1003

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

1004

x = -1.0;

1005

else

1006

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

1007

y = 2.0*y - 1.0;

1008

1009

// Compute number of derivatives

1010

unsigned int num_derivatives = 1;

1011

1012

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

1013

num_derivatives *= 2;

1014

1015

1016

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

1017

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

1018

1019

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

1020

{

1021

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

1022

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

1023

combinations[j][k] = 0;

1024

}

1025

1026

// Generate combinations of derivatives

1027

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

1028

{

1029

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

1030

{

1031

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

1032

{

1033

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

1034

combinations[row][col] = 0;

1035

else

1036

{

1037

combinations[row][col] += 1;

1038

break;

1039

}

1040

}

1041

}

1042

}

1043

1044

// Compute inverse of Jacobian

1045

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

1046

1047

// Declare transformation matrix

1048

// Declare pointer to two dimensional array and initialise

1049

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

1050

1051

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

1052

{

1053

transform[j] = new double [num_derivatives];

1054

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

1055

transform[j][k] = 1;

1056

}

1057

1058

// Construct transformation matrix

1059

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

1060

{

1061

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

1062

{

1063

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

1064

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

1065

}

1066

}

1067

1068

// Reset values

1069

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

1070

values[j] = 0;

1071

1072

// Map degree of freedom to element degree of freedom

1073

const unsigned int dof = i;

1074

1075

// Generate scalings

1076

const double scalings_y_0 = 1;

1077

1078

// Compute psitilde_a

1079

const double psitilde_a_0 = 1;

1080

1081

// Compute psitilde_bs

1082

const double psitilde_bs_0_0 = 1;

1083

1084

// Compute basisvalues

1085

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

1086

1087

// Table(s) of coefficients

1088

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

1089

{{1.41421356}};

1090

1091

// Interesting (new) part

1092

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

1093

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

1094

{{0}};

1095

1096

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

1097

{{0}};

1098

1099

// Compute reference derivatives

1100

// Declare pointer to array of derivatives on FIAT element

1101

double *derivatives = new double [num_derivatives];

1102

1103

// Declare coefficients

1104

double coeff0_0 = 0;

1105

1106

// Declare new coefficients

1107

double new_coeff0_0 = 0;

1108

1109

// Loop possible derivatives

1110

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

1111

{

1112

// Get values from coefficients array

1113

new_coeff0_0 = coefficients0[dof][0];

1114

1115

// Loop derivative order

1116

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

1117

{

1118

// Update old coefficients

1119

coeff0_0 = new_coeff0_0;

1120

1121

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

1122

{

1123

new_coeff0_0 = coeff0_0*dmats0[0][0];

1124

}

1125

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

1126

{

1127

new_coeff0_0 = coeff0_0*dmats1[0][0];

1128

}

1129

1130

}

1131

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

1132

derivatives[deriv_num] = new_coeff0_0*basisvalue0;

1133

}

1134

1135

// Transform derivatives back to physical element

1136

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

1137

{

1138

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

1139

{

1140

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

1141

}

1142

}

1143

// Delete pointer to array of derivatives on FIAT element

1144

delete [] derivatives;

1145

1146

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

1147

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

1148

{

1149

delete [] combinations[row];

1150

delete [] transform[row];

1151

}

1152

1153

delete [] combinations;

1154

delete [] transform;

1155

}

1156

1157

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

1158

virtual void evaluate_basis_derivatives_all(unsigned int n,

1159

double* values,

1160

const double* coordinates,

1161

const ufc::cell& c) const

1162

{

1163

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

1164

}

1165

1166

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

1167

virtual double evaluate_dof(unsigned int i,

1168

const ufc::function& f,

1169

const ufc::cell& c) const

1170

{

1171

// The reference points, direction and weights:

1172

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

1173

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

1174

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

1175

1176

const double * const * x = c.coordinates;

1177

double result = 0.0;

1178

// Iterate over the points:

1179

// Evaluate basis functions for affine mapping

1180

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

1181

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

1182

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

1183

1184

// Compute affine mapping y = F(X)

1185

double y[2];

1186

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

1187

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

1188

1189

// Evaluate function at physical points

1190

double values[1];

1191

f.evaluate(values, y, c);

1192

1193

// Map function values using appropriate mapping

1194

// Affine map: Do nothing

1195

1196

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

1197

1198

// Take directional components

1199

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

1200

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

1201

// Multiply by weights

1202

result *= W[i][0];

1203

1204

return result;

1205

}

1206

1207

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

1208

virtual void evaluate_dofs(double* values,

1209

const ufc::function& f,

1210

const ufc::cell& c) const

1211

{

1212

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

1213

}

1214

1215

/// Interpolate vertex values from dof values

1216

virtual void interpolate_vertex_values(double* vertex_values,

1217

const double* dof_values,

1218

const ufc::cell& c) const

1219

{

1220

// Evaluate at vertices and use affine mapping

1221

vertex_values[0] = dof_values[0];

1222

vertex_values[1] = dof_values[0];

1223

vertex_values[2] = dof_values[0];

1224

}

1225

1226

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

1227

virtual unsigned int num_sub_elements() const

1228

{

1229

return 1;

1230

}

1231

1232

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

1233

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

1234

{

1235

return new poissondg_0_finite_element_2();

1236

}

1237

1238

};

1239

1240

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

1241

/// degrees of freedom (dofs).

1242

1243

class poissondg_0_dof_map_0: public ufc::dof_map

1244

{

1245

private:

1246

1247

unsigned int __global_dimension;

1248

1249

public:

1250

1251

/// Constructor

1252

poissondg_0_dof_map_0() : ufc::dof_map()

1253

{

1254

__global_dimension = 0;

1255

}

1256

1257

/// Destructor

1258

virtual ~poissondg_0_dof_map_0()

1259

{

1260

// Do nothing

1261

}

1262

1263

/// Return a string identifying the dof map

1264

virtual const char* signature() const

1265

{

1266

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

1267

}

1268

1269

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

1270

virtual bool needs_mesh_entities(unsigned int d) const

1271

{

1272

switch ( d )

1273

{

1274

case 0:

1275

return false;

1276

break;

1277

case 1:

1278

return false;

1279

break;

1280

case 2:

1281

return true;

1282

break;

1283

}

1284

return false;

1285

}

1286

1287

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

1288

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

1289

{

1290

__global_dimension = 3*m.num_entities[2];

1291

return false;

1292

}

1293

1294

/// Initialize dof map for given cell

1295

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

1296

const ufc::cell& c)

1297

{

1298

// Do nothing

1299

}

1300

1301

/// Finish initialization of dof map for cells

1302

virtual void init_cell_finalize()

1303

{

1304

// Do nothing

1305

}

1306

1307

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

1308

virtual unsigned int global_dimension() const

1309

{

1310

return __global_dimension;

1311

}

1312

1313

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

1314

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

1315

{

1316

return 3;

1317

}

1318

1319

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

1320

virtual unsigned int max_local_dimension() const

1321

{

1322

return 3;

1323

}

1324

1325

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

1326

virtual unsigned int geometric_dimension() const

1327

{

1328

return 2;

1329

}

1330

1331

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

1332

virtual unsigned int num_facet_dofs() const

1333

{

1334

return 0;

1335

}

1336

1337

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

1338

virtual unsigned int num_entity_dofs(unsigned int d) const

1339

{

1340

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

1341

}

1342

1343

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

1344

virtual void tabulate_dofs(unsigned int* dofs,

1345

const ufc::mesh& m,

1346

const ufc::cell& c) const

1347

{

1348

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

1349

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

1350

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

1351

}

1352

1353

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

1354

virtual void tabulate_facet_dofs(unsigned int* dofs,

1355

unsigned int facet) const

1356

{

1357

switch ( facet )

1358

{

1359

case 0:

1360

1361

break;

1362

case 1:

1363

1364

break;

1365

case 2:

1366

1367

break;

1368

}

1369

}

1370

1371

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

1372

virtual void tabulate_entity_dofs(unsigned int* dofs,

1373

unsigned int d, unsigned int i) const

1374

{

1375

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

1376

}

1377

1378

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

1379

virtual void tabulate_coordinates(double** coordinates,

1380

const ufc::cell& c) const

1381

{

1382

const double * const * x = c.coordinates;

1383

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

1384

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

1385

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

1386

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

1387

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

1388

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

1389

}

1390

1391

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

1392

virtual unsigned int num_sub_dof_maps() const

1393

{

1394

return 1;

1395

}

1396

1397

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

1398

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

1399

{

1400

return new poissondg_0_dof_map_0();

1401

}

1402

1403

};

1404

1405

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

1406

/// degrees of freedom (dofs).

1407

1408

class poissondg_0_dof_map_1: public ufc::dof_map

1409

{

1410

private:

1411

1412

unsigned int __global_dimension;

1413

1414

public:

1415

1416

/// Constructor

1417

poissondg_0_dof_map_1() : ufc::dof_map()

1418

{

1419

__global_dimension = 0;

1420

}

1421

1422

/// Destructor

1423

virtual ~poissondg_0_dof_map_1()

1424

{

1425

// Do nothing

1426

}

1427

1428

/// Return a string identifying the dof map

1429

virtual const char* signature() const

1430

{

1431

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

1432

}

1433

1434

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

1435

virtual bool needs_mesh_entities(unsigned int d) const

1436

{

1437

switch ( d )

1438

{

1439

case 0:

1440

return false;

1441

break;

1442

case 1:

1443

return false;

1444

break;

1445

case 2:

1446

return true;

1447

break;

1448

}

1449

return false;

1450

}

1451

1452

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

1453

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

1454

{

1455

__global_dimension = 3*m.num_entities[2];

1456

return false;

1457

}

1458

1459

/// Initialize dof map for given cell

1460

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

1461

const ufc::cell& c)

1462

{

1463

// Do nothing

1464

}

1465

1466

/// Finish initialization of dof map for cells

1467

virtual void init_cell_finalize()

1468

{

1469

// Do nothing

1470

}

1471

1472

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

1473

virtual unsigned int global_dimension() const

1474

{

1475

return __global_dimension;

1476

}

1477

1478

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

1479

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

1480

{

1481

return 3;

1482

}

1483

1484

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

1485

virtual unsigned int max_local_dimension() const

1486

{

1487

return 3;

1488

}

1489

1490

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

1491

virtual unsigned int geometric_dimension() const

1492

{

1493

return 2;

1494

}

1495

1496

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

1497

virtual unsigned int num_facet_dofs() const

1498

{

1499

return 0;

1500

}

1501

1502

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

1503

virtual unsigned int num_entity_dofs(unsigned int d) const

1504

{

1505

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

1506

}

1507

1508

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

1509

virtual void tabulate_dofs(unsigned int* dofs,

1510

const ufc::mesh& m,

1511

const ufc::cell& c) const

1512

{

1513

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

1514

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

1515

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

1516

}

1517

1518

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

1519

virtual void tabulate_facet_dofs(unsigned int* dofs,

1520

unsigned int facet) const

1521

{

1522

switch ( facet )

1523

{

1524

case 0:

1525

1526

break;

1527

case 1:

1528

1529

break;

1530

case 2:

1531

1532

break;

1533

}

1534

}

1535

1536

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

1537

virtual void tabulate_entity_dofs(unsigned int* dofs,

1538

unsigned int d, unsigned int i) const

1539

{

1540

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

1541

}

1542

1543

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

1544

virtual void tabulate_coordinates(double** coordinates,

1545

const ufc::cell& c) const

1546

{

1547

const double * const * x = c.coordinates;

1548

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

1549

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

1550

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

1551

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

1552

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

1553

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

1554

}

1555

1556

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

1557

virtual unsigned int num_sub_dof_maps() const

1558

{

1559

return 1;

1560

}

1561

1562

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

1563

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

1564

{

1565

return new poissondg_0_dof_map_1();

1566

}

1567

1568

};

1569

1570

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

1571

/// degrees of freedom (dofs).

1572

1573

class poissondg_0_dof_map_2: public ufc::dof_map

1574

{

1575

private:

1576

1577

unsigned int __global_dimension;

1578

1579

public:

1580

1581

/// Constructor

1582

poissondg_0_dof_map_2() : ufc::dof_map()

1583

{

1584

__global_dimension = 0;

1585

}

1586

1587

/// Destructor

1588

virtual ~poissondg_0_dof_map_2()

1589

{

1590

// Do nothing

1591

}

1592

1593

/// Return a string identifying the dof map

1594

virtual const char* signature() const

1595

{

1596

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

1597

}

1598

1599

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

1600

virtual bool needs_mesh_entities(unsigned int d) const

1601

{

1602

switch ( d )

1603

{

1604

case 0:

1605

return false;

1606

break;

1607

case 1:

1608

return false;

1609

break;

1610

case 2:

1611

return true;

1612

break;

1613

}

1614

return false;

1615

}

1616

1617

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

1618

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

1619

{

1620

__global_dimension = m.num_entities[2];

1621

return false;

1622

}

1623

1624

/// Initialize dof map for given cell

1625

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

1626

const ufc::cell& c)

1627

{

1628

// Do nothing

1629

}

1630

1631

/// Finish initialization of dof map for cells

1632

virtual void init_cell_finalize()

1633

{

1634

// Do nothing

1635

}

1636

1637

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

1638

virtual unsigned int global_dimension() const

1639

{

1640

return __global_dimension;

1641

}

1642

1643

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

1644

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

1645

{

1646

return 1;

1647

}

1648

1649

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

1650

virtual unsigned int max_local_dimension() const

1651

{

1652

return 1;

1653

}

1654

1655

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

1656

virtual unsigned int geometric_dimension() const

1657

{

1658

return 2;

1659

}

1660

1661

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

1662

virtual unsigned int num_facet_dofs() const

1663

{

1664

return 0;

1665

}

1666

1667

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

1668

virtual unsigned int num_entity_dofs(unsigned int d) const

1669

{

1670

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

1671

}

1672

1673

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

1674

virtual void tabulate_dofs(unsigned int* dofs,

1675

const ufc::mesh& m,

1676

const ufc::cell& c) const

1677

{

1678

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

1679

}

1680

1681

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

1682

virtual void tabulate_facet_dofs(unsigned int* dofs,

1683

unsigned int facet) const

1684

{

1685

switch ( facet )

1686

{

1687

case 0:

1688

1689

break;

1690

case 1:

1691

1692

break;

1693

case 2:

1694

1695

break;

1696

}

1697

}

1698

1699

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

1700

virtual void tabulate_entity_dofs(unsigned int* dofs,

1701

unsigned int d, unsigned int i) const

1702

{

1703

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

1704

}

1705

1706

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

1707

virtual void tabulate_coordinates(double** coordinates,

1708

const ufc::cell& c) const

1709

{

1710

const double * const * x = c.coordinates;

1711

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

1712

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

1713

}

1714

1715

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

1716

virtual unsigned int num_sub_dof_maps() const

1717

{

1718

return 1;

1719

}

1720

1721

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

1722

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

1723

{

1724

return new poissondg_0_dof_map_2();

1725

}

1726

1727

};

1728

1729

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

1730

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

1731

/// the integral over a cell.

1732

1733

class poissondg_0_cell_integral_0_quadrature: public ufc::cell_integral

1734

{

1735

public:

1736

1737

/// Constructor

1738

poissondg_0_cell_integral_0_quadrature() : ufc::cell_integral()

1739

{

1740

// Do nothing

1741

}

1742

1743

/// Destructor

1744

virtual ~poissondg_0_cell_integral_0_quadrature()

1745

{

1746

// Do nothing

1747

}

1748

1749

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

1750

virtual void tabulate_tensor(double* A,

1751

const double * const * w,

1752

const ufc::cell& c) const

1753

{

1754

// Extract vertex coordinates

1755

const double * const * x = c.coordinates;

1756

1757

// Compute Jacobian of affine map from reference cell

1758

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

1759

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

1760

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

1761

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

1762

1763

// Compute determinant of Jacobian

1764

double detJ = J_00*J_11 - J_01*J_10;

1765

1766

// Compute inverse of Jacobian

1767

const double Jinv_00 = J_11 / detJ;

1768

const double Jinv_01 = -J_01 / detJ;

1769

const double Jinv_10 = -J_10 / detJ;

1770

const double Jinv_11 = J_00 / detJ;

1771

1772

// Set scale factor

1773

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

1774

1775

1776

// Array of quadrature weights

1777

static const double W1 = 0.5;

1778

// Quadrature points on the UFC reference element: (0.333333333, 0.333333333)

1779

1780

// Value of basis functions at quadrature points.

1781

static const double FE0_D01[1][3] = \

1782

{{-1, 0, 1}};

1783

1784

static const double FE0_D10[1][3] = \

1785

{{-1, 1, 0}};

1786

1787

1788

// Compute element tensor using UFL quadrature representation

1789

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

1790

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

1791

1792

// Loop quadrature points for integral

1793

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

1794

// Only 1 integration point, omitting IP loop.

1795

1796

// Number of operations for primary indices: 162

1797

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

1798

{

1799

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

1800

{

1801

// Number of operations to compute entry: 18

1802

A[j*3 + k] += ((Jinv_01*FE0_D10[0][j] + Jinv_11*FE0_D01[0][j])*(Jinv_01*FE0_D10[0][k] + Jinv_11*FE0_D01[0][k]) + (Jinv_00*FE0_D10[0][j] + Jinv_10*FE0_D01[0][j])*(Jinv_00*FE0_D10[0][k] + Jinv_10*FE0_D01[0][k]))*W1*det;

1803

}// end loop over 'k'

1804

}// end loop over 'j'

1805

}

1806

1807

};

1808

1809

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

1810

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

1811

/// the integral over a cell.

1812

1813

class poissondg_0_cell_integral_0: public ufc::cell_integral

1814

{

1815

private:

1816

1817

poissondg_0_cell_integral_0_quadrature integral_0_quadrature;

1818

1819

public:

1820

1821

/// Constructor

1822

poissondg_0_cell_integral_0() : ufc::cell_integral()

1823

{

1824

// Do nothing

1825

}

1826

1827

/// Destructor

1828

virtual ~poissondg_0_cell_integral_0()

1829

{

1830

// Do nothing

1831

}

1832

1833

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

1834

virtual void tabulate_tensor(double* A,

1835

const double * const * w,

1836

const ufc::cell& c) const

1837

{

1838

// Reset values of the element tensor block

1839

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

1840

A[j] = 0;

1841

1842

// Add all contributions to element tensor

1843

integral_0_quadrature.tabulate_tensor(A, w, c);

1844

}

1845

1846

};

1847

1848

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

1849

/// exterior facet tensor corresponding to the local contribution to

1850

/// a form from the integral over an exterior facet.

1851

1852

class poissondg_0_exterior_facet_integral_0_quadrature: public ufc::exterior_facet_integral

1853

{

1854

public:

1855

1856

/// Constructor

1857

poissondg_0_exterior_facet_integral_0_quadrature() : ufc::exterior_facet_integral()

1858

{

1859

// Do nothing

1860

}

1861

1862

/// Destructor

1863

virtual ~poissondg_0_exterior_facet_integral_0_quadrature()

1864

{

1865

// Do nothing

1866

}

1867

1868

/// Tabulate the tensor for the contribution from a local exterior facet

1869

virtual void tabulate_tensor(double* A,

1870

const double * const * w,

1871

const ufc::cell& c,

1872

unsigned int facet) const

1873

{

1874

// Extract vertex coordinates

1875

const double * const * x = c.coordinates;

1876

1877

// Compute Jacobian of affine map from reference cell

1878

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

1879

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

1880

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

1881

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

1882

1883

// Compute determinant of Jacobian

1884

double detJ = J_00*J_11 - J_01*J_10;

1885

1886

// Compute inverse of Jacobian

1887

const double Jinv_00 = J_11 / detJ;

1888

const double Jinv_01 = -J_01 / detJ;

1889

const double Jinv_10 = -J_10 / detJ;

1890

const double Jinv_11 = J_00 / detJ;

1891

1892

// Vertices on edges

1893

static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};

1894

1895

// Get vertices

1896

const unsigned int v0 = edge_vertices[facet][0];

1897

const unsigned int v1 = edge_vertices[facet][1];

1898

1899

// Compute scale factor (length of edge scaled by length of reference interval)

1900

const double dx0 = x[v1][0] - x[v0][0];

1901

const double dx1 = x[v1][1] - x[v0][1];

1902

const double det = std::sqrt(dx0*dx0 + dx1*dx1);

1903

1904

const bool direction = dx1*(x[facet][0] - x[v0][0]) - dx0*(x[facet][1] - x[v0][1]) < 0;

1905

1906

// Compute facet normals from the facet scale factor constants

1907

const double n0 = direction ? dx1 / det : -dx1 / det;

1908

const double n1 = direction ? -dx0 / det : dx0 / det;

1909

1910

1911

// Array of quadrature weights

1912

static const double W2[2] = {0.5, 0.5};

1913

// Quadrature points on the UFC reference element: (0.211324865), (0.788675135)

1914

1915

// Value of basis functions at quadrature points.

1916

static const double FE1_f0[2][3] = \

1917

{{0, 0.788675135, 0.211324865},

1918

{0, 0.211324865, 0.788675135}};

1919

1920

static const double FE1_f0_D01[2][3] = \

1921

{{-1, 0, 1},

1922

{-1, 0, 1}};

1923

1924

static const double FE1_f0_D10[2][3] = \

1925

{{-1, 1, 0},

1926

{-1, 1, 0}};

1927

1928

static const double FE1_f1[2][3] = \

1929

{{0.788675135, 0, 0.211324865},

1930

{0.211324865, 0, 0.788675135}};

1931

1932

static const double FE1_f2[2][3] = \

1933

{{0.788675135, 0.211324865, 0},

1934

{0.211324865, 0.788675135, 0}};

1935

1936

1937

// Compute element tensor using UFL quadrature representation

1938

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

1939

switch ( facet )

1940

{

1941

case 0:

1942

{

1943

// Total number of operations to compute element tensor (from this point): 558

1944

1945

// Loop quadrature points for integral

1946

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

1947

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

1948

{

1949

1950

// Number of operations for primary indices: 279

1951

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

1952

{

1953

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

1954

{

1955

// Number of operations to compute entry: 31

1956

A[j*3 + k] += (FE1_f0[ip][k]*FE1_f0[ip][j]*8/(w[0][0]) + ((FE1_f0[ip][j]*n0*(Jinv_00*FE1_f0_D10[ip][k] + Jinv_10*FE1_f0_D01[ip][k]) + FE1_f0[ip][j]*n1*(Jinv_01*FE1_f0_D10[ip][k] + Jinv_11*FE1_f0_D01[ip][k]))*-1 + (FE1_f0[ip][k]*n0*(Jinv_00*FE1_f0_D10[ip][j] + Jinv_10*FE1_f0_D01[ip][j]) + FE1_f0[ip][k]*n1*(Jinv_01*FE1_f0_D10[ip][j] + Jinv_11*FE1_f0_D01[ip][j]))*-1))*W2[ip]*det;

1957

}// end loop over 'k'

1958

}// end loop over 'j'

1959

}// end loop over 'ip'

1960

}

1961

break;

1962

case 1:

1963

{

1964

// Total number of operations to compute element tensor (from this point): 558

1965

1966

// Loop quadrature points for integral

1967

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

1968

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

1969

{

1970

1971

// Number of operations for primary indices: 279

1972

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

1973

{

1974

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

1975

{

1976

// Number of operations to compute entry: 31

1977

A[j*3 + k] += (FE1_f1[ip][k]*FE1_f1[ip][j]*8/(w[0][0]) + ((FE1_f1[ip][j]*n0*(Jinv_00*FE1_f0_D10[ip][k] + Jinv_10*FE1_f0_D01[ip][k]) + FE1_f1[ip][j]*n1*(Jinv_01*FE1_f0_D10[ip][k] + Jinv_11*FE1_f0_D01[ip][k]))*-1 + (FE1_f1[ip][k]*n1*(Jinv_01*FE1_f0_D10[ip][j] + Jinv_11*FE1_f0_D01[ip][j]) + FE1_f1[ip][k]*n0*(Jinv_00*FE1_f0_D10[ip][j] + Jinv_10*FE1_f0_D01[ip][j]))*-1))*W2[ip]*det;

1978

}// end loop over 'k'

1979

}// end loop over 'j'

1980

}// end loop over 'ip'

1981

}

1982

break;

1983

case 2:

1984

{

1985

// Total number of operations to compute element tensor (from this point): 558

1986

1987

// Loop quadrature points for integral

1988

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

1989

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

1990

{

1991

1992

// Number of operations for primary indices: 279

1993

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

1994

{

1995

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

1996

{

1997

// Number of operations to compute entry: 31

1998

A[j*3 + k] += (FE1_f2[ip][k]*FE1_f2[ip][j]*8/(w[0][0]) + ((FE1_f2[ip][j]*n1*(Jinv_01*FE1_f0_D10[ip][k] + Jinv_11*FE1_f0_D01[ip][k]) + FE1_f2[ip][j]*n0*(Jinv_00*FE1_f0_D10[ip][k] + Jinv_10*FE1_f0_D01[ip][k]))*-1 + (FE1_f2[ip][k]*n0*(Jinv_00*FE1_f0_D10[ip][j] + Jinv_10*FE1_f0_D01[ip][j]) + FE1_f2[ip][k]*n1*(Jinv_01*FE1_f0_D10[ip][j] + Jinv_11*FE1_f0_D01[ip][j]))*-1))*W2[ip]*det;

1999

}// end loop over 'k'

2000

}// end loop over 'j'

2001

}// end loop over 'ip'

2002

}

2003

break;

2004

}

2005

}

2006

2007

};

2008

2009

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

2010

/// exterior facet tensor corresponding to the local contribution to

2011

/// a form from the integral over an exterior facet.

2012

2013

class poissondg_0_exterior_facet_integral_0: public ufc::exterior_facet_integral

2014

{

2015

private:

2016

2017

poissondg_0_exterior_facet_integral_0_quadrature integral_0_quadrature;

2018

2019

public:

2020

2021

/// Constructor

2022

poissondg_0_exterior_facet_integral_0() : ufc::exterior_facet_integral()

2023

{

2024

// Do nothing

2025

}

2026

2027

/// Destructor

2028

virtual ~poissondg_0_exterior_facet_integral_0()

2029

{

2030

// Do nothing

2031

}

2032

2033

/// Tabulate the tensor for the contribution from a local exterior facet

2034

virtual void tabulate_tensor(double* A,

2035

const double * const * w,

2036

const ufc::cell& c,

2037

unsigned int facet) const

2038

{

2039

// Reset values of the element tensor block

2040

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

2041

A[j] = 0;

2042

2043

// Add all contributions to element tensor

2044

integral_0_quadrature.tabulate_tensor(A, w, c, facet);

2045

}

2046

2047

};

2048

2049

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

2050

/// interior facet tensor corresponding to the local contribution to

2051

/// a form from the integral over an interior facet.

2052

2053

class poissondg_0_interior_facet_integral_0_quadrature: public ufc::interior_facet_integral

2054

{

2055

public:

2056

2057

/// Constructor

2058

poissondg_0_interior_facet_integral_0_quadrature() : ufc::interior_facet_integral()

2059

{

2060

// Do nothing

2061

}

2062

2063

/// Destructor

2064

virtual ~poissondg_0_interior_facet_integral_0_quadrature()

2065

{

2066

// Do nothing

2067

}

2068

2069

/// Tabulate the tensor for the contribution from a local interior facet

2070

virtual void tabulate_tensor(double* A,

2071

const double * const * w,

2072

const ufc::cell& c0,

2073

const ufc::cell& c1,

2074

unsigned int facet0,

2075

unsigned int facet1) const

2076

{

2077

// Extract vertex coordinates

2078

const double * const * x0 = c0.coordinates;

2079

2080

// Compute Jacobian of affine map from reference cell

2081

const double J0_00 = x0[1][0] - x0[0][0];

2082

const double J0_01 = x0[2][0] - x0[0][0];

2083

const double J0_10 = x0[1][1] - x0[0][1];

2084

const double J0_11 = x0[2][1] - x0[0][1];

2085

2086

// Compute determinant of Jacobian

2087

double detJ0 = J0_00*J0_11 - J0_01*J0_10;

2088

2089

// Compute inverse of Jacobian

2090

const double Jinv0_00 = J0_11 / detJ0;

2091

const double Jinv0_01 = -J0_01 / detJ0;

2092

const double Jinv0_10 = -J0_10 / detJ0;

2093

const double Jinv0_11 = J0_00 / detJ0;

2094

2095

// Extract vertex coordinates

2096

const double * const * x1 = c1.coordinates;

2097

2098

// Compute Jacobian of affine map from reference cell

2099

const double J1_00 = x1[1][0] - x1[0][0];

2100

const double J1_01 = x1[2][0] - x1[0][0];

2101

const double J1_10 = x1[1][1] - x1[0][1];

2102

const double J1_11 = x1[2][1] - x1[0][1];

2103

2104

// Compute determinant of Jacobian

2105

double detJ1 = J1_00*J1_11 - J1_01*J1_10;

2106

2107

// Compute inverse of Jacobian

2108

const double Jinv1_00 = J1_11 / detJ1;

2109

const double Jinv1_01 = -J1_01 / detJ1;

2110

const double Jinv1_10 = -J1_10 / detJ1;

2111

const double Jinv1_11 = J1_00 / detJ1;

2112

2113

// Vertices on edges

2114

static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};

2115

2116

// Get vertices

2117

const unsigned int v0 = edge_vertices[facet0][0];

2118

const unsigned int v1 = edge_vertices[facet0][1];

2119

2120

// Compute scale factor (length of edge scaled by length of reference interval)

2121

const double dx0 = x0[v1][0] - x0[v0][0];

2122

const double dx1 = x0[v1][1] - x0[v0][1];

2123

const double det = std::sqrt(dx0*dx0 + dx1*dx1);

2124

2125

const bool direction = dx1*(x0[facet0][0] - x0[v0][0]) - dx0*(x0[facet0][1] - x0[v0][1]) < 0;

2126

2127

// Compute facet normals from the facet scale factor constants

2128

const double n00 = direction ? dx1 / det : -dx1 / det;

2129

const double n01 = direction ? -dx0 / det : dx0 / det;

2130

2131

// Compute facet normals from the facet scale factor constants

2132

const double n10 = !direction ? dx1 / det : -dx1 / det;

2133

const double n11 = !direction ? -dx0 / det : dx0 / det;

2134

2135

2136

// Array of quadrature weights

2137

static const double W2[2] = {0.5, 0.5};

2138

// Quadrature points on the UFC reference element: (0.211324865), (0.788675135)

2139

2140

// Value of basis functions at quadrature points.

2141

static const double FE1_f0[2][3] = \

2142

{{0, 0.788675135, 0.211324865},

2143

{0, 0.211324865, 0.788675135}};

2144

2145

static const double FE1_f0_D01[2][3] = \

2146

{{-1, 0, 1},

2147

{-1, 0, 1}};

2148

2149

static const double FE1_f0_D10[2][3] = \

2150

{{-1, 1, 0},

2151

{-1, 1, 0}};

2152

2153

static const double FE1_f1[2][3] = \

2154

{{0.788675135, 0, 0.211324865},

2155

{0.211324865, 0, 0.788675135}};

2156

2157

static const double FE1_f2[2][3] = \

2158

{{0.788675135, 0.211324865, 0},

2159

{0.211324865, 0.788675135, 0}};

2160

2161

2162

// Compute element tensor using UFL quadrature representation

2163

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

2164

switch ( facet0 )

2165

{

2166

case 0:

2167

switch ( facet1 )

2168

{

2169

case 0:

2170

{

2171

// Total number of operations to compute element tensor (from this point): 3024

2172

2173

// Loop quadrature points for integral

2174

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

2175

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

2176

{

2177

2178

// Number of operations for primary indices: 1512

2179

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

2180

{

2181

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

2182

{

2183

// Number of operations to compute entry: 42

2184

A[(j + 3)*6 + k] += ((FE1_f0[ip][j]*n10*FE1_f0[ip][k]*n00 + FE1_f0[ip][j]*n11*FE1_f0[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n00 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n01)*-1 + ((Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n10 + (Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n11)*-1))*W2[ip]*det;

2185

// Number of operations to compute entry: 42

2186

A[j*6 + (k + 3)] += ((((Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n10 + (Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n11)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n00 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n01)*-1) + (FE1_f0[ip][j]*n00*FE1_f0[ip][k]*n10 + FE1_f0[ip][j]*n01*FE1_f0[ip][k]*n11)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2187

// Number of operations to compute entry: 42

2188

A[j*6 + k] += ((FE1_f0[ip][j]*n00*FE1_f0[ip][k]*n00 + FE1_f0[ip][j]*n01*FE1_f0[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n01 + (Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n00)*-1 + ((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n01 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n00)*-1))*W2[ip]*det;

2189

// Number of operations to compute entry: 42

2190

A[(j + 3)*6 + (k + 3)] += ((FE1_f0[ip][j]*n11*FE1_f0[ip][k]*n11 + FE1_f0[ip][j]*n10*FE1_f0[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n10 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n11)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n10 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n11)*-1))*W2[ip]*det;

2191

}// end loop over 'k'

2192

}// end loop over 'j'

2193

}// end loop over 'ip'

2194

}

2195

break;

2196

case 1:

2197

{

2198

// Total number of operations to compute element tensor (from this point): 3024

2199

2200

// Loop quadrature points for integral

2201

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

2202

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

2203

{

2204

2205

// Number of operations for primary indices: 1512

2206

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

2207

{

2208

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

2209

{

2210

// Number of operations to compute entry: 42

2211

A[(j + 3)*6 + k] += ((FE1_f1[ip][j]*n10*FE1_f0[ip][k]*n00 + FE1_f1[ip][j]*n11*FE1_f0[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n11 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n10)*-1 + ((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n00 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n01)*-1))*W2[ip]*det;

2212

// Number of operations to compute entry: 42

2213

A[j*6 + (k + 3)] += ((FE1_f0[ip][j]*n00*FE1_f1[ip][k]*n10 + FE1_f0[ip][j]*n01*FE1_f1[ip][k]*n11)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n11 + (Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n10)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n00 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n01)*-1))*W2[ip]*det;

2214

// Number of operations to compute entry: 42

2215

A[j*6 + k] += ((FE1_f0[ip][j]*n00*FE1_f0[ip][k]*n00 + FE1_f0[ip][j]*n01*FE1_f0[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n01 + (Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n00)*-1 + ((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n01 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n00)*-1))*W2[ip]*det;

2216

// Number of operations to compute entry: 42

2217

A[(j + 3)*6 + (k + 3)] += ((FE1_f1[ip][j]*n10*FE1_f1[ip][k]*n10 + FE1_f1[ip][j]*n11*FE1_f1[ip][k]*n11)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n10 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n11)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n10 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n11)*-1))*W2[ip]*det;

2218

}// end loop over 'k'

2219

}// end loop over 'j'

2220

}// end loop over 'ip'

2221

}

2222

break;

2223

case 2:

2224

{

2225

// Total number of operations to compute element tensor (from this point): 3024

2226

2227

// Loop quadrature points for integral

2228

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

2229

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

2230

{

2231

2232

// Number of operations for primary indices: 1512

2233

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

2234

{

2235

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

2236

{

2237

// Number of operations to compute entry: 42

2238

A[(j + 3)*6 + k] += ((((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n00 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n01)*-1 + ((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n11 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n10)*-1) + (FE1_f2[ip][j]*n10*FE1_f0[ip][k]*n00 + FE1_f2[ip][j]*n11*FE1_f0[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2239

// Number of operations to compute entry: 42

2240

A[j*6 + (k + 3)] += ((((Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n11 + (Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n10)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n00 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n01)*-1) + (FE1_f0[ip][j]*n01*FE1_f2[ip][k]*n11 + FE1_f0[ip][j]*n00*FE1_f2[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2241

// Number of operations to compute entry: 42

2242

A[j*6 + k] += ((FE1_f0[ip][j]*n00*FE1_f0[ip][k]*n00 + FE1_f0[ip][j]*n01*FE1_f0[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n01 + (Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n00)*-1 + ((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n01 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n00)*-1))*W2[ip]*det;

2243

// Number of operations to compute entry: 42

2244

A[(j + 3)*6 + (k + 3)] += ((FE1_f2[ip][j]*n11*FE1_f2[ip][k]*n11 + FE1_f2[ip][j]*n10*FE1_f2[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n10 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n11)*-1 + ((Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n11 + (Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n10)*-1))*W2[ip]*det;

2245

}// end loop over 'k'

2246

}// end loop over 'j'

2247

}// end loop over 'ip'

2248

}

2249

break;

2250

}

2251

break;

2252

case 1:

2253

switch ( facet1 )

2254

{

2255

case 0:

2256

{

2257

// Total number of operations to compute element tensor (from this point): 3024

2258

2259

// Loop quadrature points for integral

2260

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

2261

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

2262

{

2263

2264

// Number of operations for primary indices: 1512

2265

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

2266

{

2267

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

2268

{

2269

// Number of operations to compute entry: 42

2270

A[(j + 3)*6 + k] += ((FE1_f0[ip][j]*n11*FE1_f1[ip][k]*n01 + FE1_f0[ip][j]*n10*FE1_f1[ip][k]*n00)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n00 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n01)*-1 + ((Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n10 + (Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n11)*-1))*W2[ip]*det;

2271

// Number of operations to compute entry: 42

2272

A[j*6 + (k + 3)] += ((((Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n10 + (Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n11)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n00 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n01)*-1) + (FE1_f1[ip][j]*n01*FE1_f0[ip][k]*n11 + FE1_f1[ip][j]*n00*FE1_f0[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2273

// Number of operations to compute entry: 42

2274

A[j*6 + k] += ((((Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n00 + (Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n01)*-1 + ((Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n00 + (Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n01)*-1) + (FE1_f1[ip][j]*n01*FE1_f1[ip][k]*n01 + FE1_f1[ip][j]*n00*FE1_f1[ip][k]*n00)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2275

// Number of operations to compute entry: 42

2276

A[(j + 3)*6 + (k + 3)] += ((FE1_f0[ip][j]*n11*FE1_f0[ip][k]*n11 + FE1_f0[ip][j]*n10*FE1_f0[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n10 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n11)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n10 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n11)*-1))*W2[ip]*det;

2277

}// end loop over 'k'

2278

}// end loop over 'j'

2279

}// end loop over 'ip'

2280

}

2281

break;

2282

case 1:

2283

{

2284

// Total number of operations to compute element tensor (from this point): 3024

2285

2286

// Loop quadrature points for integral

2287

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

2288

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

2289

{

2290

2291

// Number of operations for primary indices: 1512

2292

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

2293

{

2294

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

2295

{

2296

// Number of operations to compute entry: 42

2297

A[(j + 3)*6 + k] += ((((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n00 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n01)*-1 + ((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n11 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n10)*-1) + (FE1_f1[ip][j]*n10*FE1_f1[ip][k]*n00 + FE1_f1[ip][j]*n11*FE1_f1[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2298

// Number of operations to compute entry: 42

2299

A[j*6 + (k + 3)] += ((FE1_f1[ip][j]*n00*FE1_f1[ip][k]*n10 + FE1_f1[ip][j]*n01*FE1_f1[ip][k]*n11)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n11 + (Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n10)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n00 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n01)*-1))*W2[ip]*det;

2300

// Number of operations to compute entry: 42

2301

A[j*6 + k] += ((((Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n00 + (Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n01)*-1 + ((Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n00 + (Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n01)*-1) + (FE1_f1[ip][j]*n01*FE1_f1[ip][k]*n01 + FE1_f1[ip][j]*n00*FE1_f1[ip][k]*n00)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2302

// Number of operations to compute entry: 42

2303

A[(j + 3)*6 + (k + 3)] += ((FE1_f1[ip][j]*n10*FE1_f1[ip][k]*n10 + FE1_f1[ip][j]*n11*FE1_f1[ip][k]*n11)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n10 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n11)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n10 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n11)*-1))*W2[ip]*det;

2304

}// end loop over 'k'

2305

}// end loop over 'j'

2306

}// end loop over 'ip'

2307

}

2308

break;

2309

case 2:

2310

{

2311

// Total number of operations to compute element tensor (from this point): 3024

2312

2313

// Loop quadrature points for integral

2314

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

2315

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

2316

{

2317

2318

// Number of operations for primary indices: 1512

2319

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

2320

{

2321

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

2322

{

2323

// Number of operations to compute entry: 42

2324

A[(j + 3)*6 + k] += ((FE1_f2[ip][j]*n10*FE1_f1[ip][k]*n00 + FE1_f2[ip][j]*n11*FE1_f1[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n00 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n01)*-1 + ((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n11 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n10)*-1))*W2[ip]*det;

2325

// Number of operations to compute entry: 42

2326

A[j*6 + (k + 3)] += ((((Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n11 + (Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n10)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n00 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n01)*-1) + (FE1_f1[ip][j]*n00*FE1_f2[ip][k]*n10 + FE1_f1[ip][j]*n01*FE1_f2[ip][k]*n11)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2327

// Number of operations to compute entry: 42

2328

A[j*6 + k] += ((((Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n00 + (Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n01)*-1 + ((Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n00 + (Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n01)*-1) + (FE1_f1[ip][j]*n01*FE1_f1[ip][k]*n01 + FE1_f1[ip][j]*n00*FE1_f1[ip][k]*n00)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2329

// Number of operations to compute entry: 42

2330

A[(j + 3)*6 + (k + 3)] += ((FE1_f2[ip][j]*n11*FE1_f2[ip][k]*n11 + FE1_f2[ip][j]*n10*FE1_f2[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n10 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n11)*-1 + ((Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n11 + (Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n10)*-1))*W2[ip]*det;

2331

}// end loop over 'k'

2332

}// end loop over 'j'

2333

}// end loop over 'ip'

2334

}

2335

break;

2336

}

2337

break;

2338

case 2:

2339

switch ( facet1 )

2340

{

2341

case 0:

2342

{

2343

// Total number of operations to compute element tensor (from this point): 3024

2344

2345

// Loop quadrature points for integral

2346

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

2347

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

2348

{

2349

2350

// Number of operations for primary indices: 1512

2351

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

2352

{

2353

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

2354

{

2355

// Number of operations to compute entry: 42

2356

A[(j + 3)*6 + k] += ((FE1_f0[ip][j]*n11*FE1_f2[ip][k]*n01 + FE1_f0[ip][j]*n10*FE1_f2[ip][k]*n00)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n01 + (Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n00)*-1 + ((Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n10 + (Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n11)*-1))*W2[ip]*det;

2357

// Number of operations to compute entry: 42

2358

A[j*6 + (k + 3)] += ((FE1_f2[ip][j]*n01*FE1_f0[ip][k]*n11 + FE1_f2[ip][j]*n00*FE1_f0[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n10 + (Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n11)*-1 + ((Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n01 + (Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n00)*-1))*W2[ip]*det;

2359

// Number of operations to compute entry: 42

2360

A[j*6 + k] += ((((Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n00 + (Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n01)*-1 + ((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n01 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n00)*-1) + (FE1_f2[ip][j]*n00*FE1_f2[ip][k]*n00 + FE1_f2[ip][j]*n01*FE1_f2[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2361

// Number of operations to compute entry: 42

2362

A[(j + 3)*6 + (k + 3)] += ((FE1_f0[ip][j]*n11*FE1_f0[ip][k]*n11 + FE1_f0[ip][j]*n10*FE1_f0[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n10 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f0[ip][k]*n11)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n10 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f0[ip][j]*n11)*-1))*W2[ip]*det;

2363

}// end loop over 'k'

2364

}// end loop over 'j'

2365

}// end loop over 'ip'

2366

}

2367

break;

2368

case 1:

2369

{

2370

// Total number of operations to compute element tensor (from this point): 3024

2371

2372

// Loop quadrature points for integral

2373

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

2374

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

2375

{

2376

2377

// Number of operations for primary indices: 1512

2378

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

2379

{

2380

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

2381

{

2382

// Number of operations to compute entry: 42

2383

A[(j + 3)*6 + k] += ((FE1_f1[ip][j]*n11*FE1_f2[ip][k]*n01 + FE1_f1[ip][j]*n10*FE1_f2[ip][k]*n00)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n11 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n10)*-1 + ((Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n01 + (Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n00)*-1))*W2[ip]*det;

2384

// Number of operations to compute entry: 42

2385

A[j*6 + (k + 3)] += ((FE1_f2[ip][j]*n01*FE1_f1[ip][k]*n11 + FE1_f2[ip][j]*n00*FE1_f1[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n11 + (Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n10)*-1 + ((Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n01 + (Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n00)*-1))*W2[ip]*det;

2386

// Number of operations to compute entry: 42

2387

A[j*6 + k] += ((((Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n00 + (Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n01)*-1 + ((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n01 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n00)*-1) + (FE1_f2[ip][j]*n00*FE1_f2[ip][k]*n00 + FE1_f2[ip][j]*n01*FE1_f2[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2388

// Number of operations to compute entry: 42

2389

A[(j + 3)*6 + (k + 3)] += ((FE1_f1[ip][j]*n10*FE1_f1[ip][k]*n10 + FE1_f1[ip][j]*n11*FE1_f1[ip][k]*n11)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n10 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f1[ip][k]*n11)*-1 + ((Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n10 + (Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f1[ip][j]*n11)*-1))*W2[ip]*det;

2390

}// end loop over 'k'

2391

}// end loop over 'j'

2392

}// end loop over 'ip'

2393

}

2394

break;

2395

case 2:

2396

{

2397

// Total number of operations to compute element tensor (from this point): 3024

2398

2399

// Loop quadrature points for integral

2400

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

2401

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

2402

{

2403

2404

// Number of operations for primary indices: 1512

2405

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

2406

{

2407

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

2408

{

2409

// Number of operations to compute entry: 42

2410

A[(j + 3)*6 + k] += ((((Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n01 + (Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n00)*-1 + ((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n11 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n10)*-1) + (FE1_f2[ip][j]*n10*FE1_f2[ip][k]*n00 + FE1_f2[ip][j]*n11*FE1_f2[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2411

// Number of operations to compute entry: 42

2412

A[j*6 + (k + 3)] += ((FE1_f2[ip][j]*n01*FE1_f2[ip][k]*n11 + FE1_f2[ip][j]*n00*FE1_f2[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n11 + (Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n10)*-1 + ((Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n01 + (Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n00)*-1))*W2[ip]*det;

2413

// Number of operations to compute entry: 42

2414

A[j*6 + k] += ((((Jinv0_00*FE1_f0_D10[ip][j] + Jinv0_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n00 + (Jinv0_01*FE1_f0_D10[ip][j] + Jinv0_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n01)*-1 + ((Jinv0_01*FE1_f0_D10[ip][k] + Jinv0_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n01 + (Jinv0_00*FE1_f0_D10[ip][k] + Jinv0_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n00)*-1) + (FE1_f2[ip][j]*n00*FE1_f2[ip][k]*n00 + FE1_f2[ip][j]*n01*FE1_f2[ip][k]*n01)*4/((w[0][0] + w[0][1])/(2)))*W2[ip]*det;

2415

// Number of operations to compute entry: 42

2416

A[(j + 3)*6 + (k + 3)] += ((FE1_f2[ip][j]*n11*FE1_f2[ip][k]*n11 + FE1_f2[ip][j]*n10*FE1_f2[ip][k]*n10)*4/((w[0][0] + w[0][1])/(2)) + (((Jinv1_00*FE1_f0_D10[ip][j] + Jinv1_10*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n10 + (Jinv1_01*FE1_f0_D10[ip][j] + Jinv1_11*FE1_f0_D01[ip][j])*0.5*FE1_f2[ip][k]*n11)*-1 + ((Jinv1_01*FE1_f0_D10[ip][k] + Jinv1_11*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n11 + (Jinv1_00*FE1_f0_D10[ip][k] + Jinv1_10*FE1_f0_D01[ip][k])*0.5*FE1_f2[ip][j]*n10)*-1))*W2[ip]*det;

2417

}// end loop over 'k'

2418

}// end loop over 'j'

2419

}// end loop over 'ip'

2420

}

2421

break;

2422

}

2423

break;

2424

}

2425

}

2426

2427

};

2428

2429

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

2430

/// interior facet tensor corresponding to the local contribution to

2431

/// a form from the integral over an interior facet.

2432

2433

class poissondg_0_interior_facet_integral_0: public ufc::interior_facet_integral

2434

{

2435

private:

2436

2437

poissondg_0_interior_facet_integral_0_quadrature integral_0_quadrature;

2438

2439

public:

2440

2441

/// Constructor

2442

poissondg_0_interior_facet_integral_0() : ufc::interior_facet_integral()

2443

{

2444

// Do nothing

2445

}

2446

2447

/// Destructor

2448

virtual ~poissondg_0_interior_facet_integral_0()

2449

{

2450

// Do nothing

2451

}

2452

2453

/// Tabulate the tensor for the contribution from a local interior facet

2454

virtual void tabulate_tensor(double* A,

2455

const double * const * w,

2456

const ufc::cell& c0,

2457

const ufc::cell& c1,

2458

unsigned int facet0,

2459

unsigned int facet1) const

2460

{

2461

// Reset values of the element tensor block

2462

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

2463

A[j] = 0;

2464

2465

// Add all contributions to element tensor

2466

integral_0_quadrature.tabulate_tensor(A, w, c0, c1, facet0, facet1);

2467

}

2468

2469

};

2470

2471

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

2472

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

2473

/// mapping

2474

///

2475

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

2476

///

2477

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

2478

/// global tensor A is defined by

2479

///

2480

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

2481

///

2482

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

2483

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

2484

/// fixed functions (coefficients).

2485

2486

class poissondg_form_0: public ufc::form

2487

{

2488

public:

2489

2490

/// Constructor

2491

poissondg_form_0() : ufc::form()

2492

{

2493

// Do nothing

2494

}

2495

2496

/// Destructor

2497

virtual ~poissondg_form_0()

2498

{

2499

// Do nothing

2500

}

2501

2502

/// Return a string identifying the form

2503

virtual const char* signature() const

2504

{

2505

return "Form([Integral(IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})), Measure('cell', 0, None)), Integral(Sum(Product(Division(FloatValue(4.0, (), (), {}), Division(Sum(NegativeRestricted(Constant(Cell('triangle', 1, Space(2)), 0)), PositiveRestricted(Constant(Cell('triangle', 1, Space(2)), 0))), FloatValue(2.0, (), (), {}))), IndexSum(Product(Indexed(Sum(ComponentTensor(Product(Indexed(NegativeRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(3),), {Index(3): 2})), NegativeRestricted(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0))), MultiIndex((Index(3),), {Index(3): 2})), ComponentTensor(Product(Indexed(PositiveRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(4),), {Index(4): 2})), PositiveRestricted(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0))), MultiIndex((Index(4),), {Index(4): 2}))), MultiIndex((Index(5),), {Index(5): 2})), Indexed(Sum(ComponentTensor(Product(Indexed(NegativeRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(6),), {Index(6): 2})), NegativeRestricted(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1))), MultiIndex((Index(6),), {Index(6): 2})), ComponentTensor(Product(Indexed(PositiveRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(7),), {Index(7): 2})), PositiveRestricted(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1))), MultiIndex((Index(7),), {Index(7): 2}))), MultiIndex((Index(5),), {Index(5): 2}))), MultiIndex((Index(5),), {Index(5): 2}))), Sum(Product(IntValue(-1, (), (), {}), IndexSum(Product(Indexed(ComponentTensor(Product(FloatValue(0.5, (), (), {}), Indexed(Sum(NegativeRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0), MultiIndex((Index(8),), {Index(8): 2})), MultiIndex((Index(8),), {Index(8): 2}))), PositiveRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0), MultiIndex((Index(9),), {Index(9): 2})), MultiIndex((Index(9),), {Index(9): 2})))), MultiIndex((Index(10),), {Index(10): 2}))), MultiIndex((Index(10),), {Index(10): 2})), MultiIndex((Index(11),), {Index(11): 2})), Indexed(Sum(ComponentTensor(Product(Indexed(NegativeRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(12),), {Index(12): 2})), NegativeRestricted(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1))), MultiIndex((Index(12),), {Index(12): 2})), ComponentTensor(Product(Indexed(PositiveRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(13),), {Index(13): 2})), PositiveRestricted(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1))), MultiIndex((Index(13),), {Index(13): 2}))), MultiIndex((Index(11),), {Index(11): 2}))), MultiIndex((Index(11),), {Index(11): 2}))), Product(IntValue(-1, (), (), {}), IndexSum(Product(Indexed(ComponentTensor(Product(FloatValue(0.5, (), (), {}), Indexed(Sum(NegativeRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1), MultiIndex((Index(14),), {Index(14): 2})), MultiIndex((Index(14),), {Index(14): 2}))), PositiveRestricted(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1), MultiIndex((Index(15),), {Index(15): 2})), MultiIndex((Index(15),), {Index(15): 2})))), MultiIndex((Index(16),), {Index(16): 2}))), MultiIndex((Index(16),), {Index(16): 2})), MultiIndex((Index(17),), {Index(17): 2})), Indexed(Sum(ComponentTensor(Product(Indexed(NegativeRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(18),), {Index(18): 2})), NegativeRestricted(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0))), MultiIndex((Index(18),), {Index(18): 2})), ComponentTensor(Product(Indexed(PositiveRestricted(FacetNormal(Cell('triangle', 1, Space(2)))), MultiIndex((Index(19),), {Index(19): 2})), PositiveRestricted(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0))), MultiIndex((Index(19),), {Index(19): 2}))), MultiIndex((Index(17),), {Index(17): 2}))), MultiIndex((Index(17),), {Index(17): 2}))))), Measure('interior_facet', 0, None)), Integral(Sum(Product(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1), Product(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0), Division(FloatValue(8.0, (), (), {}), Constant(Cell('triangle', 1, Space(2)), 0)))), Sum(Product(IntValue(-1, (), (), {}), IndexSum(Product(Indexed(ComponentTensor(Product(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0), Indexed(FacetNormal(Cell('triangle', 1, Space(2))), MultiIndex((Index(20),), {Index(20): 2}))), MultiIndex((Index(20),), {Index(20): 2})), MultiIndex((Index(21),), {Index(21): 2})), Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1), MultiIndex((Index(22),), {Index(22): 2})), MultiIndex((Index(22),), {Index(22): 2})), MultiIndex((Index(21),), {Index(21): 2}))), MultiIndex((Index(21),), {Index(21): 2}))), Product(IntValue(-1, (), (), {}), IndexSum(Product(Indexed(ComponentTensor(Product(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1), Indexed(FacetNormal(Cell('triangle', 1, Space(2))), MultiIndex((Index(23),), {Index(23): 2}))), MultiIndex((Index(23),), {Index(23): 2})), MultiIndex((Index(24),), {Index(24): 2})), Indexed(ComponentTensor(SpatialDerivative(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0), MultiIndex((Index(25),), {Index(25): 2})), MultiIndex((Index(25),), {Index(25): 2})), MultiIndex((Index(24),), {Index(24): 2}))), MultiIndex((Index(24),), {Index(24): 2}))))), Measure('exterior_facet', 0, None))])";

2506

}

2507

2508

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

2509

virtual unsigned int rank() const

2510

{

2511

return 2;

2512

}

2513

2514

/// Return the number of coefficients (n)

2515

virtual unsigned int num_coefficients() const

2516

{

2517

return 1;

2518

}

2519

2520

/// Return the number of cell integrals

2521

virtual unsigned int num_cell_integrals() const

2522

{

2523

return 1;

2524

}

2525

2526

/// Return the number of exterior facet integrals

2527

virtual unsigned int num_exterior_facet_integrals() const

2528

{

2529

return 1;

2530

}

2531

2532

/// Return the number of interior facet integrals

2533

virtual unsigned int num_interior_facet_integrals() const

2534

{

2535

return 1;

2536

}

2537

2538

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

2539

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

2540

{

2541

switch ( i )

2542

{

2543

case 0:

2544

return new poissondg_0_finite_element_0();

2545

break;

2546

case 1:

2547

return new poissondg_0_finite_element_1();

2548

break;

2549

case 2:

2550

return new poissondg_0_finite_element_2();

2551

break;

2552

}

2553

return 0;

2554

}

2555

2556

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

2557

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

2558

{

2559

switch ( i )

2560

{

2561

case 0:

2562

return new poissondg_0_dof_map_0();

2563

break;

2564

case 1:

2565

return new poissondg_0_dof_map_1();

2566

break;

2567

case 2:

2568

return new poissondg_0_dof_map_2();

2569

break;

2570

}

2571

return 0;

2572

}

2573

2574

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

2575

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

2576

{

2577

return new poissondg_0_cell_integral_0();

2578

}

2579

2580

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

2581

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

2582

{

2583

return new poissondg_0_exterior_facet_integral_0();

2584

}

2585

2586

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

2587

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

2588

{

2589

return new poissondg_0_interior_facet_integral_0();

2590

}

2591

2592

};

2593

2594

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

2595

2596

class poissondg_1_finite_element_0: public ufc::finite_element

2597

{

2598

public:

2599

2600

/// Constructor

2601

poissondg_1_finite_element_0() : ufc::finite_element()

2602

{

2603

// Do nothing

2604

}

2605

2606

/// Destructor

2607

virtual ~poissondg_1_finite_element_0()

2608

{

2609

// Do nothing

2610

}

2611

2612

/// Return a string identifying the finite element

2613

virtual const char* signature() const

2614

{

2615

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

2616

}

2617

2618

/// Return the cell shape

2619

virtual ufc::shape cell_shape() const

2620

{

2621

return ufc::triangle;

2622

}

2623

2624

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

2625

virtual unsigned int space_dimension() const

2626

{

2627

return 3;

2628

}

2629

2630

/// Return the rank of the value space

2631

virtual unsigned int value_rank() const

2632

{

2633

return 0;

2634

}

2635

2636

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

2637

virtual unsigned int value_dimension(unsigned int i) const

2638

{

2639

return 1;

2640

}

2641

2642

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

2643

virtual void evaluate_basis(unsigned int i,

2644

double* values,

2645

const double* coordinates,

2646

const ufc::cell& c) const

2647

{

2648

// Extract vertex coordinates

2649

const double * const * element_coordinates = c.coordinates;

2650

2651

// Compute Jacobian of affine map from reference cell

2652

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

2653

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

2654

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

2655

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

2656

2657

// Compute determinant of Jacobian

2658

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

2659

2660

// Compute inverse of Jacobian

2661

2662

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

2663

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

2664

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

2665

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

2666

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

2667

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

2668

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

2669

2670

// Map coordinates to the reference square

2671

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

2672

x = -1.0;

2673

else

2674

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

2675

y = 2.0*y - 1.0;

2676

2677

// Reset values

2678

*values = 0;

2679

2680

// Map degree of freedom to element degree of freedom

2681

const unsigned int dof = i;

2682

2683

// Generate scalings

2684

const double scalings_y_0 = 1;

2685

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

2686

2687

// Compute psitilde_a

2688

const double psitilde_a_0 = 1;

2689

const double psitilde_a_1 = x;

2690

2691

// Compute psitilde_bs

2692

const double psitilde_bs_0_0 = 1;

2693

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2694

const double psitilde_bs_1_0 = 1;

2695

2696

// Compute basisvalues

2697

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

2698

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

2699

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2700

2701

// Table(s) of coefficients

2702

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

2703

{{0.471404521, -0.288675135, -0.166666667},

2704

{0.471404521, 0.288675135, -0.166666667},

2705

{0.471404521, 0, 0.333333333}};

2706

2707

// Extract relevant coefficients

2708

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

2709

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

2710

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

2711

2712

// Compute value(s)

2713

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;

2714

}

2715

2716

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

2717

virtual void evaluate_basis_all(double* values,

2718

const double* coordinates,

2719

const ufc::cell& c) const

2720

{

2721

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

2722

}

2723

2724

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

2725

virtual void evaluate_basis_derivatives(unsigned int i,

2726

unsigned int n,

2727

double* values,

2728

const double* coordinates,

2729

const ufc::cell& c) const

2730

{

2731

// Extract vertex coordinates

2732

const double * const * element_coordinates = c.coordinates;

2733

2734

// Compute Jacobian of affine map from reference cell

2735

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

2736

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

2737

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

2738

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

2739

2740

// Compute determinant of Jacobian

2741

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

2742

2743

// Compute inverse of Jacobian

2744

2745

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

2746

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

2747

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

2748

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

2749

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

2750

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

2751

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

2752

2753

// Map coordinates to the reference square

2754

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

2755

x = -1.0;

2756

else

2757

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

2758

y = 2.0*y - 1.0;

2759

2760

// Compute number of derivatives

2761

unsigned int num_derivatives = 1;

2762

2763

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

2764

num_derivatives *= 2;

2765

2766

2767

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

2768

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

2769

2770

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

2771

{

2772

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

2773

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

2774

combinations[j][k] = 0;

2775

}

2776

2777

// Generate combinations of derivatives

2778

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

2779

{

2780

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

2781

{

2782

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

2783

{

2784

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

2785

combinations[row][col] = 0;

2786

else

2787

{

2788

combinations[row][col] += 1;

2789

break;

2790

}

2791

}

2792

}

2793

}

2794

2795

// Compute inverse of Jacobian

2796

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

2797

2798

// Declare transformation matrix

2799

// Declare pointer to two dimensional array and initialise

2800

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

2801

2802

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

2803

{

2804

transform[j] = new double [num_derivatives];

2805

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

2806

transform[j][k] = 1;

2807

}

2808

2809

// Construct transformation matrix

2810

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

2811

{

2812

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

2813

{

2814

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

2815

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

2816

}

2817

}

2818

2819

// Reset values

2820

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

2821

values[j] = 0;

2822

2823

// Map degree of freedom to element degree of freedom

2824

const unsigned int dof = i;

2825

2826

// Generate scalings

2827

const double scalings_y_0 = 1;

2828

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

2829

2830

// Compute psitilde_a

2831

const double psitilde_a_0 = 1;

2832

const double psitilde_a_1 = x;

2833

2834

// Compute psitilde_bs

2835

const double psitilde_bs_0_0 = 1;

2836

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2837

const double psitilde_bs_1_0 = 1;

2838

2839

// Compute basisvalues

2840

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

2841

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

2842

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2843

2844

// Table(s) of coefficients

2845

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

2846

{{0.471404521, -0.288675135, -0.166666667},

2847

{0.471404521, 0.288675135, -0.166666667},

2848

{0.471404521, 0, 0.333333333}};

2849

2850

// Interesting (new) part

2851

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

2852

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

2853

{{0, 0, 0},

2854

{4.89897949, 0, 0},

2855

{0, 0, 0}};

2856

2857

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

2858

{{0, 0, 0},

2859

{2.44948974, 0, 0},

2860

{4.24264069, 0, 0}};

2861

2862

// Compute reference derivatives

2863

// Declare pointer to array of derivatives on FIAT element

2864

double *derivatives = new double [num_derivatives];

2865

2866

// Declare coefficients

2867

double coeff0_0 = 0;

2868

double coeff0_1 = 0;

2869

double coeff0_2 = 0;

2870

2871

// Declare new coefficients

2872

double new_coeff0_0 = 0;

2873

double new_coeff0_1 = 0;

2874

double new_coeff0_2 = 0;

2875

2876

// Loop possible derivatives

2877

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

2878

{

2879

// Get values from coefficients array

2880

new_coeff0_0 = coefficients0[dof][0];

2881

new_coeff0_1 = coefficients0[dof][1];

2882

new_coeff0_2 = coefficients0[dof][2];

2883

2884

// Loop derivative order

2885

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

2886

{

2887

// Update old coefficients

2888

coeff0_0 = new_coeff0_0;

2889

coeff0_1 = new_coeff0_1;

2890

coeff0_2 = new_coeff0_2;

2891

2892

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

2893

{

2894

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];

2895

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];

2896

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];

2897

}

2898

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

2899

{

2900

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];

2901

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];

2902

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];

2903

}

2904

2905

}

2906

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

2907

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;

2908

}

2909

2910

// Transform derivatives back to physical element

2911

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

2912

{

2913

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

2914

{

2915

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

2916

}

2917

}

2918

// Delete pointer to array of derivatives on FIAT element

2919

delete [] derivatives;

2920

2921

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

2922

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

2923

{

2924

delete [] combinations[row];

2925

delete [] transform[row];

2926

}

2927

2928

delete [] combinations;

2929

delete [] transform;

2930

}

2931

2932

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

2933

virtual void evaluate_basis_derivatives_all(unsigned int n,

2934

double* values,

2935

const double* coordinates,

2936

const ufc::cell& c) const

2937

{

2938

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

2939

}

2940

2941

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

2942

virtual double evaluate_dof(unsigned int i,

2943

const ufc::function& f,

2944

const ufc::cell& c) const

2945

{

2946

// The reference points, direction and weights:

2947

static const double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};

2948

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

2949

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

2950

2951

const double * const * x = c.coordinates;

2952

double result = 0.0;

2953

// Iterate over the points:

2954

// Evaluate basis functions for affine mapping

2955

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

2956

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

2957

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

2958

2959

// Compute affine mapping y = F(X)

2960

double y[2];

2961

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

2962

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

2963

2964

// Evaluate function at physical points

2965

double values[1];

2966

f.evaluate(values, y, c);

2967

2968

// Map function values using appropriate mapping

2969

// Affine map: Do nothing

2970

2971

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

2972

2973

// Take directional components

2974

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

2975

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

2976

// Multiply by weights

2977

result *= W[i][0];

2978

2979

return result;

2980

}

2981

2982

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

2983

virtual void evaluate_dofs(double* values,

2984

const ufc::function& f,

2985

const ufc::cell& c) const

2986

{

2987

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

2988

}

2989

2990

/// Interpolate vertex values from dof values

2991

virtual void interpolate_vertex_values(double* vertex_values,

2992

const double* dof_values,

2993

const ufc::cell& c) const

2994

{

2995

// Evaluate at vertices and use affine mapping

2996

vertex_values[0] = dof_values[0];

2997

vertex_values[1] = dof_values[1];

2998

vertex_values[2] = dof_values[2];

2999

}

3000

3001

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

3002

virtual unsigned int num_sub_elements() const

3003

{

3004

return 1;

3005

}

3006

3007

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

3008

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

3009

{

3010

return new poissondg_1_finite_element_0();

3011

}

3012

3013

};

3014

3015

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

3016

3017

class poissondg_1_finite_element_1: public ufc::finite_element

3018

{

3019

public:

3020

3021

/// Constructor

3022

poissondg_1_finite_element_1() : ufc::finite_element()

3023

{

3024

// Do nothing

3025

}

3026

3027

/// Destructor

3028

virtual ~poissondg_1_finite_element_1()

3029

{

3030

// Do nothing

3031

}

3032

3033

/// Return a string identifying the finite element

3034

virtual const char* signature() const

3035

{

3036

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

3037

}

3038

3039

/// Return the cell shape

3040

virtual ufc::shape cell_shape() const

3041

{

3042

return ufc::triangle;

3043

}

3044

3045

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

3046

virtual unsigned int space_dimension() const

3047

{

3048

return 3;

3049

}

3050

3051

/// Return the rank of the value space

3052

virtual unsigned int value_rank() const

3053

{

3054

return 0;

3055

}

3056

3057

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

3058

virtual unsigned int value_dimension(unsigned int i) const

3059

{

3060

return 1;

3061

}

3062

3063

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

3064

virtual void evaluate_basis(unsigned int i,

3065

double* values,

3066

const double* coordinates,

3067

const ufc::cell& c) const

3068

{

3069

// Extract vertex coordinates

3070

const double * const * element_coordinates = c.coordinates;

3071

3072

// Compute Jacobian of affine map from reference cell

3073

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

3074

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

3075

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

3076

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

3077

3078

// Compute determinant of Jacobian

3079

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

3080

3081

// Compute inverse of Jacobian

3082

3083

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

3084

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

3085

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

3086

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

3087

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

3088

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

3089

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

3090

3091

// Map coordinates to the reference square

3092

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

3093

x = -1.0;

3094

else

3095

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

3096

y = 2.0*y - 1.0;

3097

3098

// Reset values

3099

*values = 0;

3100

3101

// Map degree of freedom to element degree of freedom

3102

const unsigned int dof = i;

3103

3104

// Generate scalings

3105

const double scalings_y_0 = 1;

3106

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

3107

3108

// Compute psitilde_a

3109

const double psitilde_a_0 = 1;

3110

const double psitilde_a_1 = x;

3111

3112

// Compute psitilde_bs

3113

const double psitilde_bs_0_0 = 1;

3114

const double psitilde_bs_0_1 = 1.5*y + 0.5;

3115

const double psitilde_bs_1_0 = 1;

3116

3117

// Compute basisvalues

3118

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

3119

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

3120

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

3121

3122

// Table(s) of coefficients

3123

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

3124

{{0.471404521, -0.288675135, -0.166666667},

3125

{0.471404521, 0.288675135, -0.166666667},

3126

{0.471404521, 0, 0.333333333}};

3127

3128

// Extract relevant coefficients

3129

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

3130

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

3131

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

3132

3133

// Compute value(s)

3134

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;

3135

}

3136

3137

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

3138

virtual void evaluate_basis_all(double* values,

3139

const double* coordinates,

3140

const ufc::cell& c) const

3141

{

3142

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

3143

}

3144

3145

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

3146

virtual void evaluate_basis_derivatives(unsigned int i,

3147

unsigned int n,

3148

double* values,

3149

const double* coordinates,

3150

const ufc::cell& c) const

3151

{

3152

// Extract vertex coordinates

3153

const double * const * element_coordinates = c.coordinates;

3154

3155

// Compute Jacobian of affine map from reference cell

3156

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

3157

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

3158

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

3159

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

3160

3161

// Compute determinant of Jacobian

3162

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

3163

3164

// Compute inverse of Jacobian

3165

3166

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

3167

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

3168

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

3169

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

3170

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

3171

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

3172

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

3173

3174

// Map coordinates to the reference square

3175

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

3176

x = -1.0;

3177

else

3178

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

3179

y = 2.0*y - 1.0;

3180

3181

// Compute number of derivatives

3182

unsigned int num_derivatives = 1;

3183

3184

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

3185

num_derivatives *= 2;

3186

3187

3188

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

3189

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

3190

3191

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

3192

{

3193

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

3194

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

3195

combinations[j][k] = 0;

3196

}

3197

3198

// Generate combinations of derivatives

3199

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

3200

{

3201

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

3202

{

3203

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

3204

{

3205

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

3206

combinations[row][col] = 0;

3207

else

3208

{

3209

combinations[row][col] += 1;

3210

break;

3211

}

3212

}

3213

}

3214

}

3215

3216

// Compute inverse of Jacobian

3217

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

3218

3219

// Declare transformation matrix

3220

// Declare pointer to two dimensional array and initialise

3221

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

3222

3223

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

3224

{

3225

transform[j] = new double [num_derivatives];

3226

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

3227

transform[j][k] = 1;

3228

}

3229

3230

// Construct transformation matrix

3231

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

3232

{

3233

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

3234

{

3235

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

3236

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

3237

}

3238

}

3239

3240

// Reset values

3241

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

3242

values[j] = 0;

3243

3244

// Map degree of freedom to element degree of freedom

3245

const unsigned int dof = i;

3246

3247

// Generate scalings

3248

const double scalings_y_0 = 1;

3249

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

3250

3251

// Compute psitilde_a

3252

const double psitilde_a_0 = 1;

3253

const double psitilde_a_1 = x;

3254

3255

// Compute psitilde_bs

3256

const double psitilde_bs_0_0 = 1;

3257

const double psitilde_bs_0_1 = 1.5*y + 0.5;

3258

const double psitilde_bs_1_0 = 1;

3259

3260

// Compute basisvalues

3261

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

3262

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

3263

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

3264

3265

// Table(s) of coefficients

3266

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

3267

{{0.471404521, -0.288675135, -0.166666667},

3268

{0.471404521, 0.288675135, -0.166666667},

3269

{0.471404521, 0, 0.333333333}};

3270

3271

// Interesting (new) part

3272

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

3273

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

3274

{{0, 0, 0},

3275

{4.89897949, 0, 0},

3276

{0, 0, 0}};

3277

3278

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

3279

{{0, 0, 0},

3280

{2.44948974, 0, 0},

3281

{4.24264069, 0, 0}};

3282

3283

// Compute reference derivatives

3284

// Declare pointer to array of derivatives on FIAT element

3285

double *derivatives = new double [num_derivatives];

3286

3287

// Declare coefficients

3288

double coeff0_0 = 0;

3289

double coeff0_1 = 0;

3290

double coeff0_2 = 0;

3291

3292

// Declare new coefficients

3293

double new_coeff0_0 = 0;

3294

double new_coeff0_1 = 0;

3295

double new_coeff0_2 = 0;

3296

3297

// Loop possible derivatives

3298

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

3299

{

3300

// Get values from coefficients array

3301

new_coeff0_0 = coefficients0[dof][0];

3302

new_coeff0_1 = coefficients0[dof][1];

3303

new_coeff0_2 = coefficients0[dof][2];

3304

3305

// Loop derivative order

3306

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

3307

{

3308

// Update old coefficients

3309

coeff0_0 = new_coeff0_0;

3310

coeff0_1 = new_coeff0_1;

3311

coeff0_2 = new_coeff0_2;

3312

3313

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

3314

{

3315

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];

3316

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];

3317

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];

3318

}

3319

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

3320

{

3321

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];

3322

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];

3323

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];

3324

}

3325

3326

}

3327

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

3328

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;

3329

}

3330

3331

// Transform derivatives back to physical element

3332

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

3333

{

3334

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

3335

{

3336

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

3337

}

3338

}

3339

// Delete pointer to array of derivatives on FIAT element

3340

delete [] derivatives;

3341

3342

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

3343

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

3344

{

3345

delete [] combinations[row];

3346

delete [] transform[row];

3347

}

3348

3349

delete [] combinations;

3350

delete [] transform;

3351

}

3352

3353

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

3354

virtual void evaluate_basis_derivatives_all(unsigned int n,

3355

double* values,

3356

const double* coordinates,

3357

const ufc::cell& c) const

3358

{

3359

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

3360

}

3361

3362

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

3363

virtual double evaluate_dof(unsigned int i,

3364

const ufc::function& f,

3365

const ufc::cell& c) const

3366

{

3367

// The reference points, direction and weights:

3368

static const double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};

3369

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

3370

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

3371

3372

const double * const * x = c.coordinates;

3373

double result = 0.0;

3374

// Iterate over the points:

3375

// Evaluate basis functions for affine mapping

3376

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

3377

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

3378

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

3379

3380

// Compute affine mapping y = F(X)

3381

double y[2];

3382

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

3383

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

3384

3385

// Evaluate function at physical points

3386

double values[1];

3387

f.evaluate(values, y, c);

3388

3389

// Map function values using appropriate mapping

3390

// Affine map: Do nothing

3391

3392

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

3393

3394

// Take directional components

3395

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

3396

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

3397

// Multiply by weights

3398

result *= W[i][0];

3399

3400

return result;

3401

}

3402

3403

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

3404

virtual void evaluate_dofs(double* values,

3405

const ufc::function& f,

3406

const ufc::cell& c) const

3407

{

3408

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

3409

}

3410

3411

/// Interpolate vertex values from dof values

3412

virtual void interpolate_vertex_values(double* vertex_values,

3413

const double* dof_values,

3414

const ufc::cell& c) const

3415

{

3416

// Evaluate at vertices and use affine mapping

3417

vertex_values[0] = dof_values[0];

3418

vertex_values[1] = dof_values[1];

3419

vertex_values[2] = dof_values[2];

3420

}

3421

3422

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

3423

virtual unsigned int num_sub_elements() const

3424

{

3425

return 1;

3426

}

3427

3428

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

3429

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

3430

{

3431

return new poissondg_1_finite_element_1();

3432

}

3433

3434

};

3435

3436

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

3437

3438

class poissondg_1_finite_element_2: public ufc::finite_element

3439

{

3440

public:

3441

3442

/// Constructor

3443

poissondg_1_finite_element_2() : ufc::finite_element()

3444

{

3445

// Do nothing

3446

}

3447

3448

/// Destructor

3449

virtual ~poissondg_1_finite_element_2()

3450

{

3451

// Do nothing

3452

}

3453

3454

/// Return a string identifying the finite element

3455

virtual const char* signature() const

3456

{

3457

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

3458

}

3459

3460

/// Return the cell shape

3461

virtual ufc::shape cell_shape() const

3462

{

3463

return ufc::triangle;

3464

}

3465

3466

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

3467

virtual unsigned int space_dimension() const

3468

{

3469

return 3;

3470

}

3471

3472

/// Return the rank of the value space

3473

virtual unsigned int value_rank() const

3474

{

3475

return 0;

3476

}

3477

3478

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

3479

virtual unsigned int value_dimension(unsigned int i) const

3480

{

3481

return 1;

3482

}

3483

3484

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

3485

virtual void evaluate_basis(unsigned int i,

3486

double* values,

3487

const double* coordinates,

3488

const ufc::cell& c) const

3489

{

3490

// Extract vertex coordinates

3491

const double * const * element_coordinates = c.coordinates;

3492

3493

// Compute Jacobian of affine map from reference cell

3494

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

3495

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

3496

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

3497

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

3498

3499

// Compute determinant of Jacobian

3500

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

3501

3502

// Compute inverse of Jacobian

3503

3504

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

3505

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

3506

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

3507

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

3508

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

3509

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

3510

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

3511

3512

// Map coordinates to the reference square

3513

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

3514

x = -1.0;

3515

else

3516

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

3517

y = 2.0*y - 1.0;

3518

3519

// Reset values

3520

*values = 0;

3521

3522

// Map degree of freedom to element degree of freedom

3523

const unsigned int dof = i;

3524

3525

// Generate scalings

3526

const double scalings_y_0 = 1;

3527

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

3528

3529

// Compute psitilde_a

3530

const double psitilde_a_0 = 1;

3531

const double psitilde_a_1 = x;

3532

3533

// Compute psitilde_bs

3534

const double psitilde_bs_0_0 = 1;

3535

const double psitilde_bs_0_1 = 1.5*y + 0.5;

3536

const double psitilde_bs_1_0 = 1;

3537

3538

// Compute basisvalues

3539

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

3540

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

3541

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

3542

3543

// Table(s) of coefficients

3544

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

3545

{{0.471404521, -0.288675135, -0.166666667},

3546

{0.471404521, 0.288675135, -0.166666667},

3547

{0.471404521, 0, 0.333333333}};

3548

3549

// Extract relevant coefficients

3550

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

3551

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

3552

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

3553

3554

// Compute value(s)

3555

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2;

3556

}

3557

3558

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

3559

virtual void evaluate_basis_all(double* values,

3560

const double* coordinates,

3561

const ufc::cell& c) const

3562

{

3563

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

3564

}

3565

3566

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

3567

virtual void evaluate_basis_derivatives(unsigned int i,

3568

unsigned int n,

3569

double* values,

3570

const double* coordinates,

3571

const ufc::cell& c) const

3572

{

3573

// Extract vertex coordinates

3574

const double * const * element_coordinates = c.coordinates;

3575

3576

// Compute Jacobian of affine map from reference cell

3577

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

3578

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

3579

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

3580

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

3581

3582

// Compute determinant of Jacobian

3583

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

3584

3585

// Compute inverse of Jacobian

3586

3587

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

3588

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

3589

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

3590

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

3591

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

3592

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

3593

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

3594

3595

// Map coordinates to the reference square

3596

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

3597

x = -1.0;

3598

else

3599

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

3600

y = 2.0*y - 1.0;

3601

3602

// Compute number of derivatives

3603

unsigned int num_derivatives = 1;

3604

3605

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

3606

num_derivatives *= 2;

3607

3608

3609

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

3610

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

3611

3612

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

3613

{

3614

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

3615

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

3616

combinations[j][k] = 0;

3617

}

3618

3619

// Generate combinations of derivatives

3620

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

3621

{

3622

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

3623

{

3624

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

3625

{

3626

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

3627

combinations[row][col] = 0;

3628

else

3629

{

3630

combinations[row][col] += 1;

3631

break;

3632

}

3633

}

3634

}

3635

}

3636

3637

// Compute inverse of Jacobian

3638

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

3639

3640

// Declare transformation matrix

3641

// Declare pointer to two dimensional array and initialise

3642

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

3643

3644

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

3645

{

3646

transform[j] = new double [num_derivatives];

3647

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

3648

transform[j][k] = 1;

3649

}

3650

3651

// Construct transformation matrix

3652

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

3653

{

3654

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

3655

{

3656

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

3657

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

3658

}

3659

}

3660

3661

// Reset values

3662

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

3663

values[j] = 0;

3664

3665

// Map degree of freedom to element degree of freedom

3666

const unsigned int dof = i;

3667

3668

// Generate scalings

3669

const double scalings_y_0 = 1;

3670

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

3671

3672

// Compute psitilde_a

3673

const double psitilde_a_0 = 1;

3674

const double psitilde_a_1 = x;

3675

3676

// Compute psitilde_bs

3677

const double psitilde_bs_0_0 = 1;

3678

const double psitilde_bs_0_1 = 1.5*y + 0.5;

3679

const double psitilde_bs_1_0 = 1;

3680

3681

// Compute basisvalues

3682

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

3683

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

3684

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

3685

3686

// Table(s) of coefficients

3687

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

3688

{{0.471404521, -0.288675135, -0.166666667},

3689

{0.471404521, 0.288675135, -0.166666667},

3690

{0.471404521, 0, 0.333333333}};

3691

3692

// Interesting (new) part

3693

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

3694

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

3695

{{0, 0, 0},

3696

{4.89897949, 0, 0},

3697

{0, 0, 0}};

3698

3699

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

3700

{{0, 0, 0},

3701

{2.44948974, 0, 0},

3702

{4.24264069, 0, 0}};

3703

3704

// Compute reference derivatives

3705

// Declare pointer to array of derivatives on FIAT element

3706

double *derivatives = new double [num_derivatives];

3707

3708

// Declare coefficients

3709

double coeff0_0 = 0;

3710

double coeff0_1 = 0;

3711

double coeff0_2 = 0;

3712

3713

// Declare new coefficients

3714

double new_coeff0_0 = 0;

3715

double new_coeff0_1 = 0;

3716

double new_coeff0_2 = 0;

3717

3718

// Loop possible derivatives

3719

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

3720

{

3721

// Get values from coefficients array

3722

new_coeff0_0 = coefficients0[dof][0];

3723

new_coeff0_1 = coefficients0[dof][1];

3724

new_coeff0_2 = coefficients0[dof][2];

3725

3726

// Loop derivative order

3727

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

3728

{

3729

// Update old coefficients

3730

coeff0_0 = new_coeff0_0;

3731

coeff0_1 = new_coeff0_1;

3732

coeff0_2 = new_coeff0_2;

3733

3734

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

3735

{

3736

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0];

3737

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1];

3738

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2];

3739

}

3740

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

3741

{

3742

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0];

3743

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1];

3744

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2];

3745

}

3746

3747

}

3748

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

3749

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2;

3750

}

3751

3752

// Transform derivatives back to physical element

3753

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

3754

{

3755

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

3756

{

3757

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

3758

}

3759

}

3760

// Delete pointer to array of derivatives on FIAT element

3761

delete [] derivatives;

3762

3763

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

3764

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

3765

{

3766

delete [] combinations[row];

3767

delete [] transform[row];

3768

}

3769

3770

delete [] combinations;

3771

delete [] transform;

3772

}

3773

3774

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

3775

virtual void evaluate_basis_derivatives_all(unsigned int n,

3776

double* values,

3777

const double* coordinates,

3778

const ufc::cell& c) const

3779

{

3780

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

3781

}

3782

3783

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

3784

virtual double evaluate_dof(unsigned int i,

3785

const ufc::function& f,

3786

const ufc::cell& c) const

3787

{

3788

// The reference points, direction and weights:

3789

static const double X[3][1][2] = {{{0, 0}}, {{1, 0}}, {{0, 1}}};

3790

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

3791

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

3792

3793

const double * const * x = c.coordinates;

3794

double result = 0.0;

3795

// Iterate over the points:

3796

// Evaluate basis functions for affine mapping

3797

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

3798

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

3799

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

3800

3801

// Compute affine mapping y = F(X)

3802

double y[2];

3803

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

3804

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

3805

3806

// Evaluate function at physical points

3807

double values[1];

3808

f.evaluate(values, y, c);

3809

3810

// Map function values using appropriate mapping

3811

// Affine map: Do nothing

3812

3813

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

3814

3815

// Take directional components

3816

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

3817

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

3818

// Multiply by weights

3819

result *= W[i][0];

3820

3821

return result;

3822

}

3823

3824

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

3825

virtual void evaluate_dofs(double* values,

3826

const ufc::function& f,

3827

const ufc::cell& c) const

3828

{

3829

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

3830

}

3831

3832

/// Interpolate vertex values from dof values

3833

virtual void interpolate_vertex_values(double* vertex_values,

3834

const double* dof_values,

3835

const ufc::cell& c) const

3836

{

3837

// Evaluate at vertices and use affine mapping

3838

vertex_values[0] = dof_values[0];

3839

vertex_values[1] = dof_values[1];

3840

vertex_values[2] = dof_values[2];

3841

}

3842

3843

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

3844

virtual unsigned int num_sub_elements() const

3845

{

3846

return 1;

3847

}

3848

3849

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

3850

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

3851

{

3852

return new poissondg_1_finite_element_2();

3853

}

3854

3855

};

3856

3857

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

3858

/// degrees of freedom (dofs).

3859

3860

class poissondg_1_dof_map_0: public ufc::dof_map

3861

{

3862

private:

3863

3864

unsigned int __global_dimension;

3865

3866

public:

3867

3868

/// Constructor

3869

poissondg_1_dof_map_0() : ufc::dof_map()

3870

{

3871

__global_dimension = 0;

3872

}

3873

3874

/// Destructor

3875

virtual ~poissondg_1_dof_map_0()

3876

{

3877

// Do nothing

3878

}

3879

3880

/// Return a string identifying the dof map

3881

virtual const char* signature() const

3882

{

3883

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

3884

}

3885

3886

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

3887

virtual bool needs_mesh_entities(unsigned int d) const

3888

{

3889

switch ( d )

3890

{

3891

case 0:

3892

return false;

3893

break;

3894

case 1:

3895

return false;

3896

break;

3897

case 2:

3898

return true;

3899

break;

3900

}

3901

return false;

3902

}

3903

3904

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

3905

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

3906

{

3907

__global_dimension = 3*m.num_entities[2];

3908

return false;

3909

}

3910

3911

/// Initialize dof map for given cell

3912

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

3913

const ufc::cell& c)

3914

{

3915

// Do nothing

3916

}

3917

3918

/// Finish initialization of dof map for cells

3919

virtual void init_cell_finalize()

3920

{

3921

// Do nothing

3922

}

3923

3924

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

3925

virtual unsigned int global_dimension() const

3926

{

3927

return __global_dimension;

3928

}

3929

3930

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

3931

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

3932

{

3933

return 3;

3934

}

3935

3936

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

3937

virtual unsigned int max_local_dimension() const

3938

{

3939

return 3;

3940

}

3941

3942

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

3943

virtual unsigned int geometric_dimension() const

3944

{

3945

return 2;

3946

}

3947

3948

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

3949

virtual unsigned int num_facet_dofs() const

3950

{

3951

return 0;

3952

}

3953

3954

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

3955

virtual unsigned int num_entity_dofs(unsigned int d) const

3956

{

3957

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

3958

}

3959

3960

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

3961

virtual void tabulate_dofs(unsigned int* dofs,

3962

const ufc::mesh& m,

3963

const ufc::cell& c) const

3964

{

3965

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

3966

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

3967

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

3968

}

3969

3970

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

3971

virtual void tabulate_facet_dofs(unsigned int* dofs,

3972

unsigned int facet) const

3973

{

3974

switch ( facet )

3975

{

3976

case 0:

3977

3978

break;

3979

case 1:

3980

3981

break;

3982

case 2:

3983

3984

break;

3985

}

3986

}

3987

3988

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

3989

virtual void tabulate_entity_dofs(unsigned int* dofs,

3990

unsigned int d, unsigned int i) const

3991

{

3992

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

3993

}

3994

3995

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

3996

virtual void tabulate_coordinates(double** coordinates,

3997

const ufc::cell& c) const

3998

{

3999

const double * const * x = c.coordinates;

4000

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

4001

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

4002

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

4003

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

4004

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

4005

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

4006

}

4007

4008

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

4009

virtual unsigned int num_sub_dof_maps() const

4010

{

4011

return 1;

4012

}

4013

4014

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

4015

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

4016

{

4017

return new poissondg_1_dof_map_0();

4018

}

4019

4020

};

4021

4022

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

4023

/// degrees of freedom (dofs).

4024

4025

class poissondg_1_dof_map_1: public ufc::dof_map

4026

{

4027

private:

4028

4029

unsigned int __global_dimension;

4030

4031

public:

4032

4033

/// Constructor

4034

poissondg_1_dof_map_1() : ufc::dof_map()

4035

{

4036

__global_dimension = 0;

4037

}

4038

4039

/// Destructor

4040

virtual ~poissondg_1_dof_map_1()

4041

{

4042

// Do nothing

4043

}

4044

4045

/// Return a string identifying the dof map

4046

virtual const char* signature() const

4047

{

4048

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

4049

}

4050

4051

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

4052

virtual bool needs_mesh_entities(unsigned int d) const

4053

{

4054

switch ( d )

4055

{

4056

case 0:

4057

return false;

4058

break;

4059

case 1:

4060

return false;

4061

break;

4062

case 2:

4063

return true;

4064

break;

4065

}

4066

return false;

4067

}

4068

4069

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

4070

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

4071

{

4072

__global_dimension = 3*m.num_entities[2];

4073

return false;

4074

}

4075

4076

/// Initialize dof map for given cell

4077

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

4078

const ufc::cell& c)

4079

{

4080

// Do nothing

4081

}

4082

4083

/// Finish initialization of dof map for cells

4084

virtual void init_cell_finalize()

4085

{

4086

// Do nothing

4087

}

4088

4089

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

4090

virtual unsigned int global_dimension() const

4091

{

4092

return __global_dimension;

4093

}

4094

4095

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

4096

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

4097

{

4098

return 3;

4099

}

4100

4101

/// Return the maximum dimension of the local finite element function space

4102

virtual unsigned int max_local_dimension() const

4103

{

4104

return 3;

4105

}

4106

4107

// Return the geometric dimension of the coordinates this dof map provides

4108

virtual unsigned int geometric_dimension() const

4109

{

4110

return 2;

4111

}

4112

4113

/// Return the number of dofs on each cell facet

4114

virtual unsigned int num_facet_dofs() const

4115

{

4116

return 0;

4117

}

4118

4119

/// Return the number of dofs associated with each cell entity of dimension d

4120

virtual unsigned int num_entity_dofs(unsigned int d) const

4121

{

4122

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

4123

}

4124

4125

/// Tabulate the local-to-global mapping of dofs on a cell

4126

virtual void tabulate_dofs(unsigned int* dofs,

4127

const ufc::mesh& m,

4128

const ufc::cell& c) const

4129

{

4130

dofs[0] = 3*c.entity_indices[2][0];

4131

dofs[1] = 3*c.entity_indices[2][0] + 1;

4132

dofs[2] = 3*c.entity_indices[2][0] + 2;

4133

}

4134

4135

/// Tabulate the local-to-local mapping from facet dofs to cell dofs

4136

virtual void tabulate_facet_dofs(unsigned int* dofs,

4137

unsigned int facet) const

4138

{

4139

switch ( facet )

4140

{

4141

case 0:

4142

4143

break;

4144

case 1:

4145

4146

break;

4147

case 2:

4148

4149

break;

4150

}

4151

}

4152

4153

/// Tabulate the local-to-local mapping of dofs on entity (d, i)

4154

virtual void tabulate_entity_dofs(unsigned int* dofs,

4155

unsigned int d, unsigned int i) const

4156

{

4157

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

4158

}

4159

4160

/// Tabulate the coordinates of all dofs on a cell

4161

virtual void tabulate_coordinates(double** coordinates,

4162

const ufc::cell& c) const

4163

{

4164

const double * const * x = c.coordinates;

4165

coordinates[0][0] = x[0][0];

4166

coordinates[0][1] = x[0][1];

4167

coordinates[1][0] = x[1][0];

4168

coordinates[1][1] = x[1][1];

4169

coordinates[2][0] = x[2][0];

4170

coordinates[2][1] = x[2][1];

4171

}

4172

4173

/// Return the number of sub dof maps (for a mixed element)

4174

virtual unsigned int num_sub_dof_maps() const

4175

{

4176

return 1;

4177

}

4178

4179

/// Create a new dof_map for sub dof map i (for a mixed element)

4180

virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const

4181

{

4182

return new poissondg_1_dof_map_1();

4183

}

4184

4185

};

4186

4187

/// This class defines the interface for a local-to-global mapping of

4188

/// degrees of freedom (dofs).

4189

4190

class poissondg_1_dof_map_2: public ufc::dof_map

4191

{

4192

private:

4193

4194

unsigned int __global_dimension;

4195

4196

public:

4197

4198

/// Constructor

4199

poissondg_1_dof_map_2() : ufc::dof_map()

4200

{

4201

__global_dimension = 0;

4202

}

4203

4204

/// Destructor

4205

virtual ~poissondg_1_dof_map_2()

4206

{

4207

// Do nothing

4208

}

4209

4210

/// Return a string identifying the dof map

4211

virtual const char* signature() const

4212

{

4213

return "FFC dof map for FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1)";

4214

}

4215

4216

/// Return true iff mesh entities of topological dimension d are needed

4217

virtual bool needs_mesh_entities(unsigned int d) const

4218

{

4219

switch ( d )

4220

{

4221

case 0:

4222

return false;

4223

break;

4224

case 1:

4225

return false;

4226

break;

4227

case 2:

4228

return true;

4229

break;

4230

}

4231

return false;

4232

}

4233

4234

/// Initialize dof map for mesh (return true iff init_cell() is needed)

4235

virtual bool init_mesh(const ufc::mesh& m)

4236

{

4237

__global_dimension = 3*m.num_entities[2];

4238

return false;

4239

}

4240

4241

/// Initialize dof map for given cell

4242

virtual void init_cell(const ufc::mesh& m,

4243

const ufc::cell& c)

4244

{

4245

// Do nothing

4246

}

4247

4248

/// Finish initialization of dof map for cells

4249

virtual void init_cell_finalize()

4250

{

4251

// Do nothing

4252

}

4253

4254

/// Return the dimension of the global finite element function space

4255

virtual unsigned int global_dimension() const

4256

{

4257

return __global_dimension;

4258

}

4259

4260

/// Return the dimension of the local finite element function space for a cell

4261

virtual unsigned int local_dimension(const ufc::cell& c) const

4262

{

4263

return 3;

4264

}

4265

4266

/// Return the maximum dimension of the local finite element function space

4267

virtual unsigned int max_local_dimension() const

4268

{

4269

return 3;

4270

}

4271

4272

// Return the geometric dimension of the coordinates this dof map provides

4273

virtual unsigned int geometric_dimension() const

4274

{

4275

return 2;

4276

}

4277

4278

/// Return the number of dofs on each cell facet

4279

virtual unsigned int num_facet_dofs() const

4280

{

4281

return 0;

4282

}

4283

4284

/// Return the number of dofs associated with each cell entity of dimension d

4285

virtual unsigned int num_entity_dofs(unsigned int d) const

4286

{

4287

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

4288

}

4289

4290

/// Tabulate the local-to-global mapping of dofs on a cell

4291

virtual void tabulate_dofs(unsigned int* dofs,

4292

const ufc::mesh& m,

4293

const ufc::cell& c) const

4294

{

4295

dofs[0] = 3*c.entity_indices[2][0];

4296

dofs[1] = 3*c.entity_indices[2][0] + 1;

4297

dofs[2] = 3*c.entity_indices[2][0] + 2;

4298

}

4299

4300

/// Tabulate the local-to-local mapping from facet dofs to cell dofs

4301

virtual void tabulate_facet_dofs(unsigned int* dofs,

4302

unsigned int facet) const

4303

{

4304

switch ( facet )

4305

{

4306

case 0:

4307

4308

break;

4309

case 1:

4310

4311

break;

4312

case 2:

4313

4314

break;

4315

}

4316

}

4317

4318

/// Tabulate the local-to-local mapping of dofs on entity (d, i)

4319

virtual void tabulate_entity_dofs(unsigned int* dofs,

4320

unsigned int d, unsigned int i) const

4321

{

4322

throw std::runtime_error("Not implemented (introduced in UFC v1.1).");

4323

}

4324

4325

/// Tabulate the coordinates of all dofs on a cell

4326

virtual void tabulate_coordinates(double** coordinates,

4327

const ufc::cell& c) const

4328

{

4329

const double * const * x = c.coordinates;

4330

coordinates[0][0] = x[0][0];

4331

coordinates[0][1] = x[0][1];

4332

coordinates[1][0] = x[1][0];

4333

coordinates[1][1] = x[1][1];

4334

coordinates[2][0] = x[2][0];

4335

coordinates[2][1] = x[2][1];

4336

}

4337

4338

/// Return the number of sub dof maps (for a mixed element)

4339

virtual unsigned int num_sub_dof_maps() const

4340

{

4341

return 1;

4342

}

4343

4344

/// Create a new dof_map for sub dof map i (for a mixed element)

4345

virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const

4346

{

4347

return new poissondg_1_dof_map_2();

4348

}

4349

4350

};

4351

4352

/// This class defines the interface for the tabulation of the cell

4353

/// tensor corresponding to the local contribution to a form from

4354

/// the integral over a cell.

4355

4356

class poissondg_1_cell_integral_0_quadrature: public ufc::cell_integral

4357

{

4358

public:

4359

4360

/// Constructor

4361

poissondg_1_cell_integral_0_quadrature() : ufc::cell_integral()

4362

{

4363

// Do nothing

4364

}

4365

4366

/// Destructor

4367

virtual ~poissondg_1_cell_integral_0_quadrature()

4368

{

4369

// Do nothing

4370

}

4371

4372

/// Tabulate the tensor for the contribution from a local cell

4373

virtual void tabulate_tensor(double* A,

4374

const double * const * w,

4375

const ufc::cell& c) const

4376

{

4377

// Extract vertex coordinates

4378

const double * const * x = c.coordinates;

4379

4380

// Compute Jacobian of affine map from reference cell

4381

const double J_00 = x[1][0] - x[0][0];

4382

const double J_01 = x[2][0] - x[0][0];

4383

const double J_10 = x[1][1] - x[0][1];

4384

const double J_11 = x[2][1] - x[0][1];

4385

4386

// Compute determinant of Jacobian

4387

double detJ = J_00*J_11 - J_01*J_10;

4388

4389

// Compute inverse of Jacobian

4390

4391

// Set scale factor

4392

const double det = std::abs(detJ);

4393

4394

4395

// Array of quadrature weights

4396

static const double W4[4] = {0.159020691, 0.0909793091, 0.159020691, 0.0909793091};

4397

// Quadrature points on the UFC reference element: (0.178558728, 0.155051026), (0.0750311102, 0.644948974), (0.666390246, 0.155051026), (0.280019915, 0.644948974)

4398

4399

// Value of basis functions at quadrature points.

4400

static const double FE0[4][3] = \

4401

{{0.666390246, 0.178558728, 0.155051026},

4402

{0.280019915, 0.0750311102, 0.644948974},

4403

{0.178558728, 0.666390246, 0.155051026},

4404

{0.0750311102, 0.280019915, 0.644948974}};

4405

4406

4407

// Compute element tensor using UFL quadrature representation

4408

// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)

4409

// Total number of operations to compute element tensor: 72

4410

4411

// Loop quadrature points for integral

4412

// Number of operations to compute element tensor for following IP loop = 72

4413

for (unsigned int ip = 0; ip < 4; ip++)

4414

{

4415

4416

// Function declarations

4417

double F0 = 0;

4418

4419

// Total number of operations to compute function values = 6

4420

for (unsigned int r = 0; r < 3; r++)

4421

{

4422

F0 += FE0[ip][r]*w[0][r];

4423

}// end loop over 'r'

4424

4425

// Number of operations for primary indices: 12

4426

for (unsigned int j = 0; j < 3; j++)

4427

{

4428

// Number of operations to compute entry: 4

4429

A[j] += FE0[ip][j]*F0*W4[ip]*det;

4430

}// end loop over 'j'

4431

}// end loop over 'ip'

4432

}

4433

4434

};

4435

4436

/// This class defines the interface for the tabulation of the cell

4437

/// tensor corresponding to the local contribution to a form from

4438

/// the integral over a cell.

4439

4440

class poissondg_1_cell_integral_0: public ufc::cell_integral

4441

{

4442

private:

4443

4444

poissondg_1_cell_integral_0_quadrature integral_0_quadrature;

4445

4446

public:

4447

4448

/// Constructor

4449

poissondg_1_cell_integral_0() : ufc::cell_integral()

4450

{

4451

// Do nothing

4452

}

4453

4454

/// Destructor

4455

virtual ~poissondg_1_cell_integral_0()

4456

{

4457

// Do nothing

4458

}

4459

4460

/// Tabulate the tensor for the contribution from a local cell

4461

virtual void tabulate_tensor(double* A,

4462

const double * const * w,

4463

const ufc::cell& c) const

4464

{

4465

// Reset values of the element tensor block

4466

for (unsigned int j = 0; j < 3; j++)

4467

A[j] = 0;

4468

4469

// Add all contributions to element tensor

4470

integral_0_quadrature.tabulate_tensor(A, w, c);

4471

}

4472

4473

};

4474

4475

/// This class defines the interface for the tabulation of the

4476

/// exterior facet tensor corresponding to the local contribution to

4477

/// a form from the integral over an exterior facet.

4478

4479

class poissondg_1_exterior_facet_integral_0_quadrature: public ufc::exterior_facet_integral

4480

{

4481

public:

4482

4483

/// Constructor

4484

poissondg_1_exterior_facet_integral_0_quadrature() : ufc::exterior_facet_integral()

4485

{

4486

// Do nothing

4487

}

4488

4489

/// Destructor

4490

virtual ~poissondg_1_exterior_facet_integral_0_quadrature()

4491

{

4492

// Do nothing

4493

}

4494

4495

/// Tabulate the tensor for the contribution from a local exterior facet

4496

virtual void tabulate_tensor(double* A,

4497

const double * const * w,

4498

const ufc::cell& c,

4499

unsigned int facet) const

4500

{

4501

// Extract vertex coordinates

4502

const double * const * x = c.coordinates;

4503

4504

// Compute Jacobian of affine map from reference cell

4505

4506

// Compute determinant of Jacobian

4507

4508

// Compute inverse of Jacobian

4509

4510

// Vertices on edges

4511

static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};

4512

4513

// Get vertices

4514

const unsigned int v0 = edge_vertices[facet][0];

4515

const unsigned int v1 = edge_vertices[facet][1];

4516

4517

// Compute scale factor (length of edge scaled by length of reference interval)

4518

const double dx0 = x[v1][0] - x[v0][0];

4519

const double dx1 = x[v1][1] - x[v0][1];

4520

const double det = std::sqrt(dx0*dx0 + dx1*dx1);

4521

4522

4523

// Compute facet normals from the facet scale factor constants

4524

4525

4526

// Array of quadrature weights

4527

static const double W2[2] = {0.5, 0.5};

4528

// Quadrature points on the UFC reference element: (0.211324865), (0.788675135)

4529

4530

// Value of basis functions at quadrature points.

4531

static const double FE0_f0[2][3] = \

4532

{{0, 0.788675135, 0.211324865},

4533

{0, 0.211324865, 0.788675135}};

4534

4535

static const double FE0_f1[2][3] = \

4536

{{0.788675135, 0, 0.211324865},

4537

{0.211324865, 0, 0.788675135}};

4538

4539

static const double FE0_f2[2][3] = \

4540

{{0.788675135, 0.211324865, 0},

4541

{0.211324865, 0.788675135, 0}};

4542

4543

4544

// Compute element tensor using UFL quadrature representation

4545

// Optimisations: ('simplify expressions', False), ('ignore zero tables', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False)

4546

switch ( facet )

4547

{

4548

case 0:

4549

{

4550

// Total number of operations to compute element tensor (from this point): 36

4551

4552

// Loop quadrature points for integral

4553

// Number of operations to compute element tensor for following IP loop = 36

4554

for (unsigned int ip = 0; ip < 2; ip++)

4555

{

4556

4557

// Function declarations

4558

double F0 = 0;

4559

4560

// Total number of operations to compute function values = 6

4561

for (unsigned int r = 0; r < 3; r++)

4562

{

4563

F0 += FE0_f0[ip][r]*w[1][r];

4564

}// end loop over 'r'

4565

4566

// Number of operations for primary indices: 12

4567

for (unsigned int j = 0; j < 3; j++)

4568

{

4569

// Number of operations to compute entry: 4

4570

A[j] += FE0_f0[ip][j]*F0*W2[ip]*det;

4571

}// end loop over 'j'

4572

}// end loop over 'ip'

4573

}

4574

break;

4575

case 1:

4576

{

4577

// Total number of operations to compute element tensor (from this point): 36

4578

4579

// Loop quadrature points for integral

4580

// Number of operations to compute element tensor for following IP loop = 36

4581

for (unsigned int ip = 0; ip < 2; ip++)

4582

{

4583

4584

// Function declarations

4585

double F0 = 0;

4586

4587

// Total number of operations to compute function values = 6

4588

for (unsigned int r = 0; r < 3; r++)

4589

{

4590

F0 += FE0_f1[ip][r]*w[1][r];

4591

}// end loop over 'r'

4592

4593

// Number of operations for primary indices: 12

4594

for (unsigned int j = 0; j < 3; j++)

4595

{

4596

// Number of operations to compute entry: 4

4597

A[j] += FE0_f1[ip][j]*F0*W2[ip]*det;

4598

}// end loop over 'j'

4599

}// end loop over 'ip'

4600

}

4601

break;

4602

case 2:

4603

{

4604

// Total number of operations to compute element tensor (from this point): 36

4605

4606

// Loop quadrature points for integral

4607

// Number of operations to compute element tensor for following IP loop = 36

4608

for (unsigned int ip = 0; ip < 2; ip++)

4609

{

4610

4611

// Function declarations

4612

double F0 = 0;

4613

4614

// Total number of operations to compute function values = 6

4615

for (unsigned int r = 0; r < 3; r++)

4616

{

4617

F0 += FE0_f2[ip][r]*w[1][r];

4618

}// end loop over 'r'

4619

4620

// Number of operations for primary indices: 12

4621

for (unsigned int j = 0; j < 3; j++)

4622

{

4623

// Number of operations to compute entry: 4

4624

A[j] += FE0_f2[ip][j]*F0*W2[ip]*det;

4625

}// end loop over 'j'

4626

}// end loop over 'ip'

4627

}

4628

break;

4629

}

4630

}

4631

4632

};

4633

4634

/// This class defines the interface for the tabulation of the

4635

/// exterior facet tensor corresponding to the local contribution to

4636

/// a form from the integral over an exterior facet.

4637

4638

class poissondg_1_exterior_facet_integral_0: public ufc::exterior_facet_integral

4639

{

4640

private:

4641

4642

poissondg_1_exterior_facet_integral_0_quadrature integral_0_quadrature;

4643

4644

public:

4645

4646

/// Constructor

4647

poissondg_1_exterior_facet_integral_0() : ufc::exterior_facet_integral()

4648

{

4649

// Do nothing

4650

}

4651

4652

/// Destructor

4653

virtual ~poissondg_1_exterior_facet_integral_0()

4654

{

4655

// Do nothing

4656

}

4657

4658

/// Tabulate the tensor for the contribution from a local exterior facet

4659

virtual void tabulate_tensor(double* A,

4660

const double * const * w,

4661

const ufc::cell& c,

4662

unsigned int facet) const

4663

{

4664

// Reset values of the element tensor block

4665

for (unsigned int j = 0; j < 3; j++)

4666

A[j] = 0;

4667

4668

// Add all contributions to element tensor

4669

integral_0_quadrature.tabulate_tensor(A, w, c, facet);

4670

}

4671

4672

};

4673

4674

/// This class defines the interface for the assembly of the global

4675

/// tensor corresponding to a form with r + n arguments, that is, a

4676

/// mapping

4677

///

4678

/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R

4679

///

4680

/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r

4681

/// global tensor A is defined by

4682

///

4683

/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),

4684

///

4685

/// where each argument Vj represents the application to the

4686

/// sequence of basis functions of Vj and w1, w2, ..., wn are given

4687

/// fixed functions (coefficients).

4688

4689

class poissondg_form_1: public ufc::form

4690

{

4691

public:

4692

4693

/// Constructor

4694

poissondg_form_1() : ufc::form()

4695

{

4696

// Do nothing

4697

}

4698

4699

/// Destructor

4700

virtual ~poissondg_form_1()

4701

{

4702

// Do nothing

4703

}

4704

4705

/// Return a string identifying the form

4706

virtual const char* signature() const

4707

{

4708

return "Form([Integral(Product(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0), Function(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0)), Measure('cell', 0, None)), Integral(Product(BasisFunction(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 0), Function(FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 1), 1)), Measure('exterior_facet', 0, None))])";

4709

}

4710

4711

/// Return the rank of the global tensor (r)

4712

virtual unsigned int rank() const

4713

{

4714

return 1;

4715

}

4716

4717

/// Return the number of coefficients (n)

4718

virtual unsigned int num_coefficients() const

4719

{

4720

return 2;

4721

}

4722

4723

/// Return the number of cell integrals

4724

virtual unsigned int num_cell_integrals() const

4725

{

4726

return 1;

4727

}

4728

4729

/// Return the number of exterior facet integrals

4730

virtual unsigned int num_exterior_facet_integrals() const

4731

{

4732

return 1;

4733

}

4734

4735

/// Return the number of interior facet integrals

4736

virtual unsigned int num_interior_facet_integrals() const

4737

{

4738

return 0;

4739

}

4740

4741

/// Create a new finite element for argument function i

4742

virtual ufc::finite_element* create_finite_element(unsigned int i) const

4743

{

4744

switch ( i )

4745

{

4746

case 0:

4747

return new poissondg_1_finite_element_0();

4748

break;

4749

case 1:

4750

return new poissondg_1_finite_element_1();

4751

break;

4752

case 2:

4753

return new poissondg_1_finite_element_2();

4754

break;

4755

}

4756

return 0;

4757

}

4758

4759

/// Create a new dof map for argument function i

4760

virtual ufc::dof_map* create_dof_map(unsigned int i) const

4761

{

4762

switch ( i )

4763

{

4764

case 0:

4765

return new poissondg_1_dof_map_0();

4766

break;

4767

case 1:

4768

return new poissondg_1_dof_map_1();

4769

break;

4770

case 2:

4771

return new poissondg_1_dof_map_2();

4772

break;

4773

}

4774

return 0;

4775

}

4776

4777

/// Create a new cell integral on sub domain i

4778

virtual ufc::cell_integral* create_cell_integral(unsigned int i) const

4779

{

4780

return new poissondg_1_cell_integral_0();

4781

}

4782

4783

/// Create a new exterior facet integral on sub domain i

4784

virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const

4785

{

4786

return new poissondg_1_exterior_facet_integral_0();

4787

}

4788

4789

/// Create a new interior facet integral on sub domain i

4790

virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const

4791

{

4792

return 0;

4793

}

4794

4795

};

4796

4797

#endif