What is a clean way to generate interval,triangle, and tetra elements from one form file?

Question

I am trying to avoid copy pasting my UFL files for each element type, and then dealing with the different forms in my solver.

For example, I'd rather not carry Pressure1D, Pressure2D, and Pressure3D. Instead I would like to carry a 'base' Pressure which I can assign as one of the three. Is there a way to do this, perhaps a base class that I can use and assign appropriately?

EDIT

An example of a UFL file for divergence:

vector = VectorElement("Lagrange", etype, hvector)
scalar = FiniteElement("Lagrange", etype, hscalar)
q = Coefficient(vector)  
a = u*w*dx
L = (Dx(q[i],i))*w*dx  

I take that UFL file and copy it with sed to replace etype, hvector and hscalar with the element type (interval,triangle,tetrahedral), and hvector and hscalar with the interpolation order for the vector and scalar respectively. This templates out the UFL files for me nicely which I then run FFC upon to get cpp/h files.

What I would like to do is code the solver agnostic to which particular FunctionSpace is used. I found I can do a fair amount of it using boost::optional and use the base dolfin:: types. So I might have:

boost::optional<dolfin::Function> u;
boost::optional<dolfin::FunctionSpace> V;
switch(dimensions){
  case 1:
    V = DivergenceOneDim::FunctionSpace(mesh); break;
  case 2:
    V = DivergenceTwoDim::FunctionSpace(mesh); break;
  case 3:
    V = DivergenceThreeDim::FunctionSpace(mesh); break;
}
u = Function(*u)

This doesnt solve the problem of access 'q' without boxing/wrapping, though.

Comments

If UFL could include a file that would help a lot as well.

---

Please post a code example, and say whether you're using UFL via the DOLFIN Python interface, or if you're working directly with UFL input.

---

I updated the question, thank you!

Answer 1

Use:

Mesh UnitCubeMesh(8, 8 ,8);
const std::size_t gdim = mesh.geometry().dim();

boost::shared_ptr<FunctionSpace> V;
if (gdim == 1)
  V.reset(new DivergenceOneDim::FunctionSpace(mesh));
else if (gdim == 2)
  V.reset(new DivergenceTwoDim::FunctionSpace(mesh));
else if (gdim == 3)
  V.reset(new DivergenceThreeDim::FunctionSpace(mesh));

Function u(V);

Comments

Could you explain a bit about form? I would have expected V.reset(). Thanks again!

---

A bit of a minor thing, but if I were then to plot u, the plotter has some trouble. Could this be due to the use of dolfin::FunctionSpace in the constructor?

*** -------------------------------------------------------------------------
*** Error: Unable to plot().
*** Reason: The plottable is not compatible with the data.
*** Where: This error was encountered inside VTKPlotter.cpp.
*** Process: 0
*** -------------------------------------------------------------------------

---

Sorry, I changed syntax while answering and didn't update. I've edited my answer to change form to V.

---

Thank you again!

---

If I may piggy back a question on this answer... I am using this approach in a constructor for a class containing several functions. Once I leave the scope of the constructor I get a number of errors (the plotting error being one such error). I suspect its a scoping issue, but, I would think a dynamically allocated object to a shared_ptr would survive the constructor. The pointer is a member of the class. Any thoughts? Thanks in advance.

---

The mesh is probably going out of scope. Try:

boost::shared_ptr<Mesh> mesh(new UnitCubeMesh(8, 8 ,8));
const std::size_t gdim = mesh.geometry().dim();

boost::shared_ptr<FunctionSpace> V;
if (gdim == 1)
  V.reset(new DivergenceOneDim::FunctionSpace(mesh));
else if (gdim == 2)
  V.reset(new DivergenceTwoDim::FunctionSpace(mesh));
else if (gdim == 3)
  V.reset(new DivergenceThreeDim::FunctionSpace(mesh));

---

Sorry, please disregard I found it shortly after asking (seems to always be when I figure it out :)).

I had forgotten to use a reference to the Mesh and was doing a copy constructor on it. When the scope terminated that copied Mesh was disposed leaving Functions and FunctionSpaces that referenced a missing Mesh.

It so happened the first few error checks reported things relating to dimensions/functionspaces, so I thought it had to do with my trying to pass 1/2/3D forms.

---

I saw your answer after posting that last one - you are good thanks! :)