~njansson/dolfin/hpc

void eval(real block[], const AffineMap& map0, const AffineMap& map1, real det, unsigned int facet0, unsigned int facet1, unsigned int alignment) const;

};

class BilinearForm::TestElement : public dolfin::FiniteElement

{

public:

TestElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)

{

tensordims = new unsigned int [1];

tensordims[0] = 3;

subelements = new FiniteElement* [2];

subelements[0] = new SubElement_0();

subelements[1] = new SubElement_1();

}

~TestElement()

{

if ( tensordims ) delete [] tensordims;

if ( subelements )

{

for (unsigned int i = 0; i < elementdim(); i++)

delete subelements[i];

delete [] subelements;

}

inline unsigned int spacedim() const

{

return 15;

}

inline unsigned int shapedim() const

{

return 2;

}

inline unsigned int tensordim(unsigned int i) const

{

dolfin_assert(i < 1);

return tensordims[i];

}

inline unsigned int elementdim() const

{

return 2;

}

inline unsigned int rank() const

{

return 1;

}

100

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

101

{

102

nodes[0] = cell.entities(0)[0];

103

nodes[1] = cell.entities(0)[1];

104

nodes[2] = cell.entities(0)[2];

105

int offset = mesh.topology().size(0);

106

nodes[3] = offset + cell.entities(1)[0];

107

nodes[4] = offset + cell.entities(1)[1];

108

nodes[5] = offset + cell.entities(1)[2];

109

offset = offset + mesh.topology().size(1);

110

nodes[6] = offset + cell.entities(0)[0];

111

nodes[7] = offset + cell.entities(0)[1];

112

nodes[8] = offset + cell.entities(0)[2];

113

offset = offset + mesh.topology().size(0);

114

nodes[9] = offset + cell.entities(1)[0];

115

nodes[10] = offset + cell.entities(1)[1];

116

nodes[11] = offset + cell.entities(1)[2];

117

offset = offset + mesh.topology().size(1);

118

nodes[12] = offset + cell.entities(0)[0];

119

nodes[13] = offset + cell.entities(0)[1];

120

nodes[14] = offset + cell.entities(0)[2];

121

}

122

123

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

124

{

125

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

126

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

127

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

128

points[3] = map(5.000000000000000e-01, 5.000000000000000e-01);

129

points[4] = map(0.000000000000000e+00, 5.000000000000000e-01);

130

points[5] = map(5.000000000000000e-01, 0.000000000000000e+00);

131

points[6] = map(0.000000000000000e+00, 0.000000000000000e+00);

132

points[7] = map(1.000000000000000e+00, 0.000000000000000e+00);

133

points[8] = map(0.000000000000000e+00, 1.000000000000000e+00);

134

points[9] = map(5.000000000000000e-01, 5.000000000000000e-01);

135

points[10] = map(0.000000000000000e+00, 5.000000000000000e-01);

136

points[11] = map(5.000000000000000e-01, 0.000000000000000e+00);

137

components[0] = 0;

138

components[1] = 0;

139

components[2] = 0;

140

components[3] = 0;

141

components[4] = 0;

142

components[5] = 0;

143

components[6] = 1;

144

components[7] = 1;

145

components[8] = 1;

146

components[9] = 1;

147

components[10] = 1;

148

components[11] = 1;

149

points[12] = map(0.000000000000000e+00, 0.000000000000000e+00);

150

points[13] = map(1.000000000000000e+00, 0.000000000000000e+00);

151

points[14] = map(0.000000000000000e+00, 1.000000000000000e+00);

152

components[12] = 2;

153

components[13] = 2;

154

components[14] = 2;

155

}

156

157

void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const

158

{

159

// FIXME: Temporary fix for Lagrange elements

160

vertex_nodes[0] = vertex;

161

int offset = mesh.topology().size(0) + mesh.topology().size(1);

162

vertex_nodes[1] = offset + vertex;

163

offset = offset + mesh.topology().size(0) + mesh.topology().size(1);

164

vertex_nodes[2] = offset + vertex;

165

}

166

167

const FiniteElement& operator[] (unsigned int i) const

168

{

169

return *subelements[i];

170

}

171

172

FiniteElement& operator[] (unsigned int i)

173

{

174

return *subelements[i];

175

}

176

177

FiniteElementSpec spec() const

178

{

179

FiniteElementSpec s("mixed");

180

return s;

181

}

182

183

private:

184

185

class SubElement_0 : public dolfin::FiniteElement

186

{

187

public:

188

189

SubElement_0() : dolfin::FiniteElement(), tensordims(0), subelements(0)

190

{

191

tensordims = new unsigned int [1];

192

tensordims[0] = 2;

193

194

// Element is simple, don't need to initialize subelements

195

}

196

197

~SubElement_0()

198

{

199

if ( tensordims ) delete [] tensordims;

200

if ( subelements )

201

{

202

for (unsigned int i = 0; i < elementdim(); i++)

203

delete subelements[i];

204

delete [] subelements;

205

}

206

}

207

208

inline unsigned int spacedim() const

209

{

210

return 12;

211

}

212

213

inline unsigned int shapedim() const

214

{

215

return 2;

216

}

217

218

inline unsigned int tensordim(unsigned int i) const

219

{

220

dolfin_assert(i < 1);

221

return tensordims[i];

222

}

223

224

inline unsigned int elementdim() const

225

{

226

return 1;

227

}

228

229

inline unsigned int rank() const

230

{

231

return 1;

232

}

233

234

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

235

{

236

nodes[0] = cell.entities(0)[0];

237

nodes[1] = cell.entities(0)[1];

238

nodes[2] = cell.entities(0)[2];

239

int offset = mesh.topology().size(0);

240

nodes[3] = offset + cell.entities(1)[0];

241

nodes[4] = offset + cell.entities(1)[1];

242

nodes[5] = offset + cell.entities(1)[2];

243

offset = offset + mesh.topology().size(1);

244

nodes[6] = offset + cell.entities(0)[0];

245

nodes[7] = offset + cell.entities(0)[1];

246

nodes[8] = offset + cell.entities(0)[2];

247

offset = offset + mesh.topology().size(0);

248

nodes[9] = offset + cell.entities(1)[0];

249

nodes[10] = offset + cell.entities(1)[1];

250

nodes[11] = offset + cell.entities(1)[2];

251

}

252

253

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

254

{

255

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

256

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

257

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

258

points[3] = map(5.000000000000000e-01, 5.000000000000000e-01);

259

points[4] = map(0.000000000000000e+00, 5.000000000000000e-01);

260

points[5] = map(5.000000000000000e-01, 0.000000000000000e+00);

261

points[6] = map(0.000000000000000e+00, 0.000000000000000e+00);

262

points[7] = map(1.000000000000000e+00, 0.000000000000000e+00);

263

points[8] = map(0.000000000000000e+00, 1.000000000000000e+00);

264

points[9] = map(5.000000000000000e-01, 5.000000000000000e-01);

265

points[10] = map(0.000000000000000e+00, 5.000000000000000e-01);

266

points[11] = map(5.000000000000000e-01, 0.000000000000000e+00);

267

components[0] = 0;

268

components[1] = 0;

269

components[2] = 0;

270

components[3] = 0;

271

components[4] = 0;

272

components[5] = 0;

273

components[6] = 1;

274

components[7] = 1;

275

components[8] = 1;

276

components[9] = 1;

277

components[10] = 1;

278

components[11] = 1;

279

}

280

281

void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const

282

{

283

// FIXME: Temporary fix for Lagrange elements

284

vertex_nodes[0] = vertex;

285

int offset = mesh.topology().size(0) + mesh.topology().size(1);

286

vertex_nodes[1] = offset + vertex;

287

}

288

289

const FiniteElement& operator[] (unsigned int i) const

290

{

291

return *this;

292

}

293

294

FiniteElement& operator[] (unsigned int i)

295

{

296

return *this;

297

}

298

299

FiniteElementSpec spec() const

300

{

301

FiniteElementSpec s("Vector Lagrange", "triangle", 2, 2);

302

return s;

303

}

304

305

private:

306

307

unsigned int* tensordims;

308

FiniteElement** subelements;

309

310

};

311

312

class SubElement_1 : public dolfin::FiniteElement

313

{

314

public:

315

316

SubElement_1() : dolfin::FiniteElement(), tensordims(0), subelements(0)

317

{

318

// Element is scalar, don't need to initialize tensordims

319

320

// Element is simple, don't need to initialize subelements

321

}

322

323

~SubElement_1()

324

{

325

if ( tensordims ) delete [] tensordims;

326

if ( subelements )

327

{

328

for (unsigned int i = 0; i < elementdim(); i++)

329

delete subelements[i];

330

delete [] subelements;

331

}

332

}

333

334

inline unsigned int spacedim() const

335

{

336

return 3;

337

}

338

339

inline unsigned int shapedim() const

340

{

341

return 2;

342

}

343

344

inline unsigned int tensordim(unsigned int i) const

345

{

346

dolfin_error("Element is scalar.");

347

return 0;

348

}

349

350

inline unsigned int elementdim() const

351

{

352

return 1;

353

}

354

355

inline unsigned int rank() const

356

{

357

return 0;

358

}

359

360

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

361

{

362

nodes[0] = cell.entities(0)[0];

363

nodes[1] = cell.entities(0)[1];

364

nodes[2] = cell.entities(0)[2];

365

}

366

367

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

368

{

369

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

370

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

371

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

372

components[0] = 0;

373

components[1] = 0;

374

components[2] = 0;

375

}

376

377

void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const

378

{

379

// FIXME: Temporary fix for Lagrange elements

380

vertex_nodes[0] = vertex;

381

}

382

383

const FiniteElement& operator[] (unsigned int i) const

384

{

385

return *this;

386

}

387

388

FiniteElement& operator[] (unsigned int i)

389

{

390

return *this;

391

}

392

393

FiniteElementSpec spec() const

394

{

395

FiniteElementSpec s("Lagrange", "triangle", 1);

396

return s;

397

}

398

399

private:

400

401

unsigned int* tensordims;

402

FiniteElement** subelements;

403

404

};

405

406

unsigned int* tensordims;

407

FiniteElement** subelements;

408

409

};

