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

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_Poisson2D_2BilinearForm_finite_element_0();

481

}

482

483

};

484

485

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

486

487

class UFC_Poisson2D_2BilinearForm_finite_element_1: public ufc::finite_element

488

{

489

public:

490

491

/// Constructor

492

UFC_Poisson2D_2BilinearForm_finite_element_1() : ufc::finite_element()

493

{

494

// Do nothing

495

}

496

497

/// Destructor

498

virtual ~UFC_Poisson2D_2BilinearForm_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 "FiniteElement('Lagrange', 'triangle', 2)";

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.89897948556636, 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_Poisson2D_2BilinearForm_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_Poisson2D_2BilinearForm_dof_map_0: public ufc::dof_map

959

{

960

private:

961

962

unsigned int __global_dimension;

963

964

public:

965

966

/// Constructor

967

UFC_Poisson2D_2BilinearForm_dof_map_0() : ufc::dof_map()

968

{

969

__global_dimension = 0;

970

}

971

972

/// Destructor

973

virtual ~UFC_Poisson2D_2BilinearForm_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 FiniteElement('Lagrange', 'triangle', 2)";

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_Poisson2D_2BilinearForm_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_Poisson2D_2BilinearForm_dof_map_1: public ufc::dof_map

1134

{

1135

private:

1136

1137

unsigned int __global_dimension;

1138

1139

public:

1140

1141

/// Constructor

1142

UFC_Poisson2D_2BilinearForm_dof_map_1() : ufc::dof_map()

1143

{

1144

__global_dimension = 0;

1145

}

1146

1147

/// Destructor

1148

virtual ~UFC_Poisson2D_2BilinearForm_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 FiniteElement('Lagrange', 'triangle', 2)";

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_Poisson2D_2BilinearForm_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_Poisson2D_2BilinearForm_cell_integral_0: public ufc::cell_integral

1310

{

1311

public:

1312

1313

/// Constructor

1314

UFC_Poisson2D_2BilinearForm_cell_integral_0() : ufc::cell_integral()

1315

{

1316

// Do nothing

1317

}

1318

1319

/// Destructor

1320

virtual ~UFC_Poisson2D_2BilinearForm_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

// Number of operations to compute element tensor = 114

1352

// Compute geometry tensors

1353

// Number of operations to compute decalrations = 16

1354

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

1355

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

1356

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

1357

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

1358

1359

// Compute element tensor

1360

// Number of operations to compute tensor = 98

1361

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

1362

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

1363

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

1364

A[3] = 0;

1365

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

1366

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

1367

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

1368

A[7] = 0.499999999999999*G0_0_0;

1369

A[8] = -0.166666666666666*G0_0_1;

1370

A[9] = 0.666666666666665*G0_0_1;

1371

A[10] = 0;

1372

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

1373

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

1374

A[13] = -0.166666666666666*G0_1_0;

1375

A[14] = 0.499999999999999*G0_1_1;

1376

A[15] = 0.666666666666666*G0_1_0;

1377

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

1378

A[17] = 0;

1379

A[18] = 0;

1380

A[19] = 0.666666666666665*G0_1_0;

1381

A[20] = 0.666666666666666*G0_0_1;

1382

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

1383

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

1384

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

1385

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

1386

A[25] = 0;

1387

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

1388

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

1389

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

1390

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

1391

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

1392

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

1393

A[32] = 0;

1394

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

1395

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

1396

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

1397

}

1398

1399

};

1400

1401

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

1402

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

1403

/// mapping

1404

///

1405

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

1406

///

1407

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

1408

/// global tensor A is defined by

1409

///

1410

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

1411

///

1412

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

1413

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

1414

/// fixed functions (coefficients).

1415

1416

class UFC_Poisson2D_2BilinearForm: public ufc::form

1417

{

1418

public:

1419

1420

/// Constructor

1421

UFC_Poisson2D_2BilinearForm() : ufc::form()

1422

{

1423

// Do nothing

1424

}

1425

1426

/// Destructor

1427

virtual ~UFC_Poisson2D_2BilinearForm()

1428

{

1429

// Do nothing

1430

}

1431

1432

/// Return a string identifying the form

1433

virtual const char* signature() const

1434

{

1435

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)";

1436

}

1437

1438

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

1439

virtual unsigned int rank() const

1440

{

1441

return 2;

1442

}

1443

1444

/// Return the number of coefficients (n)

1445

virtual unsigned int num_coefficients() const

1446

{

1447

return 0;

1448

}

1449

1450

/// Return the number of cell integrals

1451

virtual unsigned int num_cell_integrals() const

1452

{

1453

return 1;

1454

}

1455

1456

/// Return the number of exterior facet integrals

1457

virtual unsigned int num_exterior_facet_integrals() const

1458

{

1459

return 0;

1460

}

1461

1462

/// Return the number of interior facet integrals

1463

virtual unsigned int num_interior_facet_integrals() const

1464

{

1465

return 0;

1466

}

1467

1468

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

1469

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

1470

{

1471

switch ( i )

1472

{

1473

case 0:

1474

return new UFC_Poisson2D_2BilinearForm_finite_element_0();

1475

break;

1476

case 1:

1477

return new UFC_Poisson2D_2BilinearForm_finite_element_1();

1478

break;

1479

}

1480

return 0;

1481

}

1482

1483

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

1484

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

1485

{

1486

switch ( i )

1487

{

1488

case 0:

1489

return new UFC_Poisson2D_2BilinearForm_dof_map_0();

1490

break;

1491

case 1:

1492

return new UFC_Poisson2D_2BilinearForm_dof_map_1();

1493

break;

1494

}

1495

return 0;

1496

}

1497

1498

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

1499

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

1500

{

1501

return new UFC_Poisson2D_2BilinearForm_cell_integral_0();

1502

}

1503

1504

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

1505

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

1506

{

1507

return 0;

1508

}

1509

1510

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

1511

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

1512

{

1513

return 0;

1514

}

1515

1516

};

1517

1518

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

1519

1520

class UFC_Poisson2D_2LinearForm_finite_element_0: public ufc::finite_element

1521

{

1522

public:

1523

1524

/// Constructor

1525

UFC_Poisson2D_2LinearForm_finite_element_0() : ufc::finite_element()

1526

{

1527

// Do nothing

1528

}

1529

1530

/// Destructor

1531

virtual ~UFC_Poisson2D_2LinearForm_finite_element_0()

1532

{

1533

// Do nothing

1534

}

1535

1536

/// Return a string identifying the finite element

1537

virtual const char* signature() const

1538

{

1539

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

1540

}

1541

1542

/// Return the cell shape

1543

virtual ufc::shape cell_shape() const

1544

{

1545

return ufc::triangle;

1546

}

1547

1548

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

1549

virtual unsigned int space_dimension() const

1550

{

1551

return 6;

1552

}

1553

1554

/// Return the rank of the value space

1555

virtual unsigned int value_rank() const

1556

{

1557

return 0;

1558

}

1559

1560

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

1561

virtual unsigned int value_dimension(unsigned int i) const

1562

{

1563

return 1;

1564

}

1565

1566

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

1567

virtual void evaluate_basis(unsigned int i,

1568

double* values,

1569

const double* coordinates,

1570

const ufc::cell& c) const

1571

{

1572

// Extract vertex coordinates

1573

const double * const * element_coordinates = c.coordinates;

1574

1575

// Compute Jacobian of affine map from reference cell

1576

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

1577

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

1578

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

1579

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

1580

1581

// Compute determinant of Jacobian

1582

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

1583

1584

// Compute inverse of Jacobian

1585

1586

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

1587

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

1588

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

1589

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

1590

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

1591

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

1592

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

1593

1594

// Map coordinates to the reference square

1595

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

1596

x = -1.0;

1597

else

1598

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

1599

y = 2.0*y - 1.0;

1600

1601

// Reset values

1602

*values = 0;

1603

1604

// Map degree of freedom to element degree of freedom

1605

const unsigned int dof = i;

1606

1607

// Generate scalings

1608

const double scalings_y_0 = 1;

1609

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

1610

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

1611

1612

// Compute psitilde_a

1613

const double psitilde_a_0 = 1;

1614

const double psitilde_a_1 = x;

1615

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

1616

1617

// Compute psitilde_bs

1618

const double psitilde_bs_0_0 = 1;

1619

const double psitilde_bs_0_1 = 1.5*y + 0.5;

1620

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;

1621

const double psitilde_bs_1_0 = 1;

1622

const double psitilde_bs_1_1 = 2.5*y + 1.5;

1623

const double psitilde_bs_2_0 = 1;

1624

1625

// Compute basisvalues

1626

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

1627

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

1628

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

1629

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

1630

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

1631

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

1632

1633

// Table(s) of coefficients

1634

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

1635

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

1636

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

1637

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

1638

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

1639

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

1640

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

1641

1642

// Extract relevant coefficients

1643

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

1644

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

1645

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

1646

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

1647

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

1648

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

1649

1650

// Compute value(s)

1651

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

1652

}

1653

1654

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

1655

virtual void evaluate_basis_all(double* values,

1656

const double* coordinates,

1657

const ufc::cell& c) const

1658

{

1659

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

1660

}

1661

1662

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

1663

virtual void evaluate_basis_derivatives(unsigned int i,

1664

unsigned int n,

1665

double* values,

1666

const double* coordinates,

1667

const ufc::cell& c) const

1668

{

1669

// Extract vertex coordinates

1670

const double * const * element_coordinates = c.coordinates;

1671

1672

// Compute Jacobian of affine map from reference cell

1673

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

1674

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

1675

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

1676

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

1677

1678

// Compute determinant of Jacobian

1679

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

1680

1681

// Compute inverse of Jacobian

1682

1683

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

1684

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

1685

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

1686

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

1687

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

1688

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

1689

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

1690

1691

// Map coordinates to the reference square

1692

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

1693

x = -1.0;

1694

else

1695

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

1696

y = 2.0*y - 1.0;

1697

1698

// Compute number of derivatives

1699

unsigned int num_derivatives = 1;

1700

1701

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

1702

num_derivatives *= 2;

1703

1704

1705

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

1706

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

1707

1708

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

1709

{

1710

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

1711

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

1712

combinations[j][k] = 0;

1713

}

1714

1715

// Generate combinations of derivatives

1716

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

1717

{

1718

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

1719

{

1720

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

1721

{

1722

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

1723

combinations[row][col] = 0;

1724

else

1725

{

1726

combinations[row][col] += 1;

1727

break;

1728

}

1729

}

1730

}

1731

}

1732

1733

// Compute inverse of Jacobian

1734

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

1735

1736

// Declare transformation matrix

1737

// Declare pointer to two dimensional array and initialise

1738

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

1739

1740

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

1741

{

1742

transform[j] = new double [num_derivatives];

1743

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

1744

transform[j][k] = 1;

1745

}

1746

1747

// Construct transformation matrix

1748

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

1749

{

1750

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

1751

{

1752

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

1753

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

1754

}

1755

}

1756

1757

// Reset values

1758

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

1759

values[j] = 0;

1760

1761

// Map degree of freedom to element degree of freedom

1762

const unsigned int dof = i;

1763

1764

// Generate scalings

1765

const double scalings_y_0 = 1;

1766

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

1767

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

1768

1769

// Compute psitilde_a

1770

const double psitilde_a_0 = 1;

1771

const double psitilde_a_1 = x;

1772

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

1773

1774

// Compute psitilde_bs

1775

const double psitilde_bs_0_0 = 1;

1776

const double psitilde_bs_0_1 = 1.5*y + 0.5;

1777

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;

1778

const double psitilde_bs_1_0 = 1;

1779

const double psitilde_bs_1_1 = 2.5*y + 1.5;

1780

const double psitilde_bs_2_0 = 1;

1781

1782

// Compute basisvalues

1783

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

1784

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

1785

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

1786

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

1787

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

1788

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

1789

1790

// Table(s) of coefficients

1791

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

1792

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

1793

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

1794

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

1795

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

1796

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

1797

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

1798

1799

// Interesting (new) part

1800

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

1801

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

1802

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

1803

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

1804

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

1805

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

1806

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

1807

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

1808

1809

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

1810

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

1811

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

1812

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

1813

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

1814

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

1815

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

1816

1817

// Compute reference derivatives

1818

// Declare pointer to array of derivatives on FIAT element

1819

double *derivatives = new double [num_derivatives];

1820

1821

// Declare coefficients

1822

double coeff0_0 = 0;

1823

double coeff0_1 = 0;

1824

double coeff0_2 = 0;

1825

double coeff0_3 = 0;

1826

double coeff0_4 = 0;

1827

double coeff0_5 = 0;

1828

1829

// Declare new coefficients

1830

double new_coeff0_0 = 0;

1831

double new_coeff0_1 = 0;

1832

double new_coeff0_2 = 0;

1833

double new_coeff0_3 = 0;

1834

double new_coeff0_4 = 0;

1835

double new_coeff0_5 = 0;

1836

1837

// Loop possible derivatives

1838

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

1839

{

1840

// Get values from coefficients array

1841

new_coeff0_0 = coefficients0[dof][0];

1842

new_coeff0_1 = coefficients0[dof][1];

1843

new_coeff0_2 = coefficients0[dof][2];

1844

new_coeff0_3 = coefficients0[dof][3];

1845

new_coeff0_4 = coefficients0[dof][4];

1846

new_coeff0_5 = coefficients0[dof][5];

1847

1848

// Loop derivative order

1849

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

1850

{

1851

// Update old coefficients

1852

coeff0_0 = new_coeff0_0;

1853

coeff0_1 = new_coeff0_1;

1854

coeff0_2 = new_coeff0_2;

1855

coeff0_3 = new_coeff0_3;

1856

coeff0_4 = new_coeff0_4;

1857

coeff0_5 = new_coeff0_5;

1858

1859

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

1860

{

1861

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

1862

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

1863

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

1864

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

1865

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

1866

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

1867

}

1868

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

1869

{

1870

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

1871

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

1872

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

1873

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

1874

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

1875

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

1876

}

1877

1878

}

1879

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

1880

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;

1881

}

1882

1883

// Transform derivatives back to physical element

1884

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

1885

{

1886

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

1887

{

1888

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

1889

}

1890

}

1891

// Delete pointer to array of derivatives on FIAT element

1892

delete [] derivatives;

1893

1894

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

1895

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

1896

{

1897

delete [] combinations[row];

1898

delete [] transform[row];

1899

}

1900

1901

delete [] combinations;

1902

delete [] transform;

1903

}

1904

1905

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

1906

virtual void evaluate_basis_derivatives_all(unsigned int n,

1907

double* values,

1908

const double* coordinates,

1909

const ufc::cell& c) const

1910

{

1911

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

1912

}

1913

1914

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

1915

virtual double evaluate_dof(unsigned int i,

1916

const ufc::function& f,

1917

const ufc::cell& c) const

1918

{

1919

// The reference points, direction and weights:

1920

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

1921

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

1922

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

1923

1924

const double * const * x = c.coordinates;

1925

double result = 0.0;

1926

// Iterate over the points:

1927

// Evaluate basis functions for affine mapping

1928

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

1929

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

1930

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

1931

1932

// Compute affine mapping y = F(X)

1933

double y[2];

1934

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

1935

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

1936

1937

// Evaluate function at physical points

1938

double values[1];

1939

f.evaluate(values, y, c);

1940

1941

// Map function values using appropriate mapping

1942

// Affine map: Do nothing

1943

1944

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

1945

1946

// Take directional components

1947

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

1948

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

1949

// Multiply by weights

1950

result *= W[i][0];

1951

1952

return result;

1953

}

1954

1955

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

1956

virtual void evaluate_dofs(double* values,

1957

const ufc::function& f,

1958

const ufc::cell& c) const

1959

{

1960

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

1961

}

1962

1963

/// Interpolate vertex values from dof values

1964

virtual void interpolate_vertex_values(double* vertex_values,

1965

const double* dof_values,

1966

const ufc::cell& c) const

1967

{

1968

// Evaluate at vertices and use affine mapping

1969

vertex_values[0] = dof_values[0];

1970

vertex_values[1] = dof_values[1];

1971

vertex_values[2] = dof_values[2];

1972

}

1973

1974

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

1975

virtual unsigned int num_sub_elements() const

1976

{

1977

return 1;

1978

}

1979

1980

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

1981

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

1982

{

1983

return new UFC_Poisson2D_2LinearForm_finite_element_0();

1984

}

1985

1986

};

1987

1988

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

1989

1990

class UFC_Poisson2D_2LinearForm_finite_element_1: public ufc::finite_element

1991

{

1992

public:

1993

1994

/// Constructor

1995

UFC_Poisson2D_2LinearForm_finite_element_1() : ufc::finite_element()

1996

{

1997

// Do nothing

1998

}

1999

2000

/// Destructor

2001

virtual ~UFC_Poisson2D_2LinearForm_finite_element_1()

2002

{

2003

// Do nothing

2004

}

2005

2006

/// Return a string identifying the finite element

2007

virtual const char* signature() const

2008

{

2009

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

2010

}

2011

2012

/// Return the cell shape

2013

virtual ufc::shape cell_shape() const

2014

{

2015

return ufc::triangle;

2016

}

2017

2018

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

2019

virtual unsigned int space_dimension() const

2020

{

2021

return 6;

2022

}

2023

2024

/// Return the rank of the value space

2025

virtual unsigned int value_rank() const

2026

{

2027

return 0;

2028

}

2029

2030

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

2031

virtual unsigned int value_dimension(unsigned int i) const

2032

{

2033

return 1;

2034

}

2035

2036

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

2037

virtual void evaluate_basis(unsigned int i,

2038

double* values,

2039

const double* coordinates,

2040

const ufc::cell& c) const

2041

{

2042

// Extract vertex coordinates

2043

const double * const * element_coordinates = c.coordinates;

2044

2045

// Compute Jacobian of affine map from reference cell

2046

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

2047

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

2048

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

2049

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

2050

2051

// Compute determinant of Jacobian

2052

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

2053

2054

// Compute inverse of Jacobian

2055

2056

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

2057

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

2058

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

2059

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

2060

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

2061

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

2062

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

2063

2064

// Map coordinates to the reference square

2065

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

2066

x = -1.0;

2067

else

2068

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

2069

y = 2.0*y - 1.0;

2070

2071

// Reset values

2072

*values = 0;

2073

2074

// Map degree of freedom to element degree of freedom

2075

const unsigned int dof = i;

2076

2077

// Generate scalings

2078

const double scalings_y_0 = 1;

2079

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

2080

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

2081

2082

// Compute psitilde_a

2083

const double psitilde_a_0 = 1;

2084

const double psitilde_a_1 = x;

2085

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

2086

2087

// Compute psitilde_bs

2088

const double psitilde_bs_0_0 = 1;

2089

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2090

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;

2091

const double psitilde_bs_1_0 = 1;

2092

const double psitilde_bs_1_1 = 2.5*y + 1.5;

2093

const double psitilde_bs_2_0 = 1;

2094

2095

// Compute basisvalues

2096

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

2097

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

2098

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2099

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

2100

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

2101

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

2102

2103

// Table(s) of coefficients

2104

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

2105

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

2106

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

2107

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

2108

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

2109

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

2110

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

2111

2112

// Extract relevant coefficients

2113

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

2114

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

2115

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

2116

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

2117

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

2118

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

2119

2120

// Compute value(s)

2121

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

2122

}

2123

2124

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

2125

virtual void evaluate_basis_all(double* values,

2126

const double* coordinates,

2127

const ufc::cell& c) const

2128

{

2129

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

2130

}

2131

2132

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

2133

virtual void evaluate_basis_derivatives(unsigned int i,

2134

unsigned int n,

2135

double* values,

2136

const double* coordinates,

2137

const ufc::cell& c) const

2138

{

2139

// Extract vertex coordinates

2140

const double * const * element_coordinates = c.coordinates;

2141

2142

// Compute Jacobian of affine map from reference cell

2143

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

2144

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

2145

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

2146

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

2147

2148

// Compute determinant of Jacobian

2149

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

2150

2151

// Compute inverse of Jacobian

2152

2153

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

2154

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

2155

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

2156

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

2157

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

2158

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

2159

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

2160

2161

// Map coordinates to the reference square

2162

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

2163

x = -1.0;

2164

else

2165

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

2166

y = 2.0*y - 1.0;

2167

2168

// Compute number of derivatives

2169

unsigned int num_derivatives = 1;

2170

2171

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

2172

num_derivatives *= 2;

2173

2174

2175

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

2176

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

2177

2178

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

2179

{

2180

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

2181

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

2182

combinations[j][k] = 0;

2183

}

2184

2185

// Generate combinations of derivatives

2186

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

2187

{

2188

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

2189

{

2190

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

2191

{

2192

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

2193

combinations[row][col] = 0;

2194

else

2195

{

2196

combinations[row][col] += 1;

2197

break;

2198

}

2199

}

2200

}

2201

}

2202

2203

// Compute inverse of Jacobian

2204

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

2205

2206

// Declare transformation matrix

2207

// Declare pointer to two dimensional array and initialise

2208

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

2209

2210

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

2211

{

2212

transform[j] = new double [num_derivatives];

2213

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

2214

transform[j][k] = 1;

2215

}

2216

2217

// Construct transformation matrix

2218

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

2219

{

2220

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

2221

{

2222

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

2223

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

2224

}

2225

}

2226

2227

// Reset values

2228

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

2229

values[j] = 0;

2230

2231

// Map degree of freedom to element degree of freedom

2232

const unsigned int dof = i;

2233

2234

// Generate scalings

2235

const double scalings_y_0 = 1;

2236

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

2237

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

2238

2239

// Compute psitilde_a

2240

const double psitilde_a_0 = 1;

2241

const double psitilde_a_1 = x;

2242

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

2243

2244

// Compute psitilde_bs

2245

const double psitilde_bs_0_0 = 1;

2246

const double psitilde_bs_0_1 = 1.5*y + 0.5;

2247

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;

2248

const double psitilde_bs_1_0 = 1;

2249

const double psitilde_bs_1_1 = 2.5*y + 1.5;

2250

const double psitilde_bs_2_0 = 1;

2251

2252

// Compute basisvalues

2253

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

2254

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

2255

const double basisvalue2 = psitilde_a_0*scalings_y_0*psitilde_bs_0_1;

2256

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

2257

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

2258

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

2259

2260

// Table(s) of coefficients

2261

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

2262

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

2263

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

2264

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

2265

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

2266

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

2267

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

2268

2269

// Interesting (new) part

2270

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

2271

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

2272

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

2273

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

2274

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

2275

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

2276

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

2277

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

2278

2279

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

2280

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

2281

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

2282

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

2283

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

2284

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

2285

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

2286

2287

// Compute reference derivatives

2288

// Declare pointer to array of derivatives on FIAT element

2289

double *derivatives = new double [num_derivatives];

2290

2291

// Declare coefficients

2292

double coeff0_0 = 0;

2293

double coeff0_1 = 0;

2294

double coeff0_2 = 0;

2295

double coeff0_3 = 0;

2296

double coeff0_4 = 0;

2297

double coeff0_5 = 0;

2298

2299

// Declare new coefficients

2300

double new_coeff0_0 = 0;

2301

double new_coeff0_1 = 0;

2302

double new_coeff0_2 = 0;

2303

double new_coeff0_3 = 0;

2304

double new_coeff0_4 = 0;

2305

double new_coeff0_5 = 0;

2306

2307

// Loop possible derivatives

2308

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

2309

{

2310

// Get values from coefficients array

2311

new_coeff0_0 = coefficients0[dof][0];

2312

new_coeff0_1 = coefficients0[dof][1];

2313

new_coeff0_2 = coefficients0[dof][2];

2314

new_coeff0_3 = coefficients0[dof][3];

2315

new_coeff0_4 = coefficients0[dof][4];

2316

new_coeff0_5 = coefficients0[dof][5];

2317

2318

// Loop derivative order

2319

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

2320

{

2321

// Update old coefficients

2322

coeff0_0 = new_coeff0_0;

2323

coeff0_1 = new_coeff0_1;

2324

coeff0_2 = new_coeff0_2;

2325

coeff0_3 = new_coeff0_3;

2326

coeff0_4 = new_coeff0_4;

2327

coeff0_5 = new_coeff0_5;

2328

2329

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

2330

{

2331

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

2332

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

2333

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

2334

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

2335

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

2336

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

2337

}

2338

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

2339

{

2340

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

2341

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

2342

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

2343

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

2344

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

2345

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

2346

}

2347

2348

}

2349

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

2350

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;

2351

}

2352

2353

// Transform derivatives back to physical element

2354

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

2355

{

2356

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

2357

{

2358

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

2359

}

2360

}

2361

// Delete pointer to array of derivatives on FIAT element

2362

delete [] derivatives;

2363

2364

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

2365

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

2366

{

2367

delete [] combinations[row];

2368

delete [] transform[row];

2369

}

2370

2371

delete [] combinations;

2372

delete [] transform;

2373

}

2374

2375

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

2376

virtual void evaluate_basis_derivatives_all(unsigned int n,

2377

double* values,

2378

const double* coordinates,

2379

const ufc::cell& c) const

2380

{

2381

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

2382

}

2383

2384

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

2385

virtual double evaluate_dof(unsigned int i,

2386

const ufc::function& f,

2387

const ufc::cell& c) const

2388

{

2389

// The reference points, direction and weights:

2390

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

2391

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

2392

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

2393

2394

const double * const * x = c.coordinates;

2395

double result = 0.0;

2396

// Iterate over the points:

2397

// Evaluate basis functions for affine mapping

2398

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

2399

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

2400

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

2401

2402

// Compute affine mapping y = F(X)

2403

double y[2];

2404

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

2405

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

2406

2407

// Evaluate function at physical points

2408

double values[1];

2409

f.evaluate(values, y, c);

2410

2411

// Map function values using appropriate mapping

2412

// Affine map: Do nothing

2413

2414

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

2415

2416

// Take directional components

2417

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

2418

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

2419

// Multiply by weights

2420

result *= W[i][0];

2421

2422

return result;

2423

}

2424

2425

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

2426

virtual void evaluate_dofs(double* values,

2427

const ufc::function& f,

2428

const ufc::cell& c) const

2429

{

2430

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

2431

}

2432

2433

/// Interpolate vertex values from dof values

2434

virtual void interpolate_vertex_values(double* vertex_values,

2435

const double* dof_values,

2436

const ufc::cell& c) const

2437

{

2438

// Evaluate at vertices and use affine mapping

2439

vertex_values[0] = dof_values[0];

2440

vertex_values[1] = dof_values[1];

2441

vertex_values[2] = dof_values[2];

2442

}

2443

2444

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

2445

virtual unsigned int num_sub_elements() const

2446

{

2447

return 1;

2448

}

2449

2450

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

2451

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

2452

{

2453

return new UFC_Poisson2D_2LinearForm_finite_element_1();

2454

}

2455

2456

};

2457

2458

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

2459

/// degrees of freedom (dofs).

2460

2461

class UFC_Poisson2D_2LinearForm_dof_map_0: public ufc::dof_map

2462

{

2463

private:

2464

2465

unsigned int __global_dimension;

2466

2467

public:

2468

2469

/// Constructor

2470

UFC_Poisson2D_2LinearForm_dof_map_0() : ufc::dof_map()

2471

{

2472

__global_dimension = 0;

2473

}

2474

2475

/// Destructor

2476

virtual ~UFC_Poisson2D_2LinearForm_dof_map_0()

2477

{

2478

// Do nothing

2479

}

2480

2481

/// Return a string identifying the dof map

2482

virtual const char* signature() const

2483

{

2484

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

2485

}

2486

2487

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

2488

virtual bool needs_mesh_entities(unsigned int d) const

2489

{

2490

switch ( d )

2491

{

2492

case 0:

2493

return true;

2494

break;

2495

case 1:

2496

return true;

2497

break;

2498

case 2:

2499

return false;

2500

break;

2501

}

2502

return false;

2503

}

2504

2505

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

2506

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

2507

{

2508

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

2509

return false;

2510

}

2511

2512

/// Initialize dof map for given cell

2513

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

2514

const ufc::cell& c)

2515

{

2516

// Do nothing

2517

}

2518

2519

/// Finish initialization of dof map for cells

2520

virtual void init_cell_finalize()

2521

{

2522

// Do nothing

2523

}

2524

2525

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

2526

virtual unsigned int global_dimension() const

2527

{

2528

return __global_dimension;

2529

}

2530

2531

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

2532

virtual unsigned int local_dimension() const

2533

{

2534

return 6;

2535

}

2536

2537

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

2538

virtual unsigned int geometric_dimension() const

2539

{

2540

return 2;

2541

}

2542

2543

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

2544

virtual unsigned int num_facet_dofs() const

2545

{

2546

return 3;

2547

}

2548

2549

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

2550

virtual unsigned int num_entity_dofs(unsigned int d) const

2551

{

2552

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

2553

}

2554

2555

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

2556

virtual void tabulate_dofs(unsigned int* dofs,

2557

const ufc::mesh& m,

2558

const ufc::cell& c) const

2559

{

2560

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

2561

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

2562

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

2563

unsigned int offset = m.num_entities[0];

2564

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

2565

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

2566

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

2567

}

2568

2569

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

2570

virtual void tabulate_facet_dofs(unsigned int* dofs,

2571

unsigned int facet) const

2572

{

2573

switch ( facet )

2574

{

2575

case 0:

2576

dofs[0] = 1;

2577

dofs[1] = 2;

2578

dofs[2] = 3;

2579

break;

2580

case 1:

2581

dofs[0] = 0;

2582

dofs[1] = 2;

2583

dofs[2] = 4;

2584

break;

2585

case 2:

2586

dofs[0] = 0;

2587

dofs[1] = 1;

2588

dofs[2] = 5;

2589

break;

2590

}

2591

}

2592

2593

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

2594

virtual void tabulate_entity_dofs(unsigned int* dofs,

2595

unsigned int d, unsigned int i) const

2596

{

2597

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

2598

}

2599

2600

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

2601

virtual void tabulate_coordinates(double** coordinates,

2602

const ufc::cell& c) const

2603

{

2604

const double * const * x = c.coordinates;

2605

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

2606

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

2607

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

2608

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

2609

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

2610

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

2611

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

2612

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

2613

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

2614

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

2615

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

2616

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

2617

}

2618

2619

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

2620

virtual unsigned int num_sub_dof_maps() const

2621

{

2622

return 1;

2623

}

2624

2625

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

2626

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

2627

{

2628

return new UFC_Poisson2D_2LinearForm_dof_map_0();

2629

}

2630

2631

};

