~njansson/dolfin/hpc

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];

359

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];

360

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];

361

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];

362

new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];

363

new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];

364

}

365

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

366

{

367

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];

368

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];

369

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];

370

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];

371

new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];

372

new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];

373

}

374

375

}

376

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

377

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;

378

}

379

380

// Transform derivatives back to physical element

381

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

382

{

383

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

384

{

385

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

386

}

387

}

388

// Delete pointer to array of derivatives on FIAT element

389

delete [] derivatives;

390

391

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

392

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

393

{

394

delete [] combinations[row];

395

delete [] transform[row];

396

}

397

398

delete [] combinations;

399

delete [] transform;

400

}

401

402

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

403

virtual void evaluate_basis_derivatives_all(unsigned int n,

404

double* values,

405

const double* coordinates,

406

const ufc::cell& c) const

407

{

408

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

409

}

410

411

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

412

virtual double evaluate_dof(unsigned int i,

413

const ufc::function& f,

414

const ufc::cell& c) const

415

{

416

// The reference points, direction and weights:

417

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

418

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

419

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

420

421

const double * const * x = c.coordinates;

422

double result = 0.0;

423

// Iterate over the points:

424

// Evaluate basis functions for affine mapping

425

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

426

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

427

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

428

429

// Compute affine mapping y = F(X)

430

double y[2];

431

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

432

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

433

434

// Evaluate function at physical points

435

double values[1];

436

f.evaluate(values, y, c);

437

438

// Map function values using appropriate mapping

439

// Affine map: Do nothing

440

441

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

442

443

// Take directional components

444

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

445

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

446

// Multiply by weights

447

result *= W[i][0];

448

449

return result;

450

}

451

452

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

453

virtual void evaluate_dofs(double* values,

454

const ufc::function& f,

455

const ufc::cell& c) const

456

{

457

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

458

}

459

460

/// Interpolate vertex values from dof values

461

virtual void interpolate_vertex_values(double* vertex_values,

462

const double* dof_values,

463

const ufc::cell& c) const

464

{

465

// Evaluate at vertices and use affine mapping

466

vertex_values[0] = dof_values[0];

467

vertex_values[1] = dof_values[1];

468

vertex_values[2] = dof_values[2];

469

}

470

471

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

472

virtual unsigned int num_sub_elements() const

473

{

474

return 1;

475

}

476

477

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

478

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

479

{

480

return new UFC_PoissonP2BilinearForm_finite_element_0();

481

}

482

483

};

484

485

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

486

487

class UFC_PoissonP2BilinearForm_finite_element_1: public ufc::finite_element

488

{

489

public:

490

491

/// Constructor

492

UFC_PoissonP2BilinearForm_finite_element_1() : ufc::finite_element()

493

{

494

// Do nothing

495

}

496

497

/// Destructor

498

virtual ~UFC_PoissonP2BilinearForm_finite_element_1()

499

{

500

// Do nothing

501

}

502

503

/// Return a string identifying the finite element

504

virtual const char* signature() const

505

{

506

return "Lagrange finite element of degree 2 on a triangle";

507

}

508

509

/// Return the cell shape

510

virtual ufc::shape cell_shape() const

511

{

512

return ufc::triangle;

513

}

514

515

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

516

virtual unsigned int space_dimension() const

517

{

518

return 6;

519

}

520

521

/// Return the rank of the value space

522

virtual unsigned int value_rank() const

523

{

524

return 0;

525

}

526

527

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

528

virtual unsigned int value_dimension(unsigned int i) const

529

{

530

return 1;

531

}

532

533

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

534

virtual void evaluate_basis(unsigned int i,

535

double* values,

536

const double* coordinates,

537

const ufc::cell& c) const

538

{

539

// Extract vertex coordinates

540

const double * const * element_coordinates = c.coordinates;

541

542

// Compute Jacobian of affine map from reference cell

543

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

544

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

545

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

546

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

547

548

// Compute determinant of Jacobian

549

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

550

551

// Compute inverse of Jacobian

552

553

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

554

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

555

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

556

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

557

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

558

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

559

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

560

561

// Map coordinates to the reference square

562

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

563

x = -1.0;

564

else

565

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

566

y = 2.0*y - 1.0;

567

568

// Reset values

569

*values = 0;

570

571

// Map degree of freedom to element degree of freedom

572

const unsigned int dof = i;

573

574

// Generate scalings

575

const double scalings_y_0 = 1;

576

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

577

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

578

579

// Compute psitilde_a

580

const double psitilde_a_0 = 1;

581

const double psitilde_a_1 = x;

582

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

583

584

// Compute psitilde_bs

585

const double psitilde_bs_0_0 = 1;

586

const double psitilde_bs_0_1 = 1.5*y + 0.5;

587

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

588

const double psitilde_bs_1_0 = 1;

589

const double psitilde_bs_1_1 = 2.5*y + 1.5;

590

const double psitilde_bs_2_0 = 1;

591

592

// Compute basisvalues

593

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

594

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

595

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

596

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

597

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

598

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

599

600

// Table(s) of coefficients

601

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

602

{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817},

603

{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818},

604

{0, 0, 0.2, 0, 0, 0.163299316185545},

605

{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},

606

{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},

607

{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};

608

609

// Extract relevant coefficients

610

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

611

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

612

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

613

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

614

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

615

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

616

617

// Compute value(s)

618

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;

619

}

620

621

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

622

virtual void evaluate_basis_all(double* values,

623

const double* coordinates,

624

const ufc::cell& c) const

625

{

626

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

627

}

628

629

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

630

virtual void evaluate_basis_derivatives(unsigned int i,

631

unsigned int n,

632

double* values,

633

const double* coordinates,

634

const ufc::cell& c) const

635

{

636

// Extract vertex coordinates

637

const double * const * element_coordinates = c.coordinates;

638

639

// Compute Jacobian of affine map from reference cell

640

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

641

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

642

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

643

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

644

645

// Compute determinant of Jacobian

646

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

647

648

// Compute inverse of Jacobian

649

650

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

651

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

652

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

653

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

654

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

655

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

656

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

657

658

// Map coordinates to the reference square

659

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

660

x = -1.0;

661

else

662

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

663

y = 2.0*y - 1.0;

664

665

// Compute number of derivatives

666

unsigned int num_derivatives = 1;

667

668

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

669

num_derivatives *= 2;

670

671

672

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

673

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

674

675

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

676

{

677

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

678

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

679

combinations[j][k] = 0;

680

}

681

682

// Generate combinations of derivatives

683

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

684

{

685

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

686

{

687

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

688

{

689

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

690

combinations[row][col] = 0;

691

else

692

{

693

combinations[row][col] += 1;

694

break;

695

}

696

}

697

}

698

}

699

700

// Compute inverse of Jacobian

701

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

702

703

// Declare transformation matrix

704

// Declare pointer to two dimensional array and initialise

705

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

706

707

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

708

{

709

transform[j] = new double [num_derivatives];

710

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

711

transform[j][k] = 1;

712

}

713

714

// Construct transformation matrix

715

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

716

{

717

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

718

{

719

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

720

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

721

}

722

}

723

724

// Reset values

725

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

726

values[j] = 0;

727

728

// Map degree of freedom to element degree of freedom

729

const unsigned int dof = i;

730

731

// Generate scalings

732

const double scalings_y_0 = 1;

733

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

734

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

735

736

// Compute psitilde_a

737

const double psitilde_a_0 = 1;

738

const double psitilde_a_1 = x;

739

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

740

741

// Compute psitilde_bs

742

const double psitilde_bs_0_0 = 1;

743

const double psitilde_bs_0_1 = 1.5*y + 0.5;

744

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

745

const double psitilde_bs_1_0 = 1;

746

const double psitilde_bs_1_1 = 2.5*y + 1.5;

747

const double psitilde_bs_2_0 = 1;

748

749

// Compute basisvalues

750

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

751

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

752

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

753

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

754

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

755

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

756

757

// Table(s) of coefficients

758

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

759

{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817},

760

{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818},

761

{0, 0, 0.2, 0, 0, 0.163299316185545},

762

{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},

763

{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},

764

{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};

765

766

// Interesting (new) part

767

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

768

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

769

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

770

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

771

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

772

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

773

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

774

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

775

776

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

777

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

778

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

779

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

780

{2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},

781

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

782

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

783

784

// Compute reference derivatives

785

// Declare pointer to array of derivatives on FIAT element

786

double *derivatives = new double [num_derivatives];

787

788

// Declare coefficients

789

double coeff0_0 = 0;

790

double coeff0_1 = 0;

791

double coeff0_2 = 0;

792

double coeff0_3 = 0;

793

double coeff0_4 = 0;

794

double coeff0_5 = 0;

795

796

// Declare new coefficients

797

double new_coeff0_0 = 0;

798

double new_coeff0_1 = 0;

799

double new_coeff0_2 = 0;

800

double new_coeff0_3 = 0;

801

double new_coeff0_4 = 0;

802

double new_coeff0_5 = 0;

803

804

// Loop possible derivatives

805

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

806

{

807

// Get values from coefficients array

808

new_coeff0_0 = coefficients0[dof][0];

809

new_coeff0_1 = coefficients0[dof][1];

810

new_coeff0_2 = coefficients0[dof][2];

811

new_coeff0_3 = coefficients0[dof][3];

812

new_coeff0_4 = coefficients0[dof][4];

813

new_coeff0_5 = coefficients0[dof][5];

814

815

// Loop derivative order

816

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

817

{

818

// Update old coefficients

819

coeff0_0 = new_coeff0_0;

820

coeff0_1 = new_coeff0_1;

821

coeff0_2 = new_coeff0_2;

822

coeff0_3 = new_coeff0_3;

823

coeff0_4 = new_coeff0_4;

824

coeff0_5 = new_coeff0_5;

825

826

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

827

{

828

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];

829

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];

830

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];

831

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];

832

new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];