410

411

class BilinearForm::TrialElement : public dolfin::FiniteElement

412

{

413

public:

414

415

TrialElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)

416

{

417

tensordims = new unsigned int [1];

418

tensordims[0] = 3;

419

420

subelements = new FiniteElement* [2];

421

subelements[0] = new SubElement_0();

422

subelements[1] = new SubElement_1();

423

}

424

425

~TrialElement()

426

{

427

if ( tensordims ) delete [] tensordims;

428

if ( subelements )

429

{

430

for (unsigned int i = 0; i < elementdim(); i++)

431

delete subelements[i];

432

delete [] subelements;

433

}

434

}

435

436

inline unsigned int spacedim() const

437

{

438

return 15;

439

}

440

441

inline unsigned int shapedim() const

442

{

443

return 2;

444

}

445

446

inline unsigned int tensordim(unsigned int i) const

447

{

448

dolfin_assert(i < 1);

449

return tensordims[i];

450

}

451

452

inline unsigned int elementdim() const

453

{

454

return 2;

455

}

456

457

inline unsigned int rank() const

458

{

459

return 1;

460

}

461

462

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

463

{

464

nodes[0] = cell.entities(0)[0];

465

nodes[1] = cell.entities(0)[1];

466

nodes[2] = cell.entities(0)[2];

467

int offset = mesh.topology().size(0);

468

nodes[3] = offset + cell.entities(1)[0];

469

nodes[4] = offset + cell.entities(1)[1];

470

nodes[5] = offset + cell.entities(1)[2];

471

offset = offset + mesh.topology().size(1);

472

nodes[6] = offset + cell.entities(0)[0];

473

nodes[7] = offset + cell.entities(0)[1];

474

nodes[8] = offset + cell.entities(0)[2];

475

offset = offset + mesh.topology().size(0);

476

nodes[9] = offset + cell.entities(1)[0];

477

nodes[10] = offset + cell.entities(1)[1];

478

nodes[11] = offset + cell.entities(1)[2];

479

offset = offset + mesh.topology().size(1);

480

nodes[12] = offset + cell.entities(0)[0];

481

nodes[13] = offset + cell.entities(0)[1];

482

nodes[14] = offset + cell.entities(0)[2];

483

}

484

485

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

486

{

487

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

488

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

489

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

490

points[3] = map(5.000000000000000e-01, 5.000000000000000e-01);

491

points[4] = map(0.000000000000000e+00, 5.000000000000000e-01);

492

points[5] = map(5.000000000000000e-01, 0.000000000000000e+00);

493

points[6] = map(0.000000000000000e+00, 0.000000000000000e+00);

494

points[7] = map(1.000000000000000e+00, 0.000000000000000e+00);

495

points[8] = map(0.000000000000000e+00, 1.000000000000000e+00);

496

points[9] = map(5.000000000000000e-01, 5.000000000000000e-01);

497

points[10] = map(0.000000000000000e+00, 5.000000000000000e-01);

498

points[11] = map(5.000000000000000e-01, 0.000000000000000e+00);

499

components[0] = 0;

500

components[1] = 0;

501

components[2] = 0;

502

components[3] = 0;

503

components[4] = 0;

504

components[5] = 0;

505

components[6] = 1;

506

components[7] = 1;

507

components[8] = 1;

508

components[9] = 1;

509

components[10] = 1;

510

components[11] = 1;

511

points[12] = map(0.000000000000000e+00, 0.000000000000000e+00);

512

points[13] = map(1.000000000000000e+00, 0.000000000000000e+00);

513

points[14] = map(0.000000000000000e+00, 1.000000000000000e+00);

514

components[12] = 2;

515

components[13] = 2;

516

components[14] = 2;

517

}

518

519

void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const

520

{

521

// FIXME: Temporary fix for Lagrange elements

522

vertex_nodes[0] = vertex;

523

int offset = mesh.topology().size(0) + mesh.topology().size(1);

524

vertex_nodes[1] = offset + vertex;

525

offset = offset + mesh.topology().size(0) + mesh.topology().size(1);

526

vertex_nodes[2] = offset + vertex;

527

}

528

529

const FiniteElement& operator[] (unsigned int i) const

530

{

531

return *subelements[i];

532

}

533

534

FiniteElement& operator[] (unsigned int i)

535

{

536

return *subelements[i];

537

}

538

539

FiniteElementSpec spec() const

540

{

541

FiniteElementSpec s("mixed");

542

return s;

543

}

544

545

private:

546

547

class SubElement_0 : public dolfin::FiniteElement

548

{

549

public:

550

551

SubElement_0() : dolfin::FiniteElement(), tensordims(0), subelements(0)

552

{

553

tensordims = new unsigned int [1];

554

tensordims[0] = 2;

555

556

// Element is simple, don't need to initialize subelements

557

}

558

559

~SubElement_0()

560

{

561

if ( tensordims ) delete [] tensordims;

562

if ( subelements )

563

{

564

for (unsigned int i = 0; i < elementdim(); i++)

565

delete subelements[i];

566

delete [] subelements;

567

}

568

}

569

570

inline unsigned int spacedim() const

571

{

572

return 12;

573

}

574

575

inline unsigned int shapedim() const

576

{

577

return 2;

578

}

579

580

inline unsigned int tensordim(unsigned int i) const

581

{

582

dolfin_assert(i < 1);

583

return tensordims[i];

584

}

585

586

inline unsigned int elementdim() const

587

{

588

return 1;

589

}

590

591

inline unsigned int rank() const

592

{

593

return 1;

594

}

595

596

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

597

{

598

nodes[0] = cell.entities(0)[0];

599

nodes[1] = cell.entities(0)[1];

600

nodes[2] = cell.entities(0)[2];

601

int offset = mesh.topology().size(0);

602

nodes[3] = offset + cell.entities(1)[0];

603

nodes[4] = offset + cell.entities(1)[1];

604

nodes[5] = offset + cell.entities(1)[2];

605

offset = offset + mesh.topology().size(1);

606

nodes[6] = offset + cell.entities(0)[0];

607

nodes[7] = offset + cell.entities(0)[1];

608

nodes[8] = offset + cell.entities(0)[2];

609

offset = offset + mesh.topology().size(0);

610

nodes[9] = offset + cell.entities(1)[0];

611

nodes[10] = offset + cell.entities(1)[1];

612

nodes[11] = offset + cell.entities(1)[2];

613

}

614

615

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

616

{

617

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

618

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

619

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

620

points[3] = map(5.000000000000000e-01, 5.000000000000000e-01);

621

points[4] = map(0.000000000000000e+00, 5.000000000000000e-01);

622

points[5] = map(5.000000000000000e-01, 0.000000000000000e+00);

623

points[6] = map(0.000000000000000e+00, 0.000000000000000e+00);

624

points[7] = map(1.000000000000000e+00, 0.000000000000000e+00);

625

points[8] = map(0.000000000000000e+00, 1.000000000000000e+00);

626

points[9] = map(5.000000000000000e-01, 5.000000000000000e-01);

627

points[10] = map(0.000000000000000e+00, 5.000000000000000e-01);

628

points[11] = map(5.000000000000000e-01, 0.000000000000000e+00);

629

components[0] = 0;

630

components[1] = 0;

631

components[2] = 0;

632

components[3] = 0;

633

components[4] = 0;

634

components[5] = 0;

635

components[6] = 1;

636

components[7] = 1;

637

components[8] = 1;

638

components[9] = 1;

639

components[10] = 1;

640

components[11] = 1;

641

}

642

643

void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const

644

{

645

// FIXME: Temporary fix for Lagrange elements

646

vertex_nodes[0] = vertex;

647

int offset = mesh.topology().size(0) + mesh.topology().size(1);

648

vertex_nodes[1] = offset + vertex;

649

}

650

651

const FiniteElement& operator[] (unsigned int i) const

652

{

653

return *this;

654

}

655

656

FiniteElement& operator[] (unsigned int i)

657

{

658

return *this;

659

}

660

661

FiniteElementSpec spec() const

662

{

663

FiniteElementSpec s("Vector Lagrange", "triangle", 2, 2);

664

return s;

665

}

666

667

private:

668

669

unsigned int* tensordims;

670

FiniteElement** subelements;

671

672

};

673

674

class SubElement_1 : public dolfin::FiniteElement

675

{

676

public:

677

678

SubElement_1() : dolfin::FiniteElement(), tensordims(0), subelements(0)

679

{

680

// Element is scalar, don't need to initialize tensordims

681

682

// Element is simple, don't need to initialize subelements

683

}

684

685

~SubElement_1()

686

{

687

if ( tensordims ) delete [] tensordims;

688

if ( subelements )

689

{

690

for (unsigned int i = 0; i < elementdim(); i++)

691

delete subelements[i];

692

delete [] subelements;

693

}

694

}

695

696

inline unsigned int spacedim() const

697

{

698

return 3;

699

}

700

701

inline unsigned int shapedim() const

702

{

703

return 2;

704

}

705

706

inline unsigned int tensordim(unsigned int i) const

707

{

708

dolfin_error("Element is scalar.");

709

return 0;

710

}

711

712

inline unsigned int elementdim() const

713

{

714

return 1;

715

}

716

717

inline unsigned int rank() const

718

{

719

return 0;

720

}

721

722

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

723

{

724

nodes[0] = cell.entities(0)[0];

725

nodes[1] = cell.entities(0)[1];

726

nodes[2] = cell.entities(0)[2];

727

}

728

729

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

730

{

731

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

732

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

733

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

734

components[0] = 0;

735

components[1] = 0;

736

components[2] = 0;

737

}

738

739

void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const

740

{

741

// FIXME: Temporary fix for Lagrange elements

742

vertex_nodes[0] = vertex;

743

}

744

745

const FiniteElement& operator[] (unsigned int i) const

746

{

747

return *this;

748

}

749

750

FiniteElement& operator[] (unsigned int i)

751

{

752

return *this;

753

}

754

755

FiniteElementSpec spec() const

756

{

757

FiniteElementSpec s("Lagrange", "triangle", 1);

758

return s;

759

}

760

761

private:

762

763

unsigned int* tensordims;

764

FiniteElement** subelements;

765

766

};

