dolfin team mailing list archive

[Marie Rognes <question156272@xxxxxxxxxxxxxxxxxxxxx>]

Question #156272 on DOLFIN changed:
https://answers.launchpad.net/dolfin/+question/156272

    Status: Open => Answered

Marie Rognes proposed the following answer:

Oh, you are in c++, ok. Define the form in a UFL file, compile it with FFC, and then assemble the resulting form in DOLFIN. 
Sample code included below. (Not compiled/tested, all of DOLFIN installations are currently broken)

Some more comments:

(1) For more info on the UFL/FFC/DOLFIN interplay, consider looking at some of the cpp demos or reading for instance the DOLFIN chapter in the FEniCS book.
(2) This way of doing it should be more trivial (and perhaps more generally useful) that evaluating the quadrature yourself.

--- DensityIntegral.ufl ---
# Compile using FFC with 
# >> ffc -l dolfin -O DensityIntegral.ufl

# Define some element and the to-be given coefficient for density
cell = "tetrahedron"
V = FiniteElement("CG", cell, 2)
density = Coefficient(V)

# Define DG element
Q = FiniteElement("DG", cell, 0)
v = TestFunction(Q)

# Define form to be used to evaluate \int_T density for each cell T
L = density*v*dx

--- main.cpp ---

#include <dolfin.h>
#include "DensityIntegral.h"

using namespace dolfin;

// Density term (right-hand side)
class Density : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    double dx = x[0] - 0.5;
    double dy = x[1] - 0.5;
    values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
  }
};

int main()
{
  // Create mesh
  UnitCube mesh(4, 4, 4);

  // Create density function (just an example)
  Density density;

  // Create test space
  DensityIntegral::FunctionSpace Q(mesh);

  // Create linear form
  DensityIntegral::LinearForm L(Q);
  L.density = density;

  // Assemble form into vector
  Vector b;
  assemble(b, L);
  info(b);

  // Next steps ...

  // 1. Filter vector to yield MeshFunction of bools over cells
  // ('markers') indicating which to refine or not

  // 2. Call
  // new_mesh = refine(mesh, markers)
  // to refine mesh locally

  return 0;
}

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.