833

new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];

834

}

835

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

836

{

837

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];

838

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];

839

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];

840

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];

841

new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];

842

new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];

843

}

844

845

}

846

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

847

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;

848

}

849

850

// Transform derivatives back to physical element

851

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

852

{

853

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

854

{

855

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

856

}

857

}

858

// Delete pointer to array of derivatives on FIAT element

859

delete [] derivatives;

860

861

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

862

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

863

{

864

delete [] combinations[row];

865

delete [] transform[row];

866

}

867

868

delete [] combinations;

869

delete [] transform;

870

}

871

872

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

873

virtual void evaluate_basis_derivatives_all(unsigned int n,

874

double* values,

875

const double* coordinates,

876

const ufc::cell& c) const

877

{

878

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

879

}

880

881

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

882

virtual double evaluate_dof(unsigned int i,

883

const ufc::function& f,

884

const ufc::cell& c) const

885

{

886

// The reference points, direction and weights:

887

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

888

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

889

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

890

891

const double * const * x = c.coordinates;

892

double result = 0.0;

893

// Iterate over the points:

894

// Evaluate basis functions for affine mapping

895

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

896

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

897

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

898

899

// Compute affine mapping y = F(X)

900

double y[2];

901

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

902

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

903

904

// Evaluate function at physical points

905

double values[1];

906

f.evaluate(values, y, c);

907

908

// Map function values using appropriate mapping

909

// Affine map: Do nothing

910

911

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

912

913

// Take directional components

914

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

915

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

916

// Multiply by weights

917

result *= W[i][0];

918

919

return result;

920

}

921

922

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

923

virtual void evaluate_dofs(double* values,

924

const ufc::function& f,

925

const ufc::cell& c) const

926

{

927

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

928

}

929

930

/// Interpolate vertex values from dof values

931

virtual void interpolate_vertex_values(double* vertex_values,

932

const double* dof_values,

933

const ufc::cell& c) const

934

{

935

// Evaluate at vertices and use affine mapping

936

vertex_values[0] = dof_values[0];

937

vertex_values[1] = dof_values[1];

938

vertex_values[2] = dof_values[2];

939

}

940

941

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

942

virtual unsigned int num_sub_elements() const

943

{

944

return 1;

945

}

946

947

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

948

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

949

{

950

return new UFC_PoissonP2BilinearForm_finite_element_1();

951

}

952

953

};

954

955

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

956

/// degrees of freedom (dofs).

957

958

class UFC_PoissonP2BilinearForm_dof_map_0: public ufc::dof_map

959

{

960

private:

961

962

unsigned int __global_dimension;

963

964

public:

965

966

/// Constructor

967

UFC_PoissonP2BilinearForm_dof_map_0() : ufc::dof_map()

968

{

969

__global_dimension = 0;

970

}

971

972

/// Destructor

973

virtual ~UFC_PoissonP2BilinearForm_dof_map_0()

974

{

975

// Do nothing

976

}

977

978

/// Return a string identifying the dof map

979

virtual const char* signature() const

980

{

981

return "FFC dof map for Lagrange finite element of degree 2 on a triangle";

982

}

983

984

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

985

virtual bool needs_mesh_entities(unsigned int d) const

986

{

987

switch ( d )

988

{

989

case 0:

990

return true;

991

break;

992

case 1:

993

return true;

994

break;

995

case 2:

996

return false;

997

break;

998

}

999

return false;

1000

}

1001

1002

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

1003

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

1004

{

1005

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

1006

return false;

1007

}

1008

1009

/// Initialize dof map for given cell

1010

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

1011

const ufc::cell& c)

1012

{

1013

// Do nothing

1014

}

1015

1016

/// Finish initialization of dof map for cells

1017

virtual void init_cell_finalize()

1018

{

1019

// Do nothing

1020

}

1021

1022

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

1023

virtual unsigned int global_dimension() const

1024

{

1025

return __global_dimension;

1026

}

1027

1028

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

1029

virtual unsigned int local_dimension() const

1030

{

1031

return 6;

1032

}

1033

1034

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

1035

virtual unsigned int geometric_dimension() const

1036

{

1037

return 2;

1038

}

1039

1040

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

1041

virtual unsigned int num_facet_dofs() const

1042

{

1043

return 3;

1044

}

1045

1046

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

1047

virtual unsigned int num_entity_dofs(unsigned int d) const

1048

{

1049

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

1050

}

1051

1052

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

1053

virtual void tabulate_dofs(unsigned int* dofs,

1054

const ufc::mesh& m,

1055

const ufc::cell& c) const

1056

{

1057

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

1058

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

1059

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

1060

unsigned int offset = m.num_entities[0];

1061

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

1062

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

1063

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

1064

}

1065

1066

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

1067

virtual void tabulate_facet_dofs(unsigned int* dofs,

1068

unsigned int facet) const

1069

{

1070

switch ( facet )

1071

{

1072

case 0:

1073

dofs[0] = 1;

1074

dofs[1] = 2;

1075

dofs[2] = 3;

1076

break;

1077

case 1:

1078

dofs[0] = 0;

1079

dofs[1] = 2;

1080

dofs[2] = 4;

1081

break;

1082

case 2:

1083

dofs[0] = 0;

1084

dofs[1] = 1;

1085

dofs[2] = 5;

1086

break;

1087

}

1088

}

1089

1090

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

1091

virtual void tabulate_entity_dofs(unsigned int* dofs,

1092

unsigned int d, unsigned int i) const

1093

{

1094

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

1095

}

1096

1097

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

1098

virtual void tabulate_coordinates(double** coordinates,

1099

const ufc::cell& c) const

1100

{

1101

const double * const * x = c.coordinates;

1102

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

1103

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

1104

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

1105

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

1106

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

1107

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

1108

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

1109

coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];

1110

coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];

1111

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

1112

coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];

1113

coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];

1114

}

1115

1116

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

1117

virtual unsigned int num_sub_dof_maps() const

1118

{

1119

return 1;

1120

}

1121

1122

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

1123

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

1124

{

1125

return new UFC_PoissonP2BilinearForm_dof_map_0();

1126

}

1127

1128

};

1129

1130

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

1131

/// degrees of freedom (dofs).

1132

1133

class UFC_PoissonP2BilinearForm_dof_map_1: public ufc::dof_map

1134

{

1135

private:

1136

1137

unsigned int __global_dimension;

1138

1139

public:

1140

1141

/// Constructor

1142

UFC_PoissonP2BilinearForm_dof_map_1() : ufc::dof_map()

1143

{

1144

__global_dimension = 0;

1145

}

1146

1147

/// Destructor

1148

virtual ~UFC_PoissonP2BilinearForm_dof_map_1()

1149

{

1150

// Do nothing

1151

}

1152

1153

/// Return a string identifying the dof map

1154

virtual const char* signature() const

1155

{

1156

return "FFC dof map for Lagrange finite element of degree 2 on a triangle";

1157

}

1158

1159

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

1160

virtual bool needs_mesh_entities(unsigned int d) const

1161

{

1162

switch ( d )

1163

{

1164

case 0:

1165

return true;

1166

break;

1167

case 1:

1168

return true;

1169

break;

1170

case 2:

1171

return false;

1172

break;

1173

}

1174

return false;

1175

}

1176

1177

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

1178

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

1179

{

1180

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

1181

return false;

1182

}

1183

1184

/// Initialize dof map for given cell

1185

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

1186

const ufc::cell& c)