767

768

unsigned int* tensordims;

769

FiniteElement** subelements;

770

771

};

772

773

BilinearForm::BilinearForm() : dolfin::BilinearForm(0)

774

{

775

// Create finite element for test space

776

_test = new TestElement();

777

778

// Create finite element for trial space

779

_trial = new TrialElement();

780

}

781

782

// Contribution from the interior

783

bool BilinearForm::interior_contribution() const { return true; }

784

785

void BilinearForm::eval(real block[], const AffineMap& map, real det) const

786

{

787

// Compute geometry tensors

788

const real G0_0_0 = det*map.g00*map.g00 + det*map.g01*map.g01;

789

const real G0_0_1 = det*map.g00*map.g10 + det*map.g01*map.g11;

790

const real G0_1_0 = det*map.g10*map.g00 + det*map.g11*map.g01;

791

const real G0_1_1 = det*map.g10*map.g10 + det*map.g11*map.g11;

792

const real G1_0_0 = det*map.g00*map.g00 + det*map.g01*map.g01;

793

const real G1_0_1 = det*map.g00*map.g10 + det*map.g01*map.g11;

794

const real G1_1_0 = det*map.g10*map.g00 + det*map.g11*map.g01;

795

const real G1_1_1 = det*map.g10*map.g10 + det*map.g11*map.g11;

796

const real G2_0 = det*map.g00;

797

const real G2_1 = det*map.g10;

798

const real G3_0 = det*map.g01;

799

const real G3_1 = det*map.g11;

800

const real G4_0 = det*map.g00;

801

const real G4_1 = det*map.g10;

802

const real G5_0 = det*map.g01;

803

const real G5_1 = det*map.g11;

804

805

// Compute element tensor

806

block[0] = 4.999999999999992e-01*G0_0_0 + 4.999999999999991e-01*G0_0_1 + 4.999999999999991e-01*G0_1_0 + 4.999999999999991e-01*G0_1_1;

807

block[1] = 1.666666666666663e-01*G0_0_0 + 1.666666666666664e-01*G0_1_0;

808

block[2] = 1.666666666666665e-01*G0_0_1 + 1.666666666666665e-01*G0_1_1;

809

block[3] = 0.000000000000000e+00;

810

block[4] = -6.666666666666657e-01*G0_0_1 - 6.666666666666656e-01*G0_1_1;

811

block[5] = -6.666666666666655e-01*G0_0_0 - 6.666666666666655e-01*G0_1_0;

812

block[6] = 0.000000000000000e+00;

813

block[7] = 0.000000000000000e+00;

814

block[8] = 0.000000000000000e+00;

815

block[9] = 0.000000000000000e+00;

816

block[10] = 0.000000000000000e+00;

817

block[11] = 0.000000000000000e+00;

818

block[12] = 1.666666666666664e-01*G2_0 + 1.666666666666663e-01*G2_1;

819

block[13] = 0.000000000000000e+00;

820

block[14] = 0.000000000000000e+00;

821

block[15] = 1.666666666666663e-01*G0_0_0 + 1.666666666666664e-01*G0_0_1;

822

block[16] = 4.999999999999992e-01*G0_0_0;

823

block[17] = -1.666666666666664e-01*G0_0_1;

824

block[18] = 6.666666666666654e-01*G0_0_1;

825

block[19] = 0.000000000000000e+00;

826

block[20] = -6.666666666666654e-01*G0_0_0 - 6.666666666666656e-01*G0_0_1;

827

block[21] = 0.000000000000000e+00;

828

block[22] = 0.000000000000000e+00;

829

block[23] = 0.000000000000000e+00;

830

block[24] = 0.000000000000000e+00;

831

block[25] = 0.000000000000000e+00;

832

block[26] = 0.000000000000000e+00;

833

block[27] = 0.000000000000000e+00;

834

block[28] = -1.666666666666664e-01*G2_0;

835

block[29] = 0.000000000000000e+00;

836

block[30] = 1.666666666666665e-01*G0_1_0 + 1.666666666666665e-01*G0_1_1;

837

block[31] = -1.666666666666664e-01*G0_1_0;

838

block[32] = 4.999999999999991e-01*G0_1_1;

839

block[33] = 6.666666666666653e-01*G0_1_0;

840

block[34] = -6.666666666666654e-01*G0_1_0 - 6.666666666666656e-01*G0_1_1;

841

block[35] = 0.000000000000000e+00;

842

block[36] = 0.000000000000000e+00;

843

block[37] = 0.000000000000000e+00;

844

block[38] = 0.000000000000000e+00;

845

block[39] = 0.000000000000000e+00;

846

block[40] = 0.000000000000000e+00;

847

block[41] = 0.000000000000000e+00;

848

block[42] = 0.000000000000000e+00;

849

block[43] = 0.000000000000000e+00;

850

block[44] = -1.666666666666664e-01*G2_1;

851

block[45] = 0.000000000000000e+00;

852

block[46] = 6.666666666666654e-01*G0_1_0;

853

block[47] = 6.666666666666653e-01*G0_0_1;

854

block[48] = 1.333333333333330e+00*G0_0_0 + 6.666666666666652e-01*G0_0_1 + 6.666666666666652e-01*G0_1_0 + 1.333333333333330e+00*G0_1_1;

855

block[49] = -1.333333333333331e+00*G0_0_0 - 6.666666666666656e-01*G0_0_1 - 6.666666666666653e-01*G0_1_0;

856

block[50] = -6.666666666666652e-01*G0_0_1 - 6.666666666666653e-01*G0_1_0 - 1.333333333333331e+00*G0_1_1;

857

block[51] = 0.000000000000000e+00;

858

block[52] = 0.000000000000000e+00;

859

block[53] = 0.000000000000000e+00;

860

block[54] = 0.000000000000000e+00;

861

block[55] = 0.000000000000000e+00;

862

block[56] = 0.000000000000000e+00;

863

block[57] = -1.666666666666663e-01*G2_0 - 1.666666666666664e-01*G2_1;

864

block[58] = -1.666666666666663e-01*G2_0 - 3.333333333333327e-01*G2_1;

865

block[59] = -3.333333333333326e-01*G2_0 - 1.666666666666663e-01*G2_1;

866

block[60] = -6.666666666666656e-01*G0_1_0 - 6.666666666666656e-01*G0_1_1;

867

block[61] = 0.000000000000000e+00;

868

block[62] = -6.666666666666654e-01*G0_0_1 - 6.666666666666659e-01*G0_1_1;

869

block[63] = -1.333333333333331e+00*G0_0_0 - 6.666666666666653e-01*G0_0_1 - 6.666666666666656e-01*G0_1_0;

870

block[64] = 1.333333333333331e+00*G0_0_0 + 6.666666666666656e-01*G0_0_1 + 6.666666666666657e-01*G0_1_0 + 1.333333333333331e+00*G0_1_1;

871

block[65] = 6.666666666666654e-01*G0_0_1 + 6.666666666666655e-01*G0_1_0;

872

block[66] = 0.000000000000000e+00;

873

block[67] = 0.000000000000000e+00;

874

block[68] = 0.000000000000000e+00;

875

block[69] = 0.000000000000000e+00;

876

block[70] = 0.000000000000000e+00;

877

block[71] = 0.000000000000000e+00;

878

block[72] = 1.666666666666663e-01*G2_0 - 1.666666666666664e-01*G2_1;

879

block[73] = 1.666666666666664e-01*G2_0;

880

block[74] = 3.333333333333327e-01*G2_0 + 1.666666666666664e-01*G2_1;

881

block[75] = -6.666666666666655e-01*G0_0_0 - 6.666666666666656e-01*G0_0_1;

882

block[76] = -6.666666666666654e-01*G0_0_0 - 6.666666666666656e-01*G0_1_0;

883

block[77] = 0.000000000000000e+00;

884

block[78] = -6.666666666666652e-01*G0_0_1 - 6.666666666666653e-01*G0_1_0 - 1.333333333333331e+00*G0_1_1;

885

block[79] = 6.666666666666657e-01*G0_0_1 + 6.666666666666654e-01*G0_1_0;

886

block[80] = 1.333333333333331e+00*G0_0_0 + 6.666666666666656e-01*G0_0_1 + 6.666666666666656e-01*G0_1_0 + 1.333333333333331e+00*G0_1_1;

887

block[81] = 0.000000000000000e+00;

888

block[82] = 0.000000000000000e+00;

889

block[83] = 0.000000000000000e+00;

890

block[84] = 0.000000000000000e+00;

891

block[85] = 0.000000000000000e+00;

892

block[86] = 0.000000000000000e+00;

893

block[87] = -1.666666666666664e-01*G2_0 + 1.666666666666664e-01*G2_1;

894

block[88] = 1.666666666666664e-01*G2_0 + 3.333333333333328e-01*G2_1;

895

block[89] = 1.666666666666664e-01*G2_1;

896

block[90] = 0.000000000000000e+00;

897

block[91] = 0.000000000000000e+00;

898

block[92] = 0.000000000000000e+00;

899

block[93] = 0.000000000000000e+00;

900

block[94] = 0.000000000000000e+00;

901

block[95] = 0.000000000000000e+00;

902

block[96] = 4.999999999999991e-01*G1_0_0 + 4.999999999999991e-01*G1_0_1 + 4.999999999999991e-01*G1_1_0 + 4.999999999999991e-01*G1_1_1;

903

block[97] = 1.666666666666664e-01*G1_0_0 + 1.666666666666664e-01*G1_1_0;

904

block[98] = 1.666666666666666e-01*G1_0_1 + 1.666666666666665e-01*G1_1_1;

905

block[99] = 0.000000000000000e+00;

906

block[100] = -6.666666666666656e-01*G1_0_1 - 6.666666666666656e-01*G1_1_1;

907

block[101] = -6.666666666666655e-01*G1_0_0 - 6.666666666666655e-01*G1_1_0;

908

block[102] = 1.666666666666663e-01*G3_0 + 1.666666666666663e-01*G3_1;

909

block[103] = 0.000000000000000e+00;

910

block[104] = 0.000000000000000e+00;

911

block[105] = 0.000000000000000e+00;

912

block[106] = 0.000000000000000e+00;

913

block[107] = 0.000000000000000e+00;

914

block[108] = 0.000000000000000e+00;

915

block[109] = 0.000000000000000e+00;

916

block[110] = 0.000000000000000e+00;

917

block[111] = 1.666666666666663e-01*G1_0_0 + 1.666666666666664e-01*G1_0_1;

918

block[112] = 4.999999999999992e-01*G1_0_0;

919

block[113] = -1.666666666666664e-01*G1_0_1;

920

block[114] = 6.666666666666654e-01*G1_0_1;

921

block[115] = 0.000000000000000e+00;

922

block[116] = -6.666666666666654e-01*G1_0_0 - 6.666666666666656e-01*G1_0_1;

923

block[117] = 0.000000000000000e+00;

924

block[118] = -1.666666666666664e-01*G3_0;

925

block[119] = 0.000000000000000e+00;

926

block[120] = 0.000000000000000e+00;

927

block[121] = 0.000000000000000e+00;

928

block[122] = 0.000000000000000e+00;

929

block[123] = 0.000000000000000e+00;

930

block[124] = 0.000000000000000e+00;

931

block[125] = 0.000000000000000e+00;

932

block[126] = 1.666666666666666e-01*G1_1_0 + 1.666666666666665e-01*G1_1_1;

933

block[127] = -1.666666666666664e-01*G1_1_0;

934

block[128] = 4.999999999999991e-01*G1_1_1;

935

block[129] = 6.666666666666653e-01*G1_1_0;

936

block[130] = -6.666666666666654e-01*G1_1_0 - 6.666666666666656e-01*G1_1_1;

937

block[131] = 0.000000000000000e+00;

938

block[132] = 0.000000000000000e+00;

939

block[133] = 0.000000000000000e+00;

940

block[134] = -1.666666666666664e-01*G3_1;

941

block[135] = 0.000000000000000e+00;

942

block[136] = 0.000000000000000e+00;

943

block[137] = 0.000000000000000e+00;

944

block[138] = 0.000000000000000e+00;

945

block[139] = 0.000000000000000e+00;

946

block[140] = 0.000000000000000e+00;

947

block[141] = 0.000000000000000e+00;

948

block[142] = 6.666666666666654e-01*G1_1_0;

949

block[143] = 6.666666666666653e-01*G1_0_1;

950

block[144] = 1.333333333333330e+00*G1_0_0 + 6.666666666666652e-01*G1_0_1 + 6.666666666666652e-01*G1_1_0 + 1.333333333333330e+00*G1_1_1;

951

block[145] = -1.333333333333331e+00*G1_0_0 - 6.666666666666656e-01*G1_0_1 - 6.666666666666653e-01*G1_1_0;

952

block[146] = -6.666666666666652e-01*G1_0_1 - 6.666666666666653e-01*G1_1_0 - 1.333333333333331e+00*G1_1_1;

953

block[147] = -1.666666666666663e-01*G3_0 - 1.666666666666664e-01*G3_1;

954

block[148] = -1.666666666666663e-01*G3_0 - 3.333333333333327e-01*G3_1;

955

block[149] = -3.333333333333326e-01*G3_0 - 1.666666666666663e-01*G3_1;

956

block[150] = 0.000000000000000e+00;

957

block[151] = 0.000000000000000e+00;

958

block[152] = 0.000000000000000e+00;

959

block[153] = 0.000000000000000e+00;

960

block[154] = 0.000000000000000e+00;

961

block[155] = 0.000000000000000e+00;

962

block[156] = -6.666666666666656e-01*G1_1_0 - 6.666666666666656e-01*G1_1_1;

963

block[157] = 0.000000000000000e+00;

964

block[158] = -6.666666666666654e-01*G1_0_1 - 6.666666666666659e-01*G1_1_1;

965

block[159] = -1.333333333333331e+00*G1_0_0 - 6.666666666666653e-01*G1_0_1 - 6.666666666666656e-01*G1_1_0;

966

block[160] = 1.333333333333331e+00*G1_0_0 + 6.666666666666656e-01*G1_0_1 + 6.666666666666657e-01*G1_1_0 + 1.333333333333331e+00*G1_1_1;

967

block[161] = 6.666666666666654e-01*G1_0_1 + 6.666666666666655e-01*G1_1_0;

968

block[162] = 1.666666666666663e-01*G3_0 - 1.666666666666664e-01*G3_1;

969

block[163] = 1.666666666666664e-01*G3_0;

970

block[164] = 3.333333333333327e-01*G3_0 + 1.666666666666664e-01*G3_1;

971

block[165] = 0.000000000000000e+00;

972

block[166] = 0.000000000000000e+00;

973

block[167] = 0.000000000000000e+00;

974

block[168] = 0.000000000000000e+00;

975

block[169] = 0.000000000000000e+00;

976

block[170] = 0.000000000000000e+00;

977

block[171] = -6.666666666666656e-01*G1_0_0 - 6.666666666666656e-01*G1_0_1;

978

block[172] = -6.666666666666654e-01*G1_0_0 - 6.666666666666656e-01*G1_1_0;

979

block[173] = 0.000000000000000e+00;

980

block[174] = -6.666666666666652e-01*G1_0_1 - 6.666666666666653e-01*G1_1_0 - 1.333333333333331e+00*G1_1_1;

981

block[175] = 6.666666666666657e-01*G1_0_1 + 6.666666666666654e-01*G1_1_0;

982

block[176] = 1.333333333333331e+00*G1_0_0 + 6.666666666666656e-01*G1_0_1 + 6.666666666666656e-01*G1_1_0 + 1.333333333333331e+00*G1_1_1;

983

block[177] = -1.666666666666664e-01*G3_0 + 1.666666666666664e-01*G3_1;

984

block[178] = 1.666666666666664e-01*G3_0 + 3.333333333333328e-01*G3_1;

985

block[179] = 1.666666666666664e-01*G3_1;

986

block[180] = -1.666666666666664e-01*G4_0 - 1.666666666666663e-01*G4_1;

987

block[181] = 0.000000000000000e+00;

988

block[182] = 0.000000000000000e+00;

989

block[183] = 1.666666666666663e-01*G4_0 + 1.666666666666664e-01*G4_1;

990

block[184] = -1.666666666666663e-01*G4_0 + 1.666666666666664e-01*G4_1;

991

block[185] = 1.666666666666664e-01*G4_0 - 1.666666666666664e-01*G4_1;

992

block[186] = -1.666666666666663e-01*G5_0 - 1.666666666666663e-01*G5_1;

993

block[187] = 0.000000000000000e+00;

994

block[188] = 0.000000000000000e+00;

995

block[189] = 1.666666666666663e-01*G5_0 + 1.666666666666664e-01*G5_1;

996

block[190] = -1.666666666666663e-01*G5_0 + 1.666666666666664e-01*G5_1;

997

block[191] = 1.666666666666664e-01*G5_0 - 1.666666666666664e-01*G5_1;

998

block[192] = 0.000000000000000e+00;

999

block[193] = 0.000000000000000e+00;

1000

block[194] = 0.000000000000000e+00;

1001

block[195] = 0.000000000000000e+00;

1002

block[196] = 1.666666666666664e-01*G4_0;

1003

block[197] = 0.000000000000000e+00;

1004

block[198] = 1.666666666666663e-01*G4_0 + 3.333333333333326e-01*G4_1;

1005

block[199] = -1.666666666666664e-01*G4_0;

1006

block[200] = -1.666666666666664e-01*G4_0 - 3.333333333333328e-01*G4_1;

1007

block[201] = 0.000000000000000e+00;

1008

block[202] = 1.666666666666664e-01*G5_0;

1009

block[203] = 0.000000000000000e+00;

1010

block[204] = 1.666666666666663e-01*G5_0 + 3.333333333333326e-01*G5_1;

1011

block[205] = -1.666666666666664e-01*G5_0;

1012

block[206] = -1.666666666666664e-01*G5_0 - 3.333333333333328e-01*G5_1;

1013

block[207] = 0.000000000000000e+00;

1014

block[208] = 0.000000000000000e+00;

1015

block[209] = 0.000000000000000e+00;

1016

block[210] = 0.000000000000000e+00;

1017

block[211] = 0.000000000000000e+00;

1018

block[212] = 1.666666666666663e-01*G4_1;

1019

block[213] = 3.333333333333326e-01*G4_0 + 1.666666666666663e-01*G4_1;

1020

block[214] = -3.333333333333327e-01*G4_0 - 1.666666666666664e-01*G4_1;

1021

block[215] = -1.666666666666664e-01*G4_1;

1022

block[216] = 0.000000000000000e+00;

1023

block[217] = 0.000000000000000e+00;

1024

block[218] = 1.666666666666663e-01*G5_1;

1025

block[219] = 3.333333333333326e-01*G5_0 + 1.666666666666663e-01*G5_1;

1026

block[220] = -3.333333333333327e-01*G5_0 - 1.666666666666664e-01*G5_1;

1027

block[221] = -1.666666666666664e-01*G5_1;

1028

block[222] = 0.000000000000000e+00;

1029

block[223] = 0.000000000000000e+00;

1030

block[224] = 0.000000000000000e+00;

1031

}