2632

2633

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

2634

/// degrees of freedom (dofs).

2635

2636

class UFC_Poisson2D_2LinearForm_dof_map_1: public ufc::dof_map

2637

{

2638

private:

2639

2640

unsigned int __global_dimension;

2641

2642

public:

2643

2644

/// Constructor

2645

UFC_Poisson2D_2LinearForm_dof_map_1() : ufc::dof_map()

2646

{

2647

__global_dimension = 0;

2648

}

2649

2650

/// Destructor

2651

virtual ~UFC_Poisson2D_2LinearForm_dof_map_1()

2652

{

2653

// Do nothing

2654

}

2655

2656

/// Return a string identifying the dof map

2657

virtual const char* signature() const

2658

{

2659

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

2660

}

2661

2662

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

2663

virtual bool needs_mesh_entities(unsigned int d) const

2664

{

2665

switch ( d )

2666

{

2667

case 0:

2668

return true;

2669

break;

2670

case 1:

2671

return true;

2672

break;

2673

case 2:

2674

return false;

2675

break;

2676

}

2677

return false;

2678

}

2679

2680

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

2681

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

2682

{

2683

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

2684

return false;

2685

}

2686

2687

/// Initialize dof map for given cell

2688

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

2689

const ufc::cell& c)

2690

{

2691

// Do nothing

2692

}