1187

{

1188

// Do nothing

1189

}

1190

1191

/// Finish initialization of dof map for cells

1192

virtual void init_cell_finalize()

1193

{

1194

// Do nothing

1195

}

1196

1197

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

1198

virtual unsigned int global_dimension() const

1199

{

1200

return __global_dimension;

1201

}

1202

1203

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

1204

virtual unsigned int local_dimension() const

1205

{

1206

return 6;

1207

}

1208

1209

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

1210

virtual unsigned int geometric_dimension() const

1211

{

1212

return 2;

1213

}

1214

1215

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

1216

virtual unsigned int num_facet_dofs() const

1217

{

1218

return 3;

1219

}

1220

1221

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

1222

virtual unsigned int num_entity_dofs(unsigned int d) const

1223

{

1224

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

1225

}

1226

1227

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

1228

virtual void tabulate_dofs(unsigned int* dofs,

1229

const ufc::mesh& m,

1230

const ufc::cell& c) const

1231

{

1232

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

1233

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

1234

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

1235

unsigned int offset = m.num_entities[0];

1236

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

1237

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

1238

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

1239

}

1240

1241

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

1242

virtual void tabulate_facet_dofs(unsigned int* dofs,

1243

unsigned int facet) const

1244

{

1245

switch ( facet )

1246

{

1247

case 0:

1248

dofs[0] = 1;

1249

dofs[1] = 2;

1250

dofs[2] = 3;

1251

break;

1252

case 1:

1253

dofs[0] = 0;

1254

dofs[1] = 2;

1255

dofs[2] = 4;

1256

break;

1257

case 2:

1258

dofs[0] = 0;

1259

dofs[1] = 1;

1260

dofs[2] = 5;

1261

break;

1262

}

1263

}

1264

1265

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

1266

virtual void tabulate_entity_dofs(unsigned int* dofs,

1267

unsigned int d, unsigned int i) const

1268

{

1269

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

1270

}

1271

1272

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

1273

virtual void tabulate_coordinates(double** coordinates,

1274

const ufc::cell& c) const

1275

{

1276

const double * const * x = c.coordinates;

1277

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

1278

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

1279

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

1280

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

1281

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

1282

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

1283

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

1284

coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];

1285

coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];

1286

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

1287

coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];

1288

coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];

1289

}

1290

1291

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

1292

virtual unsigned int num_sub_dof_maps() const

1293

{

1294

return 1;

1295

}

1296

1297

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

1298

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

1299

{

1300

return new UFC_PoissonP2BilinearForm_dof_map_1();

1301

}

1302

1303

};

1304

1305

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

1306

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

1307

/// the integral over a cell.

1308

1309

class UFC_PoissonP2BilinearForm_cell_integral_0: public ufc::cell_integral

1310

{

1311

public:

1312

1313

/// Constructor

1314

UFC_PoissonP2BilinearForm_cell_integral_0() : ufc::cell_integral()

1315

{

1316

// Do nothing

1317

}

1318

1319

/// Destructor

1320

virtual ~UFC_PoissonP2BilinearForm_cell_integral_0()

1321

{

1322

// Do nothing

1323

}

1324

1325

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

1326

virtual void tabulate_tensor(double* A,

1327

const double * const * w,

1328

const ufc::cell& c) const

1329

{

1330

// Extract vertex coordinates

1331

const double * const * x = c.coordinates;

1332

1333

// Compute Jacobian of affine map from reference cell

1334

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

1335

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

1336

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

1337

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

1338

1339

// Compute determinant of Jacobian

1340

double detJ = J_00*J_11 - J_01*J_10;

1341

1342

// Compute inverse of Jacobian

1343

const double Jinv_00 = J_11 / detJ;

1344

const double Jinv_01 = -J_01 / detJ;

1345

const double Jinv_10 = -J_10 / detJ;

1346

const double Jinv_11 = J_00 / detJ;

1347

1348

// Set scale factor

1349

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

1350

1351

// Compute geometry tensors

1352

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

1353

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

1354

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

1355

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

1356

1357

// Compute element tensor

1358

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

1359

A[1] = 0.166666666666667*G0_0_0 + 0.166666666666666*G0_1_0;

1360

A[2] = 0.166666666666666*G0_0_1 + 0.166666666666666*G0_1_1;

1361

A[3] = 0;

1362

A[4] = -0.666666666666665*G0_0_1 - 0.666666666666665*G0_1_1;

1363

A[5] = -0.666666666666666*G0_0_0 - 0.666666666666665*G0_1_0;

1364

A[6] = 0.166666666666667*G0_0_0 + 0.166666666666666*G0_0_1;

1365

A[7] = 0.499999999999999*G0_0_0;

1366

A[8] = -0.166666666666666*G0_0_1;

1367

A[9] = 0.666666666666665*G0_0_1;

1368

A[10] = 0;

1369

A[11] = -0.666666666666665*G0_0_0 - 0.666666666666665*G0_0_1;

1370

A[12] = 0.166666666666666*G0_1_0 + 0.166666666666666*G0_1_1;

1371

A[13] = -0.166666666666666*G0_1_0;

1372

A[14] = 0.499999999999999*G0_1_1;

1373

A[15] = 0.666666666666665*G0_1_0;

1374

A[16] = -0.666666666666665*G0_1_0 - 0.666666666666665*G0_1_1;

1375

A[17] = 0;

1376

A[18] = 0;

1377

A[19] = 0.666666666666665*G0_1_0;

1378

A[20] = 0.666666666666665*G0_0_1;

1379

A[21] = 1.33333333333333*G0_0_0 + 0.666666666666665*G0_0_1 + 0.666666666666665*G0_1_0 + 1.33333333333333*G0_1_1;

1380

A[22] = -1.33333333333333*G0_0_0 - 0.666666666666665*G0_0_1 - 0.666666666666665*G0_1_0;

1381

A[23] = -0.666666666666665*G0_0_1 - 0.666666666666665*G0_1_0 - 1.33333333333333*G0_1_1;

1382

A[24] = -0.666666666666665*G0_1_0 - 0.666666666666665*G0_1_1;

1383

A[25] = 0;

1384

A[26] = -0.666666666666665*G0_0_1 - 0.666666666666665*G0_1_1;

1385

A[27] = -1.33333333333333*G0_0_0 - 0.666666666666665*G0_0_1 - 0.666666666666665*G0_1_0;

1386

A[28] = 1.33333333333333*G0_0_0 + 0.666666666666665*G0_0_1 + 0.666666666666665*G0_1_0 + 1.33333333333333*G0_1_1;

1387

A[29] = 0.666666666666665*G0_0_1 + 0.666666666666665*G0_1_0;

1388

A[30] = -0.666666666666666*G0_0_0 - 0.666666666666665*G0_0_1;

1389

A[31] = -0.666666666666666*G0_0_0 - 0.666666666666665*G0_1_0;

1390

A[32] = 0;

1391

A[33] = -0.666666666666665*G0_0_1 - 0.666666666666665*G0_1_0 - 1.33333333333333*G0_1_1;

1392

A[34] = 0.666666666666665*G0_0_1 + 0.666666666666665*G0_1_0;

1393

A[35] = 1.33333333333333*G0_0_0 + 0.666666666666666*G0_0_1 + 0.666666666666666*G0_1_0 + 1.33333333333333*G0_1_1;

1394

}

1395

1396

};

1397

1398

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

1399

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

1400

/// mapping

1401

///

1402

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

1403

///

1404

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

1405

/// global tensor A is defined by

1406

///

1407

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

1408

///

1409

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

1410

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

1411

/// fixed functions (coefficients).

1412

1413

class UFC_PoissonP2BilinearForm: public ufc::form

1414