1032

1033

// No contribution from the boundary

1034

bool BilinearForm::boundary_contribution() const { return false; }

1035

1036

void BilinearForm::eval(real block[], const AffineMap& map, real det, unsigned int facet) const {}

1037

1038

// No contribution from interior boundaries

1039

bool BilinearForm::interior_boundary_contribution() const { return false; }

1040

1041

void BilinearForm::eval(real block[], const AffineMap& map0, const AffineMap& map1, real det, unsigned int facet0, unsigned int facet1, unsigned int alignment) const {}

1009

1042

1010

1043

/// This class contains the form to be evaluated, including

1011

1044

/// contributions from the interior and boundary of the domain.

1013

1046

class LinearForm : public dolfin::LinearForm

1014

1047

{

1015

1048

public:

1016

1017

class TestElement : public dolfin::FiniteElement

1018

{

1019

public:

1020

1021

TestElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)

1022

{

1023

tensordims = new unsigned int [1];

1024

tensordims[0] = 3;

1025

1026

subelements = new FiniteElement* [2];

1027

subelements[0] = new SubElement_0();

1028

subelements[1] = new SubElement_1();

1029

}

1030

1031

~TestElement()

1032

{

1033

if ( tensordims ) delete [] tensordims;

1034

if ( subelements )

1035

{

1036

for (unsigned int i = 0; i < elementdim(); i++)

1037

delete subelements[i];

1038

delete [] subelements;

1039

}

1040

}

1041

1042

inline unsigned int spacedim() const

1043

{

1044

return 15;

1045

}

1046

1047

inline unsigned int shapedim() const

1048

{

1049

return 2;

1050

}

1051

1052

inline unsigned int tensordim(unsigned int i) const

1053

{

1054

dolfin_assert(i < 1);

1055

return tensordims[i];

1056

}

1057

1058

inline unsigned int elementdim() const

1059

{

1060

return 2;

1061

}

1062

1063

inline unsigned int rank() const

1064

{

1065

return 1;

1066

}

1067

1068

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

1069

{

1070

nodes[0] = cell.vertexID(0);

1071

nodes[1] = cell.vertexID(1);

1072

nodes[2] = cell.vertexID(2);

1073

int offset = mesh.numVertices();

1074

nodes[3] = offset + cell.edgeID(0);

1075

nodes[4] = offset + cell.edgeID(1);

1076

nodes[5] = offset + cell.edgeID(2);

1077

offset = offset + mesh.numEdges();

1078

nodes[6] = offset + cell.vertexID(0);

1079

nodes[7] = offset + cell.vertexID(1);

1080

nodes[8] = offset + cell.vertexID(2);

1081

offset = offset + mesh.numVertices();

1082

nodes[9] = offset + cell.edgeID(0);

1083

nodes[10] = offset + cell.edgeID(1);

1084

nodes[11] = offset + cell.edgeID(2);

1085

offset = offset + mesh.numEdges();

1086

nodes[12] = offset + cell.vertexID(0);

1087

nodes[13] = offset + cell.vertexID(1);

1088

nodes[14] = offset + cell.vertexID(2);

1089

}

1090

1091

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

1092

{

1093

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

1094

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

1095

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

1096

points[3] = map(5.000000000000000e-01, 5.000000000000000e-01);

1097

points[4] = map(0.000000000000000e+00, 5.000000000000000e-01);

1098

points[5] = map(5.000000000000000e-01, 0.000000000000000e+00);

1099

points[6] = map(0.000000000000000e+00, 0.000000000000000e+00);

1100

points[7] = map(1.000000000000000e+00, 0.000000000000000e+00);

1101

points[8] = map(0.000000000000000e+00, 1.000000000000000e+00);

1102

points[9] = map(5.000000000000000e-01, 5.000000000000000e-01);

1103

points[10] = map(0.000000000000000e+00, 5.000000000000000e-01);

1104

points[11] = map(5.000000000000000e-01, 0.000000000000000e+00);

1105

components[0] = 0;

1106

components[1] = 0;

1107

components[2] = 0;

1108

components[3] = 0;

1109

components[4] = 0;

1110

components[5] = 0;

1111

components[6] = 1;

1112

components[7] = 1;

1113

components[8] = 1;

1114

components[9] = 1;

1115

components[10] = 1;

1116

components[11] = 1;

1117

points[12] = map(0.000000000000000e+00, 0.000000000000000e+00);

1118

points[13] = map(1.000000000000000e+00, 0.000000000000000e+00);

1119

points[14] = map(0.000000000000000e+00, 1.000000000000000e+00);

1120

components[12] = 2;

1121

components[13] = 2;

1122

components[14] = 2;

1123

}

1124

1125

void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const

1126

{

1127

// FIXME: Temporary fix for Lagrange elements

1128

values[0] = x[vertex];

1129

int offset = mesh.numVertices() + mesh.numEdges();

1130

values[1] = x[offset + vertex];

1131

offset = offset + mesh.numVertices() + mesh.numEdges();

1132

values[2] = x[offset + vertex];

1133

}

1134

1135

const FiniteElement& operator[] (unsigned int i) const

1136

{

1137

return *subelements[i];

1138

}

1139

1140

FiniteElement& operator[] (unsigned int i)

1141

{

1142

return *subelements[i];

1143

}

1144

1145

FiniteElementSpec spec() const

1146

{

1147

FiniteElementSpec s("mixed");

1148

return s;

1149

}

1150

1151

private:

1152

1153

class SubElement_0 : public dolfin::FiniteElement

1154

{

1155

public:

1156

1157

SubElement_0() : dolfin::FiniteElement(), tensordims(0), subelements(0)

1158

{

1159

tensordims = new unsigned int [1];

1160

tensordims[0] = 2;

1161

1162

// Element is simple, don't need to initialize subelements

1163

}

1164

1165

~SubElement_0()

1166

{

1167

if ( tensordims ) delete [] tensordims;

1168

if ( subelements )

1169

{

1170

for (unsigned int i = 0; i < elementdim(); i++)

1171

delete subelements[i];

1172

delete [] subelements;

1173

}

1174

}

1175

1176

inline unsigned int spacedim() const

1177

{

1178

return 12;

1179

}

1180

1181

inline unsigned int shapedim() const

1182

{

1183

return 2;

1184

}

1185

1186

inline unsigned int tensordim(unsigned int i) const

1187

{

1188

dolfin_assert(i < 1);

1189

return tensordims[i];

1190

}

1191

1192

inline unsigned int elementdim() const

1193

{

1194

return 1;

1195

}

1196

1197

inline unsigned int rank() const

1198

{

1199

return 1;

1200

}

1201

1202

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

1203

{

1204

nodes[0] = cell.vertexID(0);

1205

nodes[1] = cell.vertexID(1);

1206

nodes[2] = cell.vertexID(2);

1207

int offset = mesh.numVertices();

1208

nodes[3] = offset + cell.edgeID(0);

1209

nodes[4] = offset + cell.edgeID(1);

1210

nodes[5] = offset + cell.edgeID(2);

1211

offset = offset + mesh.numEdges();

1212

nodes[6] = offset + cell.vertexID(0);

1213

nodes[7] = offset + cell.vertexID(1);

1214

nodes[8] = offset + cell.vertexID(2);

1215

offset = offset + mesh.numVertices();

1216

nodes[9] = offset + cell.edgeID(0);

1217

nodes[10] = offset + cell.edgeID(1);

1218

nodes[11] = offset + cell.edgeID(2);

1219

}

1220

1221

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

1222

{

1223

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

1224

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

1225

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

1226

points[3] = map(5.000000000000000e-01, 5.000000000000000e-01);

1227

points[4] = map(0.000000000000000e+00, 5.000000000000000e-01);

1228

points[5] = map(5.000000000000000e-01, 0.000000000000000e+00);

1229

points[6] = map(0.000000000000000e+00, 0.000000000000000e+00);

1230

points[7] = map(1.000000000000000e+00, 0.000000000000000e+00);

1231

points[8] = map(0.000000000000000e+00, 1.000000000000000e+00);

1232

points[9] = map(5.000000000000000e-01, 5.000000000000000e-01);

1233

points[10] = map(0.000000000000000e+00, 5.000000000000000e-01);

1234

points[11] = map(5.000000000000000e-01, 0.000000000000000e+00);

1235

components[0] = 0;

1236

components[1] = 0;

1237

components[2] = 0;

1238

components[3] = 0;

1239

components[4] = 0;

1240

components[5] = 0;

1241

components[6] = 1;

1242

components[7] = 1;

1243

components[8] = 1;

1244

components[9] = 1;

1245

components[10] = 1;

1246

components[11] = 1;

1247

}

1248

1249

void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const

1250

{

1251

// FIXME: Temporary fix for Lagrange elements

1252

values[0] = x[vertex];

1253

int offset = mesh.numVertices() + mesh.numEdges();

1254

values[1] = x[offset + vertex];

1255

}

1256

1257

const FiniteElement& operator[] (unsigned int i) const

1258

{

1259

return *this;

1260

}

1261

1262

FiniteElement& operator[] (unsigned int i)

1263

{

1264

return *this;

1265

}

1266

1267

FiniteElementSpec spec() const

1268

{

1269

FiniteElementSpec s("Vector Lagrange", "triangle", 2, 2);

1270

return s;

1271

}

1272

1273

private:

1274

1275

unsigned int* tensordims;

1276

FiniteElement** subelements;

1277

1278

};

1279

1280

class SubElement_1 : public dolfin::FiniteElement

1281

{

1282

public:

1283

1284

SubElement_1() : dolfin::FiniteElement(), tensordims(0), subelements(0)

1285

{

1286

// Element is scalar, don't need to initialize tensordims

1287

1288

// Element is simple, don't need to initialize subelements

1289

}

1290

1291

~SubElement_1()

1292

{

1293

if ( tensordims ) delete [] tensordims;

1294

if ( subelements )

1295

{

1296

for (unsigned int i = 0; i < elementdim(); i++)

1297

delete subelements[i];

1298

delete [] subelements;

1299

}

1300

}

1301

1302

inline unsigned int spacedim() const

1303

{

1304

return 3;

1305

}

1306

1307

inline unsigned int shapedim() const

1308

{

1309

return 2;

1310

}

1311

1312

inline unsigned int tensordim(unsigned int i) const

1313

{

1314

dolfin_error("Element is scalar.");

1315

return 0;

1316

}

1317

1318

inline unsigned int elementdim() const

1319

{

1320

return 1;

1321

}

1322

1323

inline unsigned int rank() const

1324

{

1325

return 0;

1326

}

1327

1328

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

1329

{

1330

nodes[0] = cell.vertexID(0);

1331

nodes[1] = cell.vertexID(1);

1332

nodes[2] = cell.vertexID(2);

1333

}

1334

1335

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

1336

{

1337

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

1338

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

1339

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

1340

components[0] = 0;

1341

components[1] = 0;

1342

components[2] = 0;

1343

}

1344

1345

void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const

1346

{

1347

// FIXME: Temporary fix for Lagrange elements

1348

values[0] = x[vertex];

1349

}

1350

1351

const FiniteElement& operator[] (unsigned int i) const

1352

{

1353

return *this;

1354

}

1355

1356

FiniteElement& operator[] (unsigned int i)

1357

{

1358

return *this;

1359

}

1360

1361

FiniteElementSpec spec() const

1362

{

1363

FiniteElementSpec s("Lagrange", "triangle", 1);

1364

return s;

1365

}

1366

1367

private:

1368

1369

unsigned int* tensordims;

1370

FiniteElement** subelements;

1371

1372

};

1373

1374

unsigned int* tensordims;

1375

FiniteElement** subelements;

1376

1377

};