2693

2694

/// Finish initialization of dof map for cells

2695

virtual void init_cell_finalize()

2696

{

2697

// Do nothing

2698

}

2699

2700

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

2701

virtual unsigned int global_dimension() const

2702

{

2703

return __global_dimension;

2704

}

2705

2706

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

2707

virtual unsigned int local_dimension() const

2708

{

2709

return 6;

2710

}

2711

2712

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

2713

virtual unsigned int geometric_dimension() const

2714

{

2715

return 2;

2716

}

2717

2718

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

2719

virtual unsigned int num_facet_dofs() const

2720

{

2721

return 3;

2722

}

2723

2724

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

2725

virtual unsigned int num_entity_dofs(unsigned int d) const

2726

{

2727

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

2728

}

2729

2730

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

2731

virtual void tabulate_dofs(unsigned int* dofs,

2732

const ufc::mesh& m,

2733

const ufc::cell& c) const

2734

{

2735

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

2736

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

2737

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

2738

unsigned int offset = m.num_entities[0];

2739

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

2740

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

2741

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

2742

}

2743

2744

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

2745

virtual void tabulate_facet_dofs(unsigned int* dofs,

2746

unsigned int facet) const

2747

{

2748

switch ( facet )

2749

{

2750

case 0:

2751

dofs[0] = 1;

2752

dofs[1] = 2;

2753

dofs[2] = 3;

2754

break;

2755

case 1:

2756

dofs[0] = 0;

2757

dofs[1] = 2;

2758

dofs[2] = 4;

2759

break;

2760

case 2:

2761

dofs[0] = 0;

2762

dofs[1] = 1;

2763

dofs[2] = 5;

2764

break;

2765

}

2766

}