{

1415

public:

1416

1417

/// Constructor

1418

UFC_PoissonP2BilinearForm() : ufc::form()

1419

{

1420

// Do nothing

1421

}

1422

1423

/// Destructor

1424

virtual ~UFC_PoissonP2BilinearForm()

1425

{

1426

// Do nothing

1427

}

1428

1429

/// Return a string identifying the form

1430

virtual const char* signature() const

1431

{

1432

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

1433

}

1434

1435

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

1436

virtual unsigned int rank() const

1437

{

1438

return 2;

1439

}

1440

1441

/// Return the number of coefficients (n)

1442

virtual unsigned int num_coefficients() const

1443

{

1444

return 0;

1445

}

1446

1447

/// Return the number of cell integrals

1448

virtual unsigned int num_cell_integrals() const

1449

{

1450

return 1;

1451

}

1452

1453

/// Return the number of exterior facet integrals

1454

virtual unsigned int num_exterior_facet_integrals() const

1455

{

1456

return 0;

1457

}

1458

1459

/// Return the number of interior facet integrals

1460

virtual unsigned int num_interior_facet_integrals() const

1461

{

1462

return 0;

1463

}

1464

1465

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

1466

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

1467

{

1468

switch ( i )

1469

{

1470

case 0:

1471

return new UFC_PoissonP2BilinearForm_finite_element_0();

1472

break;

1473

case 1:

1474

return new UFC_PoissonP2BilinearForm_finite_element_1();

1475

break;

1476

}

1477

return 0;

1478

}

1479

1480

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

1481

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

1482

{

1483

switch ( i )

1484

{

1485

case 0:

1486

return new UFC_PoissonP2BilinearForm_dof_map_0();

1487

break;

1488

case 1:

1489

return new UFC_PoissonP2BilinearForm_dof_map_1();

1490

break;

1491

}

1492

return 0;

1493

}

1494

1495

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

1496

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

1497

{

1498

return new UFC_PoissonP2BilinearForm_cell_integral_0();

1499

}

1500

1501

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

1502

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

1503

{

1504

return 0;

1505

}

1506

1507

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

1508

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

1509

{

1510

return 0;

1511

}

1512

1513

};

1514

1515

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

1516

1517

class UFC_PoissonP2LinearForm_finite_element_0: public ufc::finite_element

1518

{

1519

public:

1520

1521

/// Constructor

1522

UFC_PoissonP2LinearForm_finite_element_0() : ufc::finite_element()

1523

{

1524

// Do nothing

1525

}

1526

1527

/// Destructor

1528

virtual ~UFC_PoissonP2LinearForm_finite_element_0()

1529

{

1530

// Do nothing

1531

}

1532

1533

/// Return a string identifying the finite element

1534

virtual const char* signature() const

1535

{

1536

return "Lagrange finite element of degree 2 on a triangle";

1537

}

1538

1539

/// Return the cell shape

1540

virtual ufc::shape cell_shape() const

1541

{

1542

return ufc::triangle;

1543

}

1544

1545

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

1546

virtual unsigned int space_dimension() const

1547

{

1548

return 6;

1549

}

1550

1551

/// Return the rank of the value space

1552

virtual unsigned int value_rank() const

1553

{

1554

return 0;

1555

}

1556

1557

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

1558

virtual unsigned int value_dimension(unsigned int i) const

1559

{

1560

return 1;

1561

}

1562

1563

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

1564

virtual void evaluate_basis(unsigned int i,

1565

double* values,

1566

const double* coordinates,

1567

const ufc::cell& c) const

1568

{

1569

// Extract vertex coordinates

1570

const double * const * element_coordinates = c.coordinates;

1571

1572

// Compute Jacobian of affine map from reference cell

1573

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

1574

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

1575

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

1576

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

1577

1578

// Compute determinant of Jacobian

1579

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

1580

1581

// Compute inverse of Jacobian

1582

1583

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

1584

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

1585

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

1586

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

1587

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

1588

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

1589

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

1590

1591

// Map coordinates to the reference square

1592

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

1593

x = -1.0;

1594

else

1595

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

1596

y = 2.0*y - 1.0;

1597

1598

// Reset values

1599

*values = 0;

1600

1601

// Map degree of freedom to element degree of freedom

1602

const unsigned int dof = i;

1603

1604

// Generate scalings

1605

const double scalings_y_0 = 1;

1606

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

1607

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

1608

1609

// Compute psitilde_a

1610

const double psitilde_a_0 = 1;

1611

const double psitilde_a_1 = x;

1612

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

1613

1614

// Compute psitilde_bs

1615

const double psitilde_bs_0_0 = 1;

1616

const double psitilde_bs_0_1 = 1.5*y + 0.5;

1617

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

1618

const double psitilde_bs_1_0 = 1;

1619

const double psitilde_bs_1_1 = 2.5*y + 1.5;

1620

const double psitilde_bs_2_0 = 1;

1621

1622

// Compute basisvalues

1623

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

1624

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

1625

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

1626

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

1627

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

1628

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

1629

1630

// Table(s) of coefficients

1631

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

1632

{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817},

1633

{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818},

1634

{0, 0, 0.2, 0, 0, 0.163299316185545},

1635

{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},

1636

{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},

1637

{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};

1638

1639

// Extract relevant coefficients

1640

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

1641

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

1642

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

1643

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

1644

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

1645

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

1646

1647

// Compute value(s)

1648

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;

1649

}

1650

1651

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

1652

virtual void evaluate_basis_all(double* values,

1653

const double* coordinates,

1654

const ufc::cell& c) const

1655

{

1656

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

1657

}

1658

1659

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

1660

virtual void evaluate_basis_derivatives(unsigned int i,

1661

unsigned int n,

1662

double* values,

1663

const double* coordinates,

1664

const ufc::cell& c) const

1665

{

1666

// Extract vertex coordinates

1667

const double * const * element_coordinates = c.coordinates;

1668

1669

// Compute Jacobian of affine map from reference cell

1670

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

1671

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

1672

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

1673

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

1674

1675

// Compute determinant of Jacobian

1676

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

1677

1678

// Compute inverse of Jacobian

1679

1680

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

1681

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

1682

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

1683

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

1684

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

1685

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

1686

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

1687

1688

// Map coordinates to the reference square

1689

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

1690

x = -1.0;

1691

else

1692

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

1693

y = 2.0*y - 1.0;

1694

1695

// Compute number of derivatives

1696

unsigned int num_derivatives = 1;

1697

1698

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

1699

num_derivatives *= 2;

1700

1701

1702

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

1703

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

1704

1705

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

1706

{

1707

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

1708

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

1709

combinations[j][k] = 0;

1710

}

1711

1712

// Generate combinations of derivatives

1713

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

1714

{

1715

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

1716

{

1717

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

1718

{

1719

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

1720

combinations[row][col] = 0;

1721

else

1722

{

1723

combinations[row][col] += 1;

1724

break;

1725

}

1726

}

1727

}

1728

}

1729

1730

// Compute inverse of Jacobian

1731

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

1732

1733

// Declare transformation matrix

1734

// Declare pointer to two dimensional array and initialise

1735

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

1736

1737

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

1738

{

1739

transform[j] = new double [num_derivatives];

1740

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

1741

transform[j][k] = 1;

1742

}

1743

1744

// Construct transformation matrix

1745

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

1746

{

1747

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

1748

{

1749

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

1750

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

1751

}

1752

}

1753

1754

// Reset values

1755

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

1756

values[j] = 0;

1757

1758

// Map degree of freedom to element degree of freedom

1759

const unsigned int dof = i;

1760

1761

// Generate scalings

1762

const double scalings_y_0 = 1;

1763

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

1764

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

1765

1766

// Compute psitilde_a

1767

const double psitilde_a_0 = 1;

1768

const double psitilde_a_1 = x;

1769

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

1770

1771

// Compute psitilde_bs

1772

const double psitilde_bs_0_0 = 1;

1773

const double psitilde_bs_0_1 = 1.5*y + 0.5;

1774

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

1775

const double psitilde_bs_1_0 = 1;

1776

const double psitilde_bs_1_1 = 2.5*y + 1.5;

1777

const double psitilde_bs_2_0 = 1;

1778

1779

// Compute basisvalues

1780

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

1781

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

1782

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

1783

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

1784

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

1785

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

1786

1787

// Table(s) of coefficients

1788

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

1789

{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817},

1790

{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818},

1791

{0, 0, 0.2, 0, 0, 0.163299316185545},

1792

{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},

1793

{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},

1794

{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};

1795

1796

// Interesting (new) part

1797

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

1798

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

1799

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

1800

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

1801

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

1802

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

1803

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

1804

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

1805

1806

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

1807

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

1808

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

1809

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

1810

{2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},

1811

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

1812

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

1813

1814

// Compute reference derivatives

1815

// Declare pointer to array of derivatives on FIAT element

1816

double *derivatives = new double [num_derivatives];

1817

1818

// Declare coefficients

1819

double coeff0_0 = 0;

1820

double coeff0_1 = 0;

1821

double coeff0_2 = 0;

1822

double coeff0_3 = 0;

1823

double coeff0_4 = 0;

1824

double coeff0_5 = 0;

1825

1826

// Declare new coefficients

1827

double new_coeff0_0 = 0;

1828

double new_coeff0_1 = 0;

1829

double new_coeff0_2 = 0;

1830

double new_coeff0_3 = 0;

1831

double new_coeff0_4 = 0;

1832

double new_coeff0_5 = 0;

1833

1834

// Loop possible derivatives

1835

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

1836

{

1837

// Get values from coefficients array

1838

new_coeff0_0 = coefficients0[dof][0];

1839

new_coeff0_1 = coefficients0[dof][1];

1840

new_coeff0_2 = coefficients0[dof][2];

1841

new_coeff0_3 = coefficients0[dof][3];

1842

new_coeff0_4 = coefficients0[dof][4];

1843

new_coeff0_5 = coefficients0[dof][5];

1844

1845

// Loop derivative order

1846

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

1847

{

1848

// Update old coefficients

1849

coeff0_0 = new_coeff0_0;

1850

coeff0_1 = new_coeff0_1;

1851

coeff0_2 = new_coeff0_2;

1852

coeff0_3 = new_coeff0_3;

1853

coeff0_4 = new_coeff0_4;

1854

coeff0_5 = new_coeff0_5;

1855

1856

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

1857

{

1858

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];

1859

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];

1860

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];

1861

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];

1862

new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];