1378

1379

class FunctionElement_0 : public dolfin::FiniteElement

1380

{

1381

public:

1382

1383

FunctionElement_0() : dolfin::FiniteElement(), tensordims(0), subelements(0)

1384

{

1385

tensordims = new unsigned int [1];

1386

tensordims[0] = 2;

1387

1388

// Element is simple, don't need to initialize subelements

1389

}

1390

1391

~FunctionElement_0()

1392

{

1393

if ( tensordims ) delete [] tensordims;

1394

if ( subelements )

1395

{

1396

for (unsigned int i = 0; i < elementdim(); i++)

1397

delete subelements[i];

1398

delete [] subelements;

1399

}

1400

}

1401

1402

inline unsigned int spacedim() const

1403

{

1404

return 12;

1405

}

1406

1407

inline unsigned int shapedim() const

1408

{

1409

return 2;

1410

}

1411

1412

inline unsigned int tensordim(unsigned int i) const

1413

{

1414

dolfin_assert(i < 1);

1415

return tensordims[i];

1416

}

1417

1418

inline unsigned int elementdim() const

1419

{

1420

return 1;

1421

}

1422

1423

inline unsigned int rank() const

1424

{

1425

return 1;

1426

}

1427

1428

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

1429

{

1430

nodes[0] = cell.vertexID(0);

1431

nodes[1] = cell.vertexID(1);

1432

nodes[2] = cell.vertexID(2);

1433

int offset = mesh.numVertices();

1434

nodes[3] = offset + cell.edgeID(0);

1435

nodes[4] = offset + cell.edgeID(1);

1436

nodes[5] = offset + cell.edgeID(2);

1437

offset = offset + mesh.numEdges();

1438

nodes[6] = offset + cell.vertexID(0);

1439

nodes[7] = offset + cell.vertexID(1);

1440

nodes[8] = offset + cell.vertexID(2);

1441

offset = offset + mesh.numVertices();

1442

nodes[9] = offset + cell.edgeID(0);

1443

nodes[10] = offset + cell.edgeID(1);

1444

nodes[11] = offset + cell.edgeID(2);

1445

}

1446

1447

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

1448

{

1449

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

1450

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

1451

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

1452

points[3] = map(5.000000000000000e-01, 5.000000000000000e-01);

1453

points[4] = map(0.000000000000000e+00, 5.000000000000000e-01);

1454

points[5] = map(5.000000000000000e-01, 0.000000000000000e+00);

1455

points[6] = map(0.000000000000000e+00, 0.000000000000000e+00);

1456

points[7] = map(1.000000000000000e+00, 0.000000000000000e+00);

1457

points[8] = map(0.000000000000000e+00, 1.000000000000000e+00);

1458

points[9] = map(5.000000000000000e-01, 5.000000000000000e-01);

1459

points[10] = map(0.000000000000000e+00, 5.000000000000000e-01);

1460

points[11] = map(5.000000000000000e-01, 0.000000000000000e+00);

1461

components[0] = 0;

1462

components[1] = 0;

1463

components[2] = 0;

1464

components[3] = 0;

1465

components[4] = 0;

1466

components[5] = 0;

1467

components[6] = 1;

1468

components[7] = 1;

1469

components[8] = 1;

1470

components[9] = 1;

1471

components[10] = 1;

1472

components[11] = 1;

1473

}

1474

1475

void vertexeval(real values[], unsigned int vertex, const real x[], const Mesh& mesh) const

1476

{

1477

// FIXME: Temporary fix for Lagrange elements

1478

values[0] = x[vertex];

1479

int offset = mesh.numVertices() + mesh.numEdges();

1480

values[1] = x[offset + vertex];

1481

}

1482

1483

const FiniteElement& operator[] (unsigned int i) const

1484

{

1485

return *this;

1486

}

1487

1488

FiniteElement& operator[] (unsigned int i)

1489

{

1490

return *this;

1491

}

1492

1493

FiniteElementSpec spec() const

1494

{

1495

FiniteElementSpec s("Vector Lagrange", "triangle", 2, 2);

1496

return s;

1497

}

1498

1499

private:

1500

1501

unsigned int* tensordims;

1502

FiniteElement** subelements;

1503

1504

};