2767

2768

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

2769

virtual void tabulate_entity_dofs(unsigned int* dofs,

2770

unsigned int d, unsigned int i) const

2771

{

2772

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

2773

}

2774

2775

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

2776

virtual void tabulate_coordinates(double** coordinates,

2777

const ufc::cell& c) const

2778

{

2779

const double * const * x = c.coordinates;

2780

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

2781

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

2782

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

2783

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

2784

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

2785

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

2786

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

2787

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

2788

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

2789

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

2790

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

2791

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

2792

}

2793

2794

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

2795

virtual unsigned int num_sub_dof_maps() const

2796

{

2797

return 1;

2798

}

2799

2800

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

2801

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

2802

{

2803

return new UFC_Poisson2D_2LinearForm_dof_map_1();

2804

}

2805

2806

};

2807

2808

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

2809

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

2810

/// the integral over a cell.

2811

2812

class UFC_Poisson2D_2LinearForm_cell_integral_0: public ufc::cell_integral

2813

{

2814

public:

2815

2816

/// Constructor

2817

UFC_Poisson2D_2LinearForm_cell_integral_0() : ufc::cell_integral()

2818

{

2819

// Do nothing

2820

}

2821

2822

/// Destructor

2823

virtual ~UFC_Poisson2D_2LinearForm_cell_integral_0()

2824

{

2825

// Do nothing

2826

}

2827

2828

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

2829

virtual void tabulate_tensor(double* A,

2830

const double * const * w,

2831

const ufc::cell& c) const

2832

{

2833

// Extract vertex coordinates

2834

const double * const * x = c.coordinates;

2835

2836

// Compute Jacobian of affine map from reference cell

2837

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

2838

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

2839

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

2840

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

2841

2842

// Compute determinant of Jacobian

2843

double detJ = J_00*J_11 - J_01*J_10;

2844

2845

// Compute inverse of Jacobian

2846

2847

// Set scale factor

2848

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

2849

2850

// Number of operations to compute element tensor = 48

2851

// Compute coefficients

2852

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

2853

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

2854

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

2855

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

2856

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

2857

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

2858

2859

// Compute geometry tensors

2860

// Number of operations to compute decalrations = 6

2861

const double G0_0 = det*c0_0_0_0;

2862

const double G0_1 = det*c0_0_0_1;

2863

const double G0_2 = det*c0_0_0_2;

2864

const double G0_3 = det*c0_0_0_3;

2865

const double G0_4 = det*c0_0_0_4;

2866

const double G0_5 = det*c0_0_0_5;

2867

2868

// Compute element tensor

2869

// Number of operations to compute tensor = 42

2870

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

2871

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

2872

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

2873

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

2874

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

2875

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

2876

}

2877

2878

};