1863

new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];

1864

}

1865

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

1866

{

1867

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];

1868

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];

1869

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];

1870

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];

1871

new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];

1872

new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];

1873

}

1874

1875

}

1876

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

1877

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;

1878

}

1879

1880

// Transform derivatives back to physical element

1881

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

1882

{

1883

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

1884

{

1885

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

1886

}

1887

}

1888

// Delete pointer to array of derivatives on FIAT element

1889

delete [] derivatives;

1890

1891

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

1892

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

1893

{

1894

delete [] combinations[row];

1895

delete [] transform[row];

1896

}

1897

1898

delete [] combinations;

1899

delete [] transform;

1900

}

1901

1902

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

1903

virtual void evaluate_basis_derivatives_all(unsigned int n,

1904

double* values,

1905

const double* coordinates,

1906

const ufc::cell& c) const

1907

{

1908

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

1909

}

1910

1911

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

1912

virtual double evaluate_dof(unsigned int i,

1913

const ufc::function& f,

1914

const ufc::cell& c) const

1915

{

1916

// The reference points, direction and weights:

1917

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

1918

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

1919

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

1920

1921

const double * const * x = c.coordinates;

1922

double result = 0.0;

1923

// Iterate over the points:

1924

// Evaluate basis functions for affine mapping

1925

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

1926

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

1927

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

1928

1929

// Compute affine mapping y = F(X)

1930

double y[2];

1931

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

1932

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

1933

1934

// Evaluate function at physical points

1935

double values[1];

1936

f.evaluate(values, y, c);

1937

1938

// Map function values using appropriate mapping

1939

// Affine map: Do nothing

1940

1941

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

1942

1943

// Take directional components

1944

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

1945

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

1946

// Multiply by weights

1947

result *= W[i][0];

1948

1949

return result;

1950

}

1951

1952

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

1953

virtual void evaluate_dofs(double* values,

1954

const ufc::function& f,

1955

const ufc::cell& c) const

1956

{

1957

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

1958

}

1959

1960

/// Interpolate vertex values from dof values

1961

virtual void interpolate_vertex_values(double* vertex_values,

1962

const double* dof_values,

1963

const ufc::cell& c) const

1964

{

1965

// Evaluate at vertices and use affine mapping

1966

vertex_values[0] = dof_values[0];

1967

vertex_values[1] = dof_values[1];

1968

vertex_values[2] = dof_values[2];

1969

}

1970

1971

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

1972

virtual unsigned int num_sub_elements() const

1973

{

1974

return 1;

1975

}

1976

1977

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

1978

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

1979

{

1980

return new UFC_PoissonP2LinearForm_finite_element_0();

1981

}

1982

1983

};

1984

1985

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

1986

1987

class UFC_PoissonP2LinearForm_finite_element_1: public ufc::finite_element

1988