1505

1506

LinearForm(Function& w0) : dolfin::LinearForm(1)

1507

{

1508

// Create finite element for test space

1509

_test = new TestElement();

1510

1511

// Add functions

1512

add(w0, new FunctionElement_0());

1513

}

1514

1515

void eval(real block[], const AffineMap& map) const

1516

{

1517

// Compute coefficients

1518

const real c0_0 = c[0][0];

1519

const real c0_1 = c[0][1];

1520

const real c0_2 = c[0][2];

1521

const real c0_3 = c[0][3];

1522

const real c0_4 = c[0][4];

1523

const real c0_5 = c[0][5];

1524

const real c0_6 = c[0][6];

1525

const real c0_7 = c[0][7];

1526

const real c0_8 = c[0][8];

1527

const real c0_9 = c[0][9];

1528

const real c0_10 = c[0][10];

1529

const real c0_11 = c[0][11];

1530

1531

// Compute geometry tensors

1532

const real G0_0 = map.det*c0_0;

1533

const real G0_1 = map.det*c0_1;

1534

const real G0_2 = map.det*c0_2;

1535

const real G0_3 = map.det*c0_3;

1536

const real G0_4 = map.det*c0_4;

1537

const real G0_5 = map.det*c0_5;

1538

const real G1_6 = map.det*c0_6;

1539

const real G1_7 = map.det*c0_7;

1540

const real G1_8 = map.det*c0_8;

1541

const real G1_9 = map.det*c0_9;

1542

const real G1_10 = map.det*c0_10;

1543

const real G1_11 = map.det*c0_11;

1544

1545

// Compute element tensor

1546

block[0] = 1.666666666666665e-02*G0_0 - 2.777777777777774e-03*G0_1 - 2.777777777777775e-03*G0_2 - 1.111111111111110e-02*G0_3;

1547

block[1] = -2.777777777777774e-03*G0_0 + 1.666666666666665e-02*G0_1 - 2.777777777777776e-03*G0_2 - 1.111111111111111e-02*G0_4;

1548

block[2] = -2.777777777777775e-03*G0_0 - 2.777777777777776e-03*G0_1 + 1.666666666666666e-02*G0_2 - 1.111111111111111e-02*G0_5;

1549

block[3] = -1.111111111111110e-02*G0_0 + 8.888888888888882e-02*G0_3 + 4.444444444444443e-02*G0_4 + 4.444444444444443e-02*G0_5;

1550

block[4] = -1.111111111111111e-02*G0_1 + 4.444444444444443e-02*G0_3 + 8.888888888888884e-02*G0_4 + 4.444444444444442e-02*G0_5;

1551

block[5] = -1.111111111111111e-02*G0_2 + 4.444444444444443e-02*G0_3 + 4.444444444444443e-02*G0_4 + 8.888888888888882e-02*G0_5;

1552

block[6] = 1.666666666666665e-02*G1_6 - 2.777777777777774e-03*G1_7 - 2.777777777777774e-03*G1_8 - 1.111111111111109e-02*G1_9;

1553

block[7] = -2.777777777777774e-03*G1_6 + 1.666666666666665e-02*G1_7 - 2.777777777777775e-03*G1_8 - 1.111111111111111e-02*G1_10;

1554

block[8] = -2.777777777777774e-03*G1_6 - 2.777777777777775e-03*G1_7 + 1.666666666666666e-02*G1_8 - 1.111111111111111e-02*G1_11;

1555

block[9] = -1.111111111111109e-02*G1_6 + 8.888888888888882e-02*G1_9 + 4.444444444444443e-02*G1_10 + 4.444444444444443e-02*G1_11;

1556

block[10] = -1.111111111111111e-02*G1_7 + 4.444444444444443e-02*G1_9 + 8.888888888888884e-02*G1_10 + 4.444444444444442e-02*G1_11;

1557

block[11] = -1.111111111111111e-02*G1_8 + 4.444444444444443e-02*G1_9 + 4.444444444444443e-02*G1_10 + 8.888888888888882e-02*G1_11;

1558

block[12] = 0.000000000000000e+00;

1559

block[13] = 0.000000000000000e+00;

1560

block[14] = 0.000000000000000e+00;

1561

}

1562

1563

};

1049

1050

class TestElement;

1051

1052

class FunctionElement_0;

1053

1054

LinearForm(Function& w0);

1055

1056

1057

bool interior_contribution() const;

1058

1059

void eval(real block[], const AffineMap& map, real det) const;

1060

1061

bool boundary_contribution() const;

1062

1063

void eval(real block[], const AffineMap& map, real det, unsigned int facet) const;

1064

1065

bool interior_boundary_contribution() const;

1066

1067

void eval(real block[], const AffineMap& map0, const AffineMap& map1, real det, unsigned int facet0, unsigned int facet1, unsigned int alignment) const;

1068

1069

};

1070

1071

class LinearForm::TestElement : public dolfin::FiniteElement

1072