2879

2880

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

2881

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

2882

/// mapping

2883

///

2884

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

2885

///

2886

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

2887

/// global tensor A is defined by

2888

///

2889

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

2890

///

2891

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

2892

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

2893

/// fixed functions (coefficients).

2894

2895

class UFC_Poisson2D_2LinearForm: public ufc::form

2896

{

2897

public:

2898

2899

/// Constructor

2900

UFC_Poisson2D_2LinearForm() : ufc::form()

2901

{

2902

// Do nothing

2903

}

2904

2905

/// Destructor

2906

virtual ~UFC_Poisson2D_2LinearForm()

2907

{

2908

// Do nothing

2909

}

2910

2911

/// Return a string identifying the form

2912

virtual const char* signature() const

2913

{

2914

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)";

2915

}

2916

2917

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

2918

virtual unsigned int rank() const

2919

{

2920

return 1;

2921

}

2922

2923

/// Return the number of coefficients (n)

2924

virtual unsigned int num_coefficients() const

2925

{

2926

return 1;

2927

}

2928

2929

/// Return the number of cell integrals

2930

virtual unsigned int num_cell_integrals() const

2931

{

2932

return 1;

2933

}

2934

2935

/// Return the number of exterior facet integrals

2936

virtual unsigned int num_exterior_facet_integrals() const