{

1989

public:

1990

1991

/// Constructor

1992

UFC_PoissonP2LinearForm_finite_element_1() : ufc::finite_element()

1993

{

1994

// Do nothing

1995

}

1996

1997

/// Destructor

1998

virtual ~UFC_PoissonP2LinearForm_finite_element_1()

1999

{

2000

// Do nothing

2001

}

2002

2003

/// Return a string identifying the finite element

2004

virtual const char* signature() const

2005

{

2006

return "Lagrange finite element of degree 2 on a triangle";

2007

}

2008

2009

/// Return the cell shape

2010

virtual ufc::shape cell_shape() const

2011

{

2012

return ufc::triangle;

2013

}

2014

2015

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

2016

virtual unsigned int space_dimension() const

2017

{

2018

return 6;

2019

}

2020

2021

/// Return the rank of the value space

2022

virtual unsigned int value_rank() const

2023

{

2024

return 0;

2025

}

2026

2027

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

2028

virtual unsigned int value_dimension(unsigned int i) const

2029

{

2030

return 1;

2031

}

2032

2033

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

2034

virtual void evaluate_basis(unsigned int i,

2035

double* values,

2036

const double* coordinates,

2037

const ufc::cell& c) const

2038

{

2039

// Extract vertex coordinates

2040

const double * const * element_coordinates = c.coordinates;

2041

2042

// Compute Jacobian of affine map from reference cell

2043

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

2044

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

2045

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

2046

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

2047

2048

// Compute determinant of Jacobian

2049

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

2050

2051

// Compute inverse of Jacobian

2052

2053

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

2054

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

2055

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

2056

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

2057

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

2058

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

2059

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

2060

2061

// Map coordinates to the reference square

2062

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

2063

x = -1.0;

2064

else

2065

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

2066

y = 2.0*y - 1.0;

2067

2068

// Reset values

2069

*values = 0;

2070

2071

// Map degree of freedom to element degree of freedom

2072

const unsigned int dof = i;

2073

2074

// Generate scalings

2075

const double scalings_y_0 = 1;

2076

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

2077

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

2078

2079

// Compute psitilde_a

2080

const double psitilde_a_0 = 1;

2081

const double psitilde_a_1 = x;

2082

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

2083

2084

// Compute psitilde_bs

2085

const double psitilde_bs_0_0 = 1;

2086

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2087

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

2088

const double psitilde_bs_1_0 = 1;

2089

const double psitilde_bs_1_1 = 2.5*y + 1.5;

2090

const double psitilde_bs_2_0 = 1;

2091

2092

// Compute basisvalues

2093

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

2094

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

2095

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2096

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

2097

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

2098

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

2099

2100

// Table(s) of coefficients

2101

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

2102

{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817},

2103

{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818},

2104

{0, 0, 0.2, 0, 0, 0.163299316185545},

2105

{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},

2106

{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},

2107

{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};

2108

2109

// Extract relevant coefficients

2110

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

2111

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

2112

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

2113

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

2114

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

2115

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

2116

2117

// Compute value(s)

2118

*values = coeff0_0*basisvalue0 + coeff0_1*basisvalue1 + coeff0_2*basisvalue2 + coeff0_3*basisvalue3 + coeff0_4*basisvalue4 + coeff0_5*basisvalue5;

2119

}

2120

2121

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

2122

virtual void evaluate_basis_all(double* values,

2123

const double* coordinates,

2124

const ufc::cell& c) const

2125

{

2126

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

2127

}

2128

2129

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

2130

virtual void evaluate_basis_derivatives(unsigned int i,

2131

unsigned int n,

2132

double* values,

2133

const double* coordinates,

2134

const ufc::cell& c) const

2135

{

2136

// Extract vertex coordinates

2137

const double * const * element_coordinates = c.coordinates;

2138

2139

// Compute Jacobian of affine map from reference cell

2140

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

2141

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

2142

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

2143

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

2144

2145

// Compute determinant of Jacobian

2146

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

2147

2148

// Compute inverse of Jacobian

2149

2150

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

2151

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

2152

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

2153

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

2154

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

2155

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

2156

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

2157

2158

// Map coordinates to the reference square

2159

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

2160

x = -1.0;

2161

else

2162

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

2163

y = 2.0*y - 1.0;

2164

2165

// Compute number of derivatives

2166

unsigned int num_derivatives = 1;

2167

2168

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

2169

num_derivatives *= 2;

2170

2171

2172

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

2173

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

2174

2175

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

2176

{

2177

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

2178

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

2179

combinations[j][k] = 0;

2180

}

2181

2182

// Generate combinations of derivatives

2183

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

2184

{

2185

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

2186

{

2187

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

2188

{

2189

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

2190

combinations[row][col] = 0;

2191

else

2192

{

2193

combinations[row][col] += 1;

2194

break;

2195

}

2196

}

2197

}

2198

}

2199

2200

// Compute inverse of Jacobian

2201

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

2202

2203

// Declare transformation matrix

2204

// Declare pointer to two dimensional array and initialise

2205

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

2206

2207

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

2208

{

2209

transform[j] = new double [num_derivatives];

2210

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

2211

transform[j][k] = 1;

2212

}

2213

2214

// Construct transformation matrix

2215

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

2216

{

2217

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

2218

{

2219

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

2220

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

2221

}

2222

}

2223

2224

// Reset values

2225

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

2226

values[j] = 0;

2227

2228

// Map degree of freedom to element degree of freedom

2229

const unsigned int dof = i;

2230

2231

// Generate scalings

2232

const double scalings_y_0 = 1;

2233

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

2234

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

2235

2236

// Compute psitilde_a

2237

const double psitilde_a_0 = 1;

2238

const double psitilde_a_1 = x;

2239

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

2240

2241

// Compute psitilde_bs

2242

const double psitilde_bs_0_0 = 1;

2243

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2244

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

2245

const double psitilde_bs_1_0 = 1;

2246

const double psitilde_bs_1_1 = 2.5*y + 1.5;

2247

const double psitilde_bs_2_0 = 1;

2248

2249

// Compute basisvalues

2250

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

2251

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

2252

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2253

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

2254

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

2255

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

2256

2257

// Table(s) of coefficients

2258

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

2259

{{0, -0.173205080756888, -0.1, 0.121716123890037, 0.0942809041582063, 0.0544331053951817},

2260

{0, 0.173205080756888, -0.1, 0.121716123890037, -0.0942809041582064, 0.0544331053951818},

2261

{0, 0, 0.2, 0, 0, 0.163299316185545},

2262

{0.471404520791032, 0.23094010767585, 0.133333333333333, 0, 0.188561808316413, -0.163299316185545},

2263

{0.471404520791032, -0.23094010767585, 0.133333333333333, 0, -0.188561808316413, -0.163299316185545},

2264

{0.471404520791032, 0, -0.266666666666667, -0.243432247780074, 0, 0.0544331053951817}};

2265

2266

// Interesting (new) part

2267

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

2268

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

2269

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

2270

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

2271

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

2272

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

2273

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

2274

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

2275

2276

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

2277

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

2278

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

2279

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

2280

{2.58198889747161, 4.74341649025257, -0.912870929175277, 0, 0, 0},

2281

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

2282

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

2283

2284

// Compute reference derivatives

2285

// Declare pointer to array of derivatives on FIAT element

2286

double *derivatives = new double [num_derivatives];

2287

2288

// Declare coefficients

2289

double coeff0_0 = 0;

2290

double coeff0_1 = 0;

2291

double coeff0_2 = 0;

2292

double coeff0_3 = 0;

2293

double coeff0_4 = 0;

2294

double coeff0_5 = 0;

2295

2296

// Declare new coefficients

2297

double new_coeff0_0 = 0;

2298

double new_coeff0_1 = 0;

2299

double new_coeff0_2 = 0;

2300

double new_coeff0_3 = 0;

2301

double new_coeff0_4 = 0;

2302

double new_coeff0_5 = 0;

2303

2304

// Loop possible derivatives

2305

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

2306

{

2307

// Get values from coefficients array

2308

new_coeff0_0 = coefficients0[dof][0];

2309

new_coeff0_1 = coefficients0[dof][1];

2310

new_coeff0_2 = coefficients0[dof][2];

2311

new_coeff0_3 = coefficients0[dof][3];

2312

new_coeff0_4 = coefficients0[dof][4];

2313

new_coeff0_5 = coefficients0[dof][5];

2314

2315

// Loop derivative order

2316

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

2317

{

2318

// Update old coefficients

2319

coeff0_0 = new_coeff0_0;

2320

coeff0_1 = new_coeff0_1;

2321

coeff0_2 = new_coeff0_2;

2322

coeff0_3 = new_coeff0_3;

2323

coeff0_4 = new_coeff0_4;

2324

coeff0_5 = new_coeff0_5;

2325

2326

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

2327

{

2328

new_coeff0_0 = coeff0_0*dmats0[0][0] + coeff0_1*dmats0[1][0] + coeff0_2*dmats0[2][0] + coeff0_3*dmats0[3][0] + coeff0_4*dmats0[4][0] + coeff0_5*dmats0[5][0];

2329

new_coeff0_1 = coeff0_0*dmats0[0][1] + coeff0_1*dmats0[1][1] + coeff0_2*dmats0[2][1] + coeff0_3*dmats0[3][1] + coeff0_4*dmats0[4][1] + coeff0_5*dmats0[5][1];

2330

new_coeff0_2 = coeff0_0*dmats0[0][2] + coeff0_1*dmats0[1][2] + coeff0_2*dmats0[2][2] + coeff0_3*dmats0[3][2] + coeff0_4*dmats0[4][2] + coeff0_5*dmats0[5][2];

2331

new_coeff0_3 = coeff0_0*dmats0[0][3] + coeff0_1*dmats0[1][3] + coeff0_2*dmats0[2][3] + coeff0_3*dmats0[3][3] + coeff0_4*dmats0[4][3] + coeff0_5*dmats0[5][3];

2332

new_coeff0_4 = coeff0_0*dmats0[0][4] + coeff0_1*dmats0[1][4] + coeff0_2*dmats0[2][4] + coeff0_3*dmats0[3][4] + coeff0_4*dmats0[4][4] + coeff0_5*dmats0[5][4];

2333

new_coeff0_5 = coeff0_0*dmats0[0][5] + coeff0_1*dmats0[1][5] + coeff0_2*dmats0[2][5] + coeff0_3*dmats0[3][5] + coeff0_4*dmats0[4][5] + coeff0_5*dmats0[5][5];

2334

}

2335

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

2336

{

2337

new_coeff0_0 = coeff0_0*dmats1[0][0] + coeff0_1*dmats1[1][0] + coeff0_2*dmats1[2][0] + coeff0_3*dmats1[3][0] + coeff0_4*dmats1[4][0] + coeff0_5*dmats1[5][0];

2338

new_coeff0_1 = coeff0_0*dmats1[0][1] + coeff0_1*dmats1[1][1] + coeff0_2*dmats1[2][1] + coeff0_3*dmats1[3][1] + coeff0_4*dmats1[4][1] + coeff0_5*dmats1[5][1];

2339

new_coeff0_2 = coeff0_0*dmats1[0][2] + coeff0_1*dmats1[1][2] + coeff0_2*dmats1[2][2] + coeff0_3*dmats1[3][2] + coeff0_4*dmats1[4][2] + coeff0_5*dmats1[5][2];

2340

new_coeff0_3 = coeff0_0*dmats1[0][3] + coeff0_1*dmats1[1][3] + coeff0_2*dmats1[2][3] + coeff0_3*dmats1[3][3] + coeff0_4*dmats1[4][3] + coeff0_5*dmats1[5][3];

2341

new_coeff0_4 = coeff0_0*dmats1[0][4] + coeff0_1*dmats1[1][4] + coeff0_2*dmats1[2][4] + coeff0_3*dmats1[3][4] + coeff0_4*dmats1[4][4] + coeff0_5*dmats1[5][4];

2342

new_coeff0_5 = coeff0_0*dmats1[0][5] + coeff0_1*dmats1[1][5] + coeff0_2*dmats1[2][5] + coeff0_3*dmats1[3][5] + coeff0_4*dmats1[4][5] + coeff0_5*dmats1[5][5];

2343

}

2344

2345

}

2346

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

2347

derivatives[deriv_num] = new_coeff0_0*basisvalue0 + new_coeff0_1*basisvalue1 + new_coeff0_2*basisvalue2 + new_coeff0_3*basisvalue3 + new_coeff0_4*basisvalue4 + new_coeff0_5*basisvalue5;

2348

}

2349

2350

// Transform derivatives back to physical element

2351

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

2352

{

2353

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

2354

{

2355

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

2356

}

2357

}

2358

// Delete pointer to array of derivatives on FIAT element

2359

delete [] derivatives;

2360

2361

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

2362

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

2363

{

2364

delete [] combinations[row];

2365

delete [] transform[row];

2366

}

2367

2368

delete [] combinations;

2369

delete [] transform;

2370

}

2371

2372

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

2373

virtual void evaluate_basis_derivatives_all(unsigned int n,

2374

double* values,

2375

const double* coordinates,

2376

const ufc::cell& c) const

2377

{

2378

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

2379

}

2380

2381

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

2382

virtual double evaluate_dof(unsigned int i,

2383

const ufc::function& f,

2384

const ufc::cell& c) const

2385

{

2386

// The reference points, direction and weights:

2387

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

2388

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

2389

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

2390

2391

const double * const * x = c.coordinates;

2392

double result = 0.0;

2393

// Iterate over the points:

2394

// Evaluate basis functions for affine mapping

2395

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

2396

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

2397

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

2398

2399

// Compute affine mapping y = F(X)

2400

double y[2];

2401

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

2402

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

2403

2404

// Evaluate function at physical points

2405

double values[1];

2406

f.evaluate(values, y, c);

2407

2408

// Map function values using appropriate mapping

2409

// Affine map: Do nothing

2410

2411

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

2412

2413

// Take directional components

2414

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

2415

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

2416

// Multiply by weights

2417

result *= W[i][0];

2418

2419

return result;

2420

}

2421

2422

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

2423

virtual void evaluate_dofs(double* values,

2424

const ufc::function& f,

2425

const ufc::cell& c) const

2426

{

2427

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

2428

}

2429

2430

/// Interpolate vertex values from dof values

2431

virtual void interpolate_vertex_values(double* vertex_values,

2432

const double* dof_values,

2433

const ufc::cell& c) const

2434

{

2435

// Evaluate at vertices and use affine mapping

2436

vertex_values[0] = dof_values[0];

2437

vertex_values[1] = dof_values[1];

2438

vertex_values[2] = dof_values[2];

2439

}

2440

2441

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

2442

virtual unsigned int num_sub_elements() const

2443

{

2444

return 1;

2445

}

2446

2447

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

2448

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

2449

{

2450

return new UFC_PoissonP2LinearForm_finite_element_1();

2451

}

2452

2453

};