{

1073

public:

1074

1075

TestElement() : dolfin::FiniteElement(), tensordims(0), subelements(0)

1076

{

1077

tensordims = new unsigned int [1];

1078

tensordims[0] = 3;

1079

1080

subelements = new FiniteElement* [2];

1081

subelements[0] = new SubElement_0();

1082

subelements[1] = new SubElement_1();

1083

}

1084

1085

~TestElement()

1086

{

1087

if ( tensordims ) delete [] tensordims;

1088

if ( subelements )

1089

{

1090

for (unsigned int i = 0; i < elementdim(); i++)

1091

delete subelements[i];

1092

delete [] subelements;

1093

}

1094

}

1095

1096

inline unsigned int spacedim() const

1097

{

1098

return 15;

1099

}

1100

1101

inline unsigned int shapedim() const

1102

{

1103

return 2;

1104

}

1105

1106

inline unsigned int tensordim(unsigned int i) const

1107

{

1108

dolfin_assert(i < 1);

1109

return tensordims[i];

1110

}

1111

1112

inline unsigned int elementdim() const

1113

{

1114

return 2;

1115

}

1116

1117

inline unsigned int rank() const

1118

{

1119

return 1;

1120

}

1121

1122

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

1123

{

1124

nodes[0] = cell.entities(0)[0];

1125

nodes[1] = cell.entities(0)[1];

1126

nodes[2] = cell.entities(0)[2];

1127

int offset = mesh.topology().size(0);

1128

nodes[3] = offset + cell.entities(1)[0];

1129

nodes[4] = offset + cell.entities(1)[1];

1130

nodes[5] = offset + cell.entities(1)[2];

1131

offset = offset + mesh.topology().size(1);

1132

nodes[6] = offset + cell.entities(0)[0];

1133

nodes[7] = offset + cell.entities(0)[1];

1134

nodes[8] = offset + cell.entities(0)[2];

1135

offset = offset + mesh.topology().size(0);

1136

nodes[9] = offset + cell.entities(1)[0];

1137

nodes[10] = offset + cell.entities(1)[1];

1138

nodes[11] = offset + cell.entities(1)[2];

1139

offset = offset + mesh.topology().size(1);

1140

nodes[12] = offset + cell.entities(0)[0];

1141

nodes[13] = offset + cell.entities(0)[1];

1142

nodes[14] = offset + cell.entities(0)[2];

1143

}

1144

1145

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

1146

{

1147

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

1148

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

1149

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

1150

points[3] = map(5.000000000000000e-01, 5.000000000000000e-01);

1151

points[4] = map(0.000000000000000e+00, 5.000000000000000e-01);

1152

points[5] = map(5.000000000000000e-01, 0.000000000000000e+00);

1153

points[6] = map(0.000000000000000e+00, 0.000000000000000e+00);

1154

points[7] = map(1.000000000000000e+00, 0.000000000000000e+00);

1155

points[8] = map(0.000000000000000e+00, 1.000000000000000e+00);

1156

points[9] = map(5.000000000000000e-01, 5.000000000000000e-01);

1157

points[10] = map(0.000000000000000e+00, 5.000000000000000e-01);

1158

points[11] = map(5.000000000000000e-01, 0.000000000000000e+00);

1159

components[0] = 0;

1160

components[1] = 0;

1161

components[2] = 0;

1162

components[3] = 0;

1163

components[4] = 0;

1164

components[5] = 0;

1165

components[6] = 1;

1166

components[7] = 1;

1167

components[8] = 1;

1168

components[9] = 1;

1169

components[10] = 1;

1170

components[11] = 1;

1171

points[12] = map(0.000000000000000e+00, 0.000000000000000e+00);

1172

points[13] = map(1.000000000000000e+00, 0.000000000000000e+00);

1173

points[14] = map(0.000000000000000e+00, 1.000000000000000e+00);

1174

components[12] = 2;

1175

components[13] = 2;

1176

components[14] = 2;

1177

}

1178

1179

void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const

1180

{

1181

// FIXME: Temporary fix for Lagrange elements

1182

vertex_nodes[0] = vertex;

1183

int offset = mesh.topology().size(0) + mesh.topology().size(1);

1184

vertex_nodes[1] = offset + vertex;

1185

offset = offset + mesh.topology().size(0) + mesh.topology().size(1);

1186

vertex_nodes[2] = offset + vertex;

1187

}

1188

1189

const FiniteElement& operator[] (unsigned int i) const

1190

{

1191

return *subelements[i];

1192

}

1193

1194

FiniteElement& operator[] (unsigned int i)

1195

{

1196

return *subelements[i];

1197

}

1198

1199

FiniteElementSpec spec() const

1200

{

1201

FiniteElementSpec s("mixed");

1202

return s;

1203

}

1204

1205

private:

1206

1207

class SubElement_0 : public dolfin::FiniteElement

1208

{

1209

public:

1210

1211

SubElement_0() : dolfin::FiniteElement(), tensordims(0), subelements(0)

1212

{

1213

tensordims = new unsigned int [1];

1214

tensordims[0] = 2;

1215

1216

// Element is simple, don't need to initialize subelements

1217

}

1218

1219

~SubElement_0()

1220

{

1221

if ( tensordims ) delete [] tensordims;

1222

if ( subelements )

1223

{

1224

for (unsigned int i = 0; i < elementdim(); i++)

1225

delete subelements[i];

1226

delete [] subelements;

1227

}

1228

}

1229

1230

inline unsigned int spacedim() const

1231

{

1232

return 12;

1233

}

1234

1235

inline unsigned int shapedim() const

1236

{

1237

return 2;

1238

}

1239

1240

inline unsigned int tensordim(unsigned int i) const

1241

{

1242

dolfin_assert(i < 1);

1243

return tensordims[i];

1244

}

1245

1246

inline unsigned int elementdim() const

1247

{

1248

return 1;

1249

}

1250

1251

inline unsigned int rank() const

1252

{

1253

return 1;

1254

}

1255

1256

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

1257

{

1258

nodes[0] = cell.entities(0)[0];

1259

nodes[1] = cell.entities(0)[1];

1260

nodes[2] = cell.entities(0)[2];

1261

int offset = mesh.topology().size(0);

1262

nodes[3] = offset + cell.entities(1)[0];

1263

nodes[4] = offset + cell.entities(1)[1];

1264

nodes[5] = offset + cell.entities(1)[2];

1265

offset = offset + mesh.topology().size(1);

1266

nodes[6] = offset + cell.entities(0)[0];

1267

nodes[7] = offset + cell.entities(0)[1];

1268

nodes[8] = offset + cell.entities(0)[2];

1269

offset = offset + mesh.topology().size(0);

1270

nodes[9] = offset + cell.entities(1)[0];

1271

nodes[10] = offset + cell.entities(1)[1];

1272

nodes[11] = offset + cell.entities(1)[2];

1273

}

1274

1275

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

1276

{

1277

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

1278

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

1279

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

1280

points[3] = map(5.000000000000000e-01, 5.000000000000000e-01);

1281

points[4] = map(0.000000000000000e+00, 5.000000000000000e-01);

1282

points[5] = map(5.000000000000000e-01, 0.000000000000000e+00);

1283

points[6] = map(0.000000000000000e+00, 0.000000000000000e+00);

1284

points[7] = map(1.000000000000000e+00, 0.000000000000000e+00);

1285

points[8] = map(0.000000000000000e+00, 1.000000000000000e+00);

1286

points[9] = map(5.000000000000000e-01, 5.000000000000000e-01);

1287

points[10] = map(0.000000000000000e+00, 5.000000000000000e-01);

1288

points[11] = map(5.000000000000000e-01, 0.000000000000000e+00);

1289

components[0] = 0;

1290

components[1] = 0;

1291

components[2] = 0;

1292

components[3] = 0;

1293

components[4] = 0;

1294

components[5] = 0;

1295

components[6] = 1;

1296

components[7] = 1;

1297

components[8] = 1;

1298

components[9] = 1;

1299

components[10] = 1;

1300

components[11] = 1;

1301

}

1302

1303

void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const

1304

{

1305

// FIXME: Temporary fix for Lagrange elements

1306

vertex_nodes[0] = vertex;

1307

int offset = mesh.topology().size(0) + mesh.topology().size(1);

1308

vertex_nodes[1] = offset + vertex;

1309

}

1310

1311

const FiniteElement& operator[] (unsigned int i) const

1312

{

1313

return *this;

1314

}

1315

1316

FiniteElement& operator[] (unsigned int i)

1317

{

1318

return *this;

1319

}

1320

1321

FiniteElementSpec spec() const

1322

{

1323

FiniteElementSpec s("Vector Lagrange", "triangle", 2, 2);

1324

return s;

1325

}

1326

1327

private:

1328

1329

unsigned int* tensordims;

1330

FiniteElement** subelements;

1331

1332

};

1333

1334

class SubElement_1 : public dolfin::FiniteElement

1335

{

1336

public:

1337

1338

SubElement_1() : dolfin::FiniteElement(), tensordims(0), subelements(0)

1339

{

1340

// Element is scalar, don't need to initialize tensordims

1341

1342

// Element is simple, don't need to initialize subelements

1343

}

1344

1345

~SubElement_1()

1346

{

1347

if ( tensordims ) delete [] tensordims;

1348

if ( subelements )

1349

{

1350

for (unsigned int i = 0; i < elementdim(); i++)

1351

delete subelements[i];

1352

delete [] subelements;

1353

}

1354

}

1355

1356

inline unsigned int spacedim() const

1357

{

1358

return 3;

1359

}

1360

1361

inline unsigned int shapedim() const

1362

{

1363

return 2;

1364

}

1365

1366

inline unsigned int tensordim(unsigned int i) const

1367

{

1368

dolfin_error("Element is scalar.");

1369

return 0;

1370

}

1371

1372

inline unsigned int elementdim() const

1373

{

1374

return 1;

1375

}

1376

1377

inline unsigned int rank() const

1378

{

1379

return 0;

1380

}

1381

1382

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

1383

{

1384

nodes[0] = cell.entities(0)[0];

1385

nodes[1] = cell.entities(0)[1];

1386

nodes[2] = cell.entities(0)[2];

1387

}

1388

1389

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

1390

{

1391

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

1392

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

1393

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

1394

components[0] = 0;

1395

components[1] = 0;

1396

components[2] = 0;

1397

}

1398

1399

void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const

1400

{

1401

// FIXME: Temporary fix for Lagrange elements

1402

vertex_nodes[0] = vertex;

1403

}

1404

1405

const FiniteElement& operator[] (unsigned int i) const

1406

{

1407

return *this;

1408

}

1409

1410

FiniteElement& operator[] (unsigned int i)

1411

{

1412

return *this;

1413

}

1414

1415

FiniteElementSpec spec() const

1416

{

1417

FiniteElementSpec s("Lagrange", "triangle", 1);

1418

return s;

1419

}

1420

1421

private:

1422

1423

unsigned int* tensordims;

1424

FiniteElement** subelements;

1425

1426

};

1427

1428

unsigned int* tensordims;

1429

FiniteElement** subelements;

1430

1431

};