2937

{

2938

return 0;

2939

}

2940

2941

/// Return the number of interior facet integrals

2942

virtual unsigned int num_interior_facet_integrals() const

2943

{

2944

return 0;

2945

}

2946

2947

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

2948

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

2949

{

2950

switch ( i )

2951

{

2952

case 0:

2953

return new UFC_Poisson2D_2LinearForm_finite_element_0();

2954

break;

2955

case 1:

2956

return new UFC_Poisson2D_2LinearForm_finite_element_1();

2957

break;

2958

}

2959

return 0;

2960

}

2961

2962

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

2963

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

2964

{

2965

switch ( i )

2966

{

2967

case 0:

2968

return new UFC_Poisson2D_2LinearForm_dof_map_0();

2969

break;

2970

case 1:

2971

return new UFC_Poisson2D_2LinearForm_dof_map_1();

2972

break;

2973

}

2974

return 0;

2975

}

2976

2977

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

2978

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

2979

{

2980

return new UFC_Poisson2D_2LinearForm_cell_integral_0();

2981

}

2982

2983

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

2984

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

2985

{

2986

return 0;

2987

}

2988

2989

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

2990

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

2991

{

2992

return 0;

2993

}

2994

2995

};

2996

2997

// DOLFIN wrappers

2998

2999

