3-d Poisson in c++

Question

Hi, I just got started with FEniCS and I'm trying to solve a simple 3-D Poisson equation. After modifying the c++ demon code in neumann-poisson, I constantly get this error:
"
Mesh cell type (tetrahedra) does not match cell type of form
"
Basically, I only changed the mesh from the original 2-D demo and the elements in the ufl file:

[Poisson.ufl]
V=FiniteElement("Lagrange", tetrahedra,1)
R=FiniteElement("R",tetrahedra,1)

[main.cpp]
BoxMesh mesh(0,0,0,10,10,10,20,20,20);

I guess there something I have missed. Would someone please help me point out how I could solve this? I apologize if this question is too naive, but I really have no previous knowledge how dolfin works. Thank you very much.

Answer 1

Hi, try changing the ufl file to

V = FiniteElement("Lagrange", tetrahedron, 1)                                    
R = FiniteElement("R", tetrahedron, 0) 

The RealSpace can't have degree > 1 and the cell type in 3d is tetrahedron.

Comments

Hi, MiroK,
Thank you for this suggestion. I tried with this change, but the same error still occurs. I found this error is raised in the following lines from AssemblerBase.cpp:

    if (mesh.type().cell_type() == CellType::tetrahedron && element->cell_shape() != ufc::tetrahedron)
{
  dolfin_error("AssemblerBase.cpp",
               "assemble form",
               "Mesh cell type (tetrahedra) does not match cell type of form");
}

It seem somehow the program does not recognize tetrahedron as the element celll_shape. Would please see what is happening here? Thank you.

---

Could you post the entire code?

---

Also, after you change the ufl file make sure to run

ffc -l dolfin Poisson.ufl

---

Hi, MiroK,
Thank you very much. This solved the problem. I thought the ffc command were included in cmake and redid cmake only each time. The code is as following. It now works without the Dirichlet boundary condition, although I haven't been able to verify the correctness of the solution yet. I'm also still trying to figure out how a Dirichlet BC in the middle can be added.

Poisson.ufl:

V = FiniteElement("Lagrange", tetrahedron, 1)
R = FiniteElement("R", tetrahedron, 0)
element = V * R
#element = FiniteElement("Lagrange", tetrahedra, 1)
#V = element
#R = element

(u, c) = TrialFunctions(element)
(v, d) = TestFunctions(element)
f = Coefficient(V)
g = Coefficient(V)

a = (inner(grad(u), grad(v)) + c*v + u*d)*dx
L = f*v*dx + g*v*ds

main.cpp:

#include <dolfin.h>
#include "Poisson.h"
//#include <cmath>

using namespace dolfin;

// Source term (right-hand side)
class Source : public Expression
{
    public:
    double zz, s0, sfluc, eeps;
    Source(double z, double si, double sf):Expression() {
        zz=z; s0=si; sfluc=sf; eeps=1e-20;
    }

  void eval(Array<double>& values, const Array<double>& x) const
  {
      if (abs(x[2]-zz)<1e-9) {
          double rv=-1.0+2.0*dolfin::rand();
          values[0]=s0+rv*sfluc;
      }
      else {
          values[0]=0;
      }
  }
};

// Boundary flux (Neumann boundary condition)
class Flux : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = 0;
  }
};
// (Dirichlet boundary condition)
class DirichletBoundary : public SubDomain {
    public:
        double zz, vol;
        DirichletBoundary(double z, double v) : SubDomain() {
            zz=z; vol=v;
        }
    bool inside(const Array<double>& x, bool on_boundary) const {
        return (x[2]-zz>=0) && (x[2]-zz<2e-9);
    }
};

int main()
{
    double epsilon=8.854e-12;
    double Lx=500e-9;
    double Ly=500e-9;
    double Lz=20e-9;
    double zTI=10e-9;
    double zele=10e-9;
    double Nx=20;
    double Ny=20;
    double Nz=20;
    double s0=0;
    double sfluc=1e13;
  // Create mesh and function space
  //BoxMesh mesh(0,0,0,Lx,Ly,Lz,Nx,Ny,Nz);
  BoxMesh mesh(0,0,0,20,20,20,20,20,20);
  //plot(mesh);
  //interactive();
  //std::getchar();
  Poisson::FunctionSpace V(mesh);
  //V.print_dofmap();
  cout << V.dim() << endl;
  /* BC */
  Constant v0(0.0);
  DirichletBoundary dBC(zele,1);
  DirichletBC dirbc(V,v0,dBC); 

  // Define variational problem
  Poisson::BilinearForm a(V, V);
  Poisson::LinearForm L(V);
  Source f(zTI,s0,sfluc/epsilon);
  Flux g;
  L.f = f;
  L.g = g;
  Function w(V);
  //LinearVariationalProblem problem(a, L, w);

  // Compute solution
  //LinearVariationalSolver solver(problem);
  //solver.solve();

    solve(a==L, w, dirbc);

  // Extract subfunction
  Function u = w[0];

  // Plot solution
  plot(u);
  interactive();

  return 0;
}