1432

1433

class LinearForm::FunctionElement_0 : public dolfin::FiniteElement

1434

{

1435

public:

1436

1437

FunctionElement_0() : dolfin::FiniteElement(), tensordims(0), subelements(0)

1438

{

1439

tensordims = new unsigned int [1];

1440

tensordims[0] = 2;

1441

1442

// Element is simple, don't need to initialize subelements

1443

}

1444

1445

~FunctionElement_0()

1446

{

1447

if ( tensordims ) delete [] tensordims;

1448

if ( subelements )

1449

{

1450

for (unsigned int i = 0; i < elementdim(); i++)

1451

delete subelements[i];

1452

delete [] subelements;

1453

}

1454

}

1455

1456

inline unsigned int spacedim() const

1457

{

1458

return 12;

1459

}

1460

1461

inline unsigned int shapedim() const

1462

{

1463

return 2;

1464

}

1465

1466

inline unsigned int tensordim(unsigned int i) const

1467

{

1468

dolfin_assert(i < 1);

1469

return tensordims[i];

1470

}

1471

1472

inline unsigned int elementdim() const

1473

{

1474

return 1;

1475

}

1476

1477

inline unsigned int rank() const

1478

{

1479

return 1;

1480

}

1481

1482

void nodemap(int nodes[], const Cell& cell, const Mesh& mesh) const

1483

{

1484

nodes[0] = cell.entities(0)[0];

1485

nodes[1] = cell.entities(0)[1];

1486

nodes[2] = cell.entities(0)[2];

1487

int offset = mesh.topology().size(0);

1488

nodes[3] = offset + cell.entities(1)[0];

1489

nodes[4] = offset + cell.entities(1)[1];

1490

nodes[5] = offset + cell.entities(1)[2];

1491

offset = offset + mesh.topology().size(1);

1492

nodes[6] = offset + cell.entities(0)[0];

1493

nodes[7] = offset + cell.entities(0)[1];

1494

nodes[8] = offset + cell.entities(0)[2];

1495

offset = offset + mesh.topology().size(0);

1496

nodes[9] = offset + cell.entities(1)[0];

1497

nodes[10] = offset + cell.entities(1)[1];

1498

nodes[11] = offset + cell.entities(1)[2];

1499

}

1500

1501

void pointmap(Point points[], unsigned int components[], const AffineMap& map) const

1502

{

1503

points[0] = map(0.000000000000000e+00, 0.000000000000000e+00);

1504

points[1] = map(1.000000000000000e+00, 0.000000000000000e+00);

1505

points[2] = map(0.000000000000000e+00, 1.000000000000000e+00);

1506

points[3] = map(5.000000000000000e-01, 5.000000000000000e-01);

1507

points[4] = map(0.000000000000000e+00, 5.000000000000000e-01);

1508

points[5] = map(5.000000000000000e-01, 0.000000000000000e+00);

1509

points[6] = map(0.000000000000000e+00, 0.000000000000000e+00);

1510

points[7] = map(1.000000000000000e+00, 0.000000000000000e+00);

1511

points[8] = map(0.000000000000000e+00, 1.000000000000000e+00);

1512

points[9] = map(5.000000000000000e-01, 5.000000000000000e-01);

1513

points[10] = map(0.000000000000000e+00, 5.000000000000000e-01);

1514

points[11] = map(5.000000000000000e-01, 0.000000000000000e+00);

1515

components[0] = 0;

1516

components[1] = 0;

1517

components[2] = 0;

1518

components[3] = 0;

1519

components[4] = 0;

1520

components[5] = 0;

1521

components[6] = 1;

1522

components[7] = 1;

1523

components[8] = 1;

1524

components[9] = 1;

1525

components[10] = 1;

1526

components[11] = 1;

1527

}

1528

1529

void vertexeval(uint vertex_nodes[], unsigned int vertex, const Mesh& mesh) const

1530

{

1531

// FIXME: Temporary fix for Lagrange elements

1532

vertex_nodes[0] = vertex;

1533

int offset = mesh.topology().size(0) + mesh.topology().size(1);

1534

vertex_nodes[1] = offset + vertex;

1535

}

1536

1537

const FiniteElement& operator[] (unsigned int i) const

1538

{

1539

return *this;

1540

}

1541

1542

FiniteElement& operator[] (unsigned int i)

1543

{

1544

return *this;

1545

}

1546

1547

FiniteElementSpec spec() const

1548

{

1549

FiniteElementSpec s("Vector Lagrange", "triangle", 2, 2);

1550

return s;

1551

}

1552

1553

private:

1554

1555

unsigned int* tensordims;

1556

FiniteElement** subelements;

1557

1558

};

1559

1560

LinearForm::LinearForm(Function& w0) : dolfin::LinearForm(1)

1561

{

1562

// Create finite element for test space

1563

_test = new TestElement();

1564

1565

// Add functions

1566

initFunction(0, w0, new FunctionElement_0());

1567

}

1568

1569

// Contribution from the interior

1570

bool LinearForm::interior_contribution() const { return true; }

1571

1572

void LinearForm::eval(real block[], const AffineMap& map, real det) const

1573

{

1574

// Compute coefficients

1575

const real c0_0 = c[0][0];

1576

const real c0_1 = c[0][1];

1577

const real c0_2 = c[0][2];

1578

const real c0_3 = c[0][3];

1579

const real c0_4 = c[0][4];

1580

const real c0_5 = c[0][5];

1581

const real c0_6 = c[0][6];

1582

const real c0_7 = c[0][7];

1583

const real c0_8 = c[0][8];

1584

const real c0_9 = c[0][9];

1585

const real c0_10 = c[0][10];

1586

const real c0_11 = c[0][11];

1587

1588

// Compute geometry tensors

1589

const real G0_0 = det*c0_0;

1590

const real G0_1 = det*c0_1;

1591

const real G0_2 = det*c0_2;

1592

const real G0_3 = det*c0_3;

1593

const real G0_4 = det*c0_4;

1594

const real G0_5 = det*c0_5;

1595

const real G1_6 = det*c0_6;

1596

const real G1_7 = det*c0_7;

1597

const real G1_8 = det*c0_8;

1598

const real G1_9 = det*c0_9;

1599

const real G1_10 = det*c0_10;

1600

const real G1_11 = det*c0_11;

1601

1602

// Compute element tensor

1603

block[0] = 1.666666666666665e-02*G0_0 - 2.777777777777774e-03*G0_1 - 2.777777777777775e-03*G0_2 - 1.111111111111110e-02*G0_3;

1604

block[1] = -2.777777777777774e-03*G0_0 + 1.666666666666665e-02*G0_1 - 2.777777777777776e-03*G0_2 - 1.111111111111111e-02*G0_4;

1605

block[2] = -2.777777777777775e-03*G0_0 - 2.777777777777776e-03*G0_1 + 1.666666666666666e-02*G0_2 - 1.111111111111111e-02*G0_5;

1606

block[3] = -1.111111111111110e-02*G0_0 + 8.888888888888882e-02*G0_3 + 4.444444444444443e-02*G0_4 + 4.444444444444443e-02*G0_5;

1607

block[4] = -1.111111111111111e-02*G0_1 + 4.444444444444443e-02*G0_3 + 8.888888888888884e-02*G0_4 + 4.444444444444442e-02*G0_5;

1608

block[5] = -1.111111111111111e-02*G0_2 + 4.444444444444443e-02*G0_3 + 4.444444444444443e-02*G0_4 + 8.888888888888882e-02*G0_5;

1609

block[6] = 1.666666666666665e-02*G1_6 - 2.777777777777774e-03*G1_7 - 2.777777777777774e-03*G1_8 - 1.111111111111109e-02*G1_9;

1610

block[7] = -2.777777777777774e-03*G1_6 + 1.666666666666665e-02*G1_7 - 2.777777777777775e-03*G1_8 - 1.111111111111111e-02*G1_10;

1611

block[8] = -2.777777777777774e-03*G1_6 - 2.777777777777775e-03*G1_7 + 1.666666666666666e-02*G1_8 - 1.111111111111111e-02*G1_11;

1612

block[9] = -1.111111111111109e-02*G1_6 + 8.888888888888882e-02*G1_9 + 4.444444444444443e-02*G1_10 + 4.444444444444443e-02*G1_11;

1613

block[10] = -1.111111111111111e-02*G1_7 + 4.444444444444443e-02*G1_9 + 8.888888888888884e-02*G1_10 + 4.444444444444442e-02*G1_11;

1614

block[11] = -1.111111111111111e-02*G1_8 + 4.444444444444443e-02*G1_9 + 4.444444444444443e-02*G1_10 + 8.888888888888882e-02*G1_11;

1615

block[12] = 0.000000000000000e+00;

1616

block[13] = 0.000000000000000e+00;

1617

block[14] = 0.000000000000000e+00;

1618

}

1619

1620

// No contribution from the boundary

1621

bool LinearForm::boundary_contribution() const { return false; }

1622

1623

void LinearForm::eval(real block[], const AffineMap& map, real det, unsigned int facet) const {}

1624

1625

// No contribution from interior boundaries

1626

bool LinearForm::interior_boundary_contribution() const { return false; }

1627

1628

void LinearForm::eval(real block[], const AffineMap& map0, const AffineMap& map1, real det, unsigned int facet0, unsigned int facet1, unsigned int alignment) const {}

1564

1629

1565

1630

} }

1566

1631

Older »