#include <dolfin/fem/Form.h>

3000

#include <dolfin/fem/FiniteElement.h>

3001

#include <dolfin/fem/DofMap.h>

3002

#include <dolfin/function/Coefficient.h>

3003

#include <dolfin/function/Function.h>

3004

#include <dolfin/function/FunctionSpace.h>

3005

3006

class Poisson2D_2BilinearFormFunctionSpace0 : public dolfin::FunctionSpace

3007

{

3008

public:

3009

3010

Poisson2D_2BilinearFormFunctionSpace0(const dolfin::Mesh& mesh)

3011

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3012

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_2LinearForm_finite_element_1()))),

3013

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_2LinearForm_dof_map_1()), mesh)))

3014

{

3015

// Do nothing

3016

}

3017

3018

};

3019

3020

class Poisson2D_2BilinearFormFunctionSpace1 : public dolfin::FunctionSpace

3021

{

3022

public:

3023

3024

Poisson2D_2BilinearFormFunctionSpace1(const dolfin::Mesh& mesh)

3025

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3026

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_2LinearForm_finite_element_1()))),

3027

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_2LinearForm_dof_map_1()), mesh)))

3028

{

3029

// Do nothing

3030

}

3031

3032

};

3033

3034

class Poisson2D_2LinearFormFunctionSpace0 : public dolfin::FunctionSpace

3035

{

3036

public:

3037

3038

Poisson2D_2LinearFormFunctionSpace0(const dolfin::Mesh& mesh)

3039

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3040

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_2LinearForm_finite_element_1()))),

3041

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_2LinearForm_dof_map_1()), mesh)))

3042

{

3043

// Do nothing

3044

}

3045

3046

};

3047

3048

class Poisson2D_2LinearFormCoefficientSpace0 : public dolfin::FunctionSpace

3049

{

3050

public:

3051

3052

Poisson2D_2LinearFormCoefficientSpace0(const dolfin::Mesh& mesh)

3053

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3054

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_2LinearForm_finite_element_1()))),

3055

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_2LinearForm_dof_map_1()), mesh)))

3056

{

3057

// Do nothing

3058

}

3059

3060

};

3061

3062

class Poisson2D_2TestSpace : public dolfin::FunctionSpace

3063