2454

2455

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

2456

/// degrees of freedom (dofs).

2457

2458

class UFC_PoissonP2LinearForm_dof_map_0: public ufc::dof_map

2459

{

2460

private:

2461

2462

unsigned int __global_dimension;

2463

2464

public:

2465

2466

/// Constructor

2467

UFC_PoissonP2LinearForm_dof_map_0() : ufc::dof_map()

2468

{

2469

__global_dimension = 0;

2470

}

2471

2472

/// Destructor

2473

virtual ~UFC_PoissonP2LinearForm_dof_map_0()

2474

{

2475

// Do nothing

2476

}

2477

2478

/// Return a string identifying the dof map

2479

virtual const char* signature() const

2480

{

2481

return "FFC dof map for Lagrange finite element of degree 2 on a triangle";

2482

}

2483

2484

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

2485

virtual bool needs_mesh_entities(unsigned int d) const

2486

{

2487

switch ( d )

2488

{

2489

case 0:

2490

return true;

2491

break;

2492

case 1:

2493

return true;

2494

break;

2495

case 2:

2496

return false;

2497

break;

2498

}

2499

return false;

2500

}

2501

2502

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

2503

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

2504

{

2505

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

2506

return false;

2507

}

2508

2509

/// Initialize dof map for given cell

2510

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

2511

const ufc::cell& c)

2512

{

2513

// Do nothing

2514

}

2515

2516

/// Finish initialization of dof map for cells

2517

virtual void init_cell_finalize()

2518

{

2519

// Do nothing

2520

}

2521

2522

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

2523

virtual unsigned int global_dimension() const

2524

{

2525

return __global_dimension;

2526

}

2527

2528

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

2529

virtual unsigned int local_dimension() const

2530

{

2531

return 6;

2532

}

2533

2534

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

2535

virtual unsigned int geometric_dimension() const

2536

{

2537

return 2;

2538

}

2539

2540

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

2541

virtual unsigned int num_facet_dofs() const

2542

{

2543

return 3;

2544

}

2545

2546

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

2547

virtual unsigned int num_entity_dofs(unsigned int d) const

2548

{

2549

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

2550

}

2551

2552

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

2553

virtual void tabulate_dofs(unsigned int* dofs,

2554

const ufc::mesh& m,

2555

const ufc::cell& c) const

2556

{

2557

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

2558

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

2559

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

2560

unsigned int offset = m.num_entities[0];

2561

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

2562

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

2563

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

2564

}

2565

2566

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

2567

virtual void tabulate_facet_dofs(unsigned int* dofs,

2568

unsigned int facet) const

2569

{

2570

switch ( facet )

2571

{

2572

case 0:

2573

dofs[0] = 1;

2574

dofs[1] = 2;

2575

dofs[2] = 3;

2576

break;

2577

case 1:

2578

dofs[0] = 0;

2579

dofs[1] = 2;

2580

dofs[2] = 4;

2581

break;

2582

case 2:

2583

dofs[0] = 0;

2584

dofs[1] = 1;

2585

dofs[2] = 5;

2586

break;

2587

}

2588

}

2589

2590

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

2591

virtual void tabulate_entity_dofs(unsigned int* dofs,

2592

unsigned int d, unsigned int i) const

2593

{

2594

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

2595

}

2596

2597

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

2598

virtual void tabulate_coordinates(double** coordinates,

2599

const ufc::cell& c) const

2600

{

2601

const double * const * x = c.coordinates;

2602

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

2603

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

2604

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

2605

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

2606

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

2607

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

2608

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

2609

coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];

2610

coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];

2611

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

2612

coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];

2613

coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];

2614

}

2615

2616

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

2617

virtual unsigned int num_sub_dof_maps() const

2618

{

2619

return 1;

2620

}

2621

2622

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

2623

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

2624

{

2625

return new UFC_PoissonP2LinearForm_dof_map_0();

2626

}

2627

2628

};

2629

2630

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

2631

/// degrees of freedom (dofs).

2632

2633

class UFC_PoissonP2LinearForm_dof_map_1: public ufc::dof_map

2634

{

2635

private:

2636

2637

unsigned int __global_dimension;

2638

2639

public:

2640

2641

/// Constructor

2642

UFC_PoissonP2LinearForm_dof_map_1() : ufc::dof_map()

2643

{

2644

__global_dimension = 0;

2645

}

2646

2647

/// Destructor

2648

virtual ~UFC_PoissonP2LinearForm_dof_map_1()

2649

{

2650

// Do nothing

2651

}

2652

2653

/// Return a string identifying the dof map

2654

virtual const char* signature() const

2655

{

2656

return "FFC dof map for Lagrange finite element of degree 2 on a triangle";

2657

}

2658

2659

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

2660

virtual bool needs_mesh_entities(unsigned int d) const

2661

{

2662

switch ( d )

2663

{

2664

case 0:

2665

return true;

2666

break;

2667

case 1:

2668

return true;

2669

break;

2670

case 2:

2671

return false;

2672

break;

2673

}

2674

return false;

2675

}

2676

2677

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

2678

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

2679

{

2680

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

2681

return false;

2682

}

2683

2684

/// Initialize dof map for given cell

2685

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

2686

const ufc::cell& c)

2687

{

2688

// Do nothing

2689

}

2690

2691

/// Finish initialization of dof map for cells

2692

virtual void init_cell_finalize()

2693

{

2694

// Do nothing

2695

}

2696

2697

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

2698

virtual unsigned int global_dimension() const

2699

{

2700

return __global_dimension;

2701

}

2702

2703

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

2704

virtual unsigned int local_dimension() const

2705

{

2706

return 6;

2707

}

2708

2709

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

2710

virtual unsigned int geometric_dimension() const

2711

{

2712

return 2;

2713

}

2714

2715

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

2716

virtual unsigned int num_facet_dofs() const

2717

{

2718

return 3;

2719

}

2720

2721

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

2722

virtual unsigned int num_entity_dofs(unsigned int d) const

2723

{

2724

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

2725

}

2726

2727

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

2728

virtual void tabulate_dofs(unsigned int* dofs,

2729

const ufc::mesh& m,

2730

const ufc::cell& c) const

2731

{

2732

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

2733

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

2734

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

2735

unsigned int offset = m.num_entities[0];

2736

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

2737

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

2738

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

2739

}

2740

2741

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

2742

virtual void tabulate_facet_dofs(unsigned int* dofs,

2743

unsigned int facet) const

2744

{

2745

switch ( facet )

2746

{

2747

case 0:

2748

dofs[0] = 1;

2749

dofs[1] = 2;

2750

dofs[2] = 3;

2751

break;

2752

case 1:

2753

dofs[0] = 0;

2754

dofs[1] = 2;

2755

dofs[2] = 4;

2756

break;

2757

case 2:

2758

dofs[0] = 0;

2759

dofs[1] = 1;

2760

dofs[2] = 5;

2761

break;

2762

}

2763

}