{

3064

public:

3065

3066

Poisson2D_2TestSpace(const dolfin::Mesh& mesh)

3067

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3068

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_2LinearForm_finite_element_1()))),

3069

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_2LinearForm_dof_map_1()), mesh)))

3070

{

3071

// Do nothing

3072

}

3073

3074

};

3075

3076

class Poisson2D_2TrialSpace : public dolfin::FunctionSpace

3077

{

3078

public:

3079

3080

Poisson2D_2TrialSpace(const dolfin::Mesh& mesh)

3081

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3082

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_2LinearForm_finite_element_1()))),

3083

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_2LinearForm_dof_map_1()), mesh)))

3084

{

3085

// Do nothing

3086

}

3087

3088

};

3089

3090

class Poisson2D_2CoefficientSpace : public dolfin::FunctionSpace

3091

{

3092

public:

3093

3094

Poisson2D_2CoefficientSpace(const dolfin::Mesh& mesh)

3095

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3096

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_2LinearForm_finite_element_1()))),

3097

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_2LinearForm_dof_map_1()), mesh)))

3098

{

3099

// Do nothing

3100

}

3101

3102

};

3103

3104

class Poisson2D_2FunctionSpace : public dolfin::FunctionSpace

3105

{

3106

public:

3107

3108

Poisson2D_2FunctionSpace(const dolfin::Mesh& mesh)

3109

: dolfin::FunctionSpace(boost::shared_ptr<const dolfin::Mesh>(&mesh, dolfin::NoDeleter<const dolfin::Mesh>()),

3110

boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new UFC_Poisson2D_2LinearForm_finite_element_1()))),

3111

boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new UFC_Poisson2D_2LinearForm_dof_map_1()), mesh)))

3112

{

3113

// Do nothing

3114

}

3115

3116

};

3117

3118

class Poisson2D_2BilinearForm : public dolfin::Form

3119

{

3120

public:

3121

3122

// Create form on given function space(s)

3123

Poisson2D_2BilinearForm(const dolfin::FunctionSpace& V0, const dolfin::FunctionSpace& V1) : dolfin::Form()

3124

{

3125

boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());

3126

_function_spaces.push_back(_V0);

3127

boost::shared_ptr<const dolfin::FunctionSpace> _V1(&V1, dolfin::NoDeleter<const dolfin::FunctionSpace>());

3128

_function_spaces.push_back(_V1);

3129

3130

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_2BilinearForm());

3131

}

3132

3133

// Create form on given function space(s) (shared data)

3134

Poisson2D_2BilinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0, boost::shared_ptr<const dolfin::FunctionSpace> V1) : dolfin::Form()

3135

{

3136

_function_spaces.push_back(V0);

3137

_function_spaces.push_back(V1);

3138

3139

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_2BilinearForm());

3140

}

3141

3142

// Destructor

3143

~Poisson2D_2BilinearForm() {}

3144

3145

};

3146

3147

class Poisson2D_2LinearFormCoefficient0 : public dolfin::Coefficient

3148

{

3149

public:

3150

3151

// Constructor

3152

Poisson2D_2LinearFormCoefficient0(dolfin::Form& form) : dolfin::Coefficient(form) {}

3153

3154

// Destructor

3155

~Poisson2D_2LinearFormCoefficient0() {}

3156

3157

// Attach function to coefficient

3158

const Poisson2D_2LinearFormCoefficient0& operator= (dolfin::Function& v)

3159

{

3160

attach(v);

3161

return *this;

3162

}

3163

3164

/// Create function space for coefficient

3165

const dolfin::FunctionSpace* create_function_space() const

3166

{

3167

return new Poisson2D_2LinearFormCoefficientSpace0(form.mesh());

3168

}

3169

3170

/// Return coefficient number

3171

dolfin::uint number() const

3172

{

3173

return 0;

3174

}

3175

3176

/// Return coefficient name

3177

virtual std::string name() const

3178

{

3179

return "f";

3180

}

3181

3182

};

3183

class Poisson2D_2LinearForm : public dolfin::Form

3184

{

3185

public:

3186

3187

// Create form on given function space(s)

3188

Poisson2D_2LinearForm(const dolfin::FunctionSpace& V0) : dolfin::Form(), f(*this)

3189

{

3190

boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());

3191

_function_spaces.push_back(_V0);

3192

3193

_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));

3194

3195

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_2LinearForm());

3196

}

3197

3198

// Create form on given function space(s) (shared data)

3199

Poisson2D_2LinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0) : dolfin::Form(), f(*this)

3200

{

3201

_function_spaces.push_back(V0);

3202

3203

_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));

3204

3205

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_2LinearForm());

3206

}

3207

3208

// Create form on given function space(s) with given coefficient(s)

3209

Poisson2D_2LinearForm(const dolfin::FunctionSpace& V0, dolfin::Function& w0) : dolfin::Form(), f(*this)

3210

{

3211

boost::shared_ptr<const dolfin::FunctionSpace> _V0(&V0, dolfin::NoDeleter<const dolfin::FunctionSpace>());

3212

_function_spaces.push_back(_V0);

3213

3214

_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));

3215

3216

this->f = w0;

3217

3218

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_2LinearForm());

3219

}

3220

3221

// Create form on given function space(s) with given coefficient(s) (shared data)

3222

Poisson2D_2LinearForm(boost::shared_ptr<const dolfin::FunctionSpace> V0, dolfin::Function& w0) : dolfin::Form(), f(*this)

3223

{

3224

_function_spaces.push_back(V0);

3225

3226

_coefficients.push_back(boost::shared_ptr<const dolfin::Function>(static_cast<const dolfin::Function*>(0)));

3227

3228

this->f = w0;

3229

3230

_ufc_form = boost::shared_ptr<const ufc::form>(new UFC_Poisson2D_2LinearForm());

3231

}

3232

3233

// Destructor

3234

~Poisson2D_2LinearForm() {}

3235

3236

// Coefficients

3237

Poisson2D_2LinearFormCoefficient0 f;

3238

3239

};

3240

3241

#endif

Older »