2764

2765

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

2766

virtual void tabulate_entity_dofs(unsigned int* dofs,

2767

unsigned int d, unsigned int i) const

2768

{

2769

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

2770

}

2771

2772

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

2773

virtual void tabulate_coordinates(double** coordinates,

2774

const ufc::cell& c) const

2775

{

2776

const double * const * x = c.coordinates;

2777

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

2778

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

2779

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

2780

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

2781

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

2782

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

2783

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

2784

coordinates[3][1] = 0.5*x[1][1] + 0.5*x[2][1];

2785

coordinates[4][0] = 0.5*x[0][0] + 0.5*x[2][0];

2786

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

2787

coordinates[5][0] = 0.5*x[0][0] + 0.5*x[1][0];

2788

coordinates[5][1] = 0.5*x[0][1] + 0.5*x[1][1];

2789

}

2790

2791

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

2792

virtual unsigned int num_sub_dof_maps() const

2793

{

2794

return 1;

2795

}

2796

2797

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

2798

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

2799

{

2800

return new UFC_PoissonP2LinearForm_dof_map_1();

2801

}

2802

2803

};

2804

2805

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

2806

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

2807

/// the integral over a cell.

2808

2809

class UFC_PoissonP2LinearForm_cell_integral_0: public ufc::cell_integral

2810

{

2811

public:

2812

2813

/// Constructor

2814

UFC_PoissonP2LinearForm_cell_integral_0() : ufc::cell_integral()

2815

{

2816

// Do nothing

2817

}

2818

2819

/// Destructor

2820

virtual ~UFC_PoissonP2LinearForm_cell_integral_0()

2821

{

2822

// Do nothing

2823

}

2824

2825

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

2826

virtual void tabulate_tensor(double* A,

2827

const double * const * w,

2828

const ufc::cell& c) const

2829

{

2830

// Extract vertex coordinates

2831

const double * const * x = c.coordinates;

2832

2833

// Compute Jacobian of affine map from reference cell

2834

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

2835

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

2836

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

2837

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

2838

2839

// Compute determinant of Jacobian

2840

double detJ = J_00*J_11 - J_01*J_10;

2841

2842

// Compute inverse of Jacobian

2843

2844

// Set scale factor

2845

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

2846

2847

// Compute coefficients

2848

const double c0_0_0_0 = w[0][0];

2849

const double c0_0_0_1 = w[0][1];

2850

const double c0_0_0_2 = w[0][2];

2851

const double c0_0_0_3 = w[0][3];

2852

const double c0_0_0_4 = w[0][4];

2853

const double c0_0_0_5 = w[0][5];

2854

2855

// Compute geometry tensors

2856

const double G0_0 = det*c0_0_0_0;

2857

const double G0_1 = det*c0_0_0_1;

2858

const double G0_2 = det*c0_0_0_2;

2859

const double G0_3 = det*c0_0_0_3;

2860

const double G0_4 = det*c0_0_0_4;

2861

const double G0_5 = det*c0_0_0_5;

2862

2863

// Compute element tensor

2864

A[0] = 0.0166666666666666*G0_0 - 0.00277777777777777*G0_1 - 0.00277777777777778*G0_2 - 0.0111111111111111*G0_3;

2865

A[1] = -0.00277777777777777*G0_0 + 0.0166666666666666*G0_1 - 0.00277777777777777*G0_2 - 0.0111111111111111*G0_4;

2866

A[2] = -0.00277777777777778*G0_0 - 0.00277777777777777*G0_1 + 0.0166666666666666*G0_2 - 0.0111111111111111*G0_5;

2867

A[3] = -0.0111111111111111*G0_0 + 0.0888888888888888*G0_3 + 0.0444444444444444*G0_4 + 0.0444444444444444*G0_5;

2868

A[4] = -0.0111111111111111*G0_1 + 0.0444444444444444*G0_3 + 0.0888888888888888*G0_4 + 0.0444444444444444*G0_5;

2869

A[5] = -0.0111111111111111*G0_2 + 0.0444444444444444*G0_3 + 0.0444444444444444*G0_4 + 0.0888888888888888*G0_5;

2870

}

2871

2872

};

2873

2874

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

2875

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

2876

/// mapping

2877

///

2878

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

2879

///

2880

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

2881

/// global tensor A is defined by

2882

///

2883

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

2884

///

2885

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

2886

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

2887

/// fixed functions (coefficients).

2888

2889

class UFC_PoissonP2LinearForm: public ufc::form

2890

{

2891

public:

2892

2893

/// Constructor

2894

UFC_PoissonP2LinearForm() : ufc::form()

2895

{

2896

// Do nothing

2897

}

2898

2899

/// Destructor

2900

virtual ~UFC_PoissonP2LinearForm()

2901

{

2902

// Do nothing

2903

}

2904

2905

/// Return a string identifying the form

2906

virtual const char* signature() const

2907

{

2908

return "w0_a0[0, 1, 2, 3, 4, 5] | vi0[0, 1, 2, 3, 4, 5]*va0[0, 1, 2, 3, 4, 5]*dX(0)";

2909

}

2910

2911

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

2912

virtual unsigned int rank() const

2913

{

2914

return 1;

2915

}

2916

2917

/// Return the number of coefficients (n)

2918

virtual unsigned int num_coefficients() const

2919

{

2920

return 1;

2921

}

2922

2923

/// Return the number of cell integrals

2924

virtual unsigned int num_cell_integrals() const

2925

{

2926

return 1;

2927

}

2928

2929

/// Return the number of exterior facet integrals

2930

virtual unsigned int num_exterior_facet_integrals() const

2931

{

2932

return 0;

2933

}

2934

2935

/// Return the number of interior facet integrals

2936

virtual unsigned int num_interior_facet_integrals() const

2937

{

2938

return 0;

2939

}

2940

2941

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

2942

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

2943

{

2944

switch ( i )

2945

{

2946

case 0:

2947

return new UFC_PoissonP2LinearForm_finite_element_0();

2948

break;

2949

case 1:

2950

return new UFC_PoissonP2LinearForm_finite_element_1();

2951

break;

2952

}

2953

return 0;

2954

}

2955

2956

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

2957

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

2958

{

2959

switch ( i )

2960

{

2961

case 0:

2962

return new UFC_PoissonP2LinearForm_dof_map_0();

2963

break;

2964

case 1:

2965

return new UFC_PoissonP2LinearForm_dof_map_1();

2966

break;

2967

}

2968

return 0;

2969

}

2970

2971

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

2972

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

2973

{

2974

return new UFC_PoissonP2LinearForm_cell_integral_0();

2975

}

2976

2977

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

2978

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

2979

{

2980

return 0;

2981

}

2982

2983

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

2984

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

2985

{

2986

return 0;

2987

}

2988

2989

};

2990

2991

// DOLFIN wrappers

2992

2993

#include <dolfin/fem/Form.h>

2994

2995

class PoissonP2BilinearForm : public dolfin::Form

2996

{

2997

public:

2998

2999

PoissonP2BilinearForm() : dolfin::Form()

3000

{

3001

// Do nothing

3002

}

3003

3004

/// Return UFC form

3005

virtual const ufc::form& form() const

3006

{

3007

return __form;

3008

}

3009

3010

/// Return array of coefficients

3011

virtual const dolfin::Array<dolfin::Function*>& coefficients() const

3012

{

3013

return __coefficients;

3014

}

3015

3016

private:

3017

3018

// UFC form

3019

UFC_PoissonP2BilinearForm __form;

3020

3021

/// Array of coefficients

3022

dolfin::Array<dolfin::Function*> __coefficients;

3023

3024

};

3025

3026

class PoissonP2LinearForm : public dolfin::Form

3027

{

3028

public:

3029

3030

PoissonP2LinearForm(dolfin::Function& w0) : dolfin::Form()

3031

{

3032

__coefficients.push_back(&w0);

3033

}

3034

3035

/// Return UFC form

3036

virtual const ufc::form& form() const

3037

{

3038

return __form;

3039

}

3040

3041

/// Return array of coefficients

3042

virtual const dolfin::Array<dolfin::Function*>& coefficients() const

3043

{

3044

return __coefficients;

3045

}

3046

3047

private:

3048

3049

// UFC form

3050

UFC_PoissonP2LinearForm __form;

3051

3052

/// Array of coefficients

3053

dolfin::Array<dolfin::Function*> __coefficients;

3054

3055

};

3056

3057

#endif

Older »