C++ interface - Coefficient has not been set.


315
views
0
12 months ago by
Hi,
I been working with the python interface for some time now but just recently I started to use the C++ interface. I ran into a problem with assigning a value to a Coefficient of a BilinearForm defined in a ufl file. Here is the example of the main.cpp and PeriodicProblem.ufl
vector = VectorElement("Lagrange", tetrahedron, 1)
scalar = FiniteElement("Lagrange", tetrahedron, 1)

u = TrialFunction(vector)
v = TestFunction(vector)

mu = Coefficient(scalar)
lmbda = 1e4
A = as_tensor([[[[2*mu+lmbda,0,0],[0,lmbda,0],[0,0,lmbda]],\
                [[0,mu,0],[mu,0,0],[0,0,0]],\
            [[0,0,mu],[0,0,0],[mu,0,0]]],\
            [[[0,mu,0],[mu,0,0],[0,0,0]],\
            [[lmbda,0,0],[0,2*mu+lmbda,0],[0,0,lmbda]],\
            [[0,0,0],[0,0,mu],[0,mu,0]]],\
            [[[0,0,mu],[0,0,0],[mu,0,0]],\
            [[0,0,0],[0,0,mu],[0,mu,0]],\
            [[lmbda,0,0],[0,lmbda,0],[0,0,2*mu+lmbda]]]])


i,j,k,l = indices(4)

a = A[i,j,k,l]*Dx(u[i],j)*Dx(v[k],l)*dx
L = A[0,0,i,j]*Dx(v[i],j)*dx​


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

using namespace dolfin;

class PeriodicBoundary : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return (near(x[0],0) or near(x[1],0) or near(x[2],0))
            and (not (
                    (near(x[0],0) and (near(x[1],1) or near(x[2],2))) or
                    (near(x[1],0) and (near(x[0],1.2) or near(x[2],2))) or
                    (near(x[2],0) and (near(x[0],1.2) or near(x[1],1)))
                    )) and on_boundary;
  }

  void map(const Array<double>& x, Array<double>& y) const
  {
    if(near(x[0],1.2) and near(x[1],1) and near(x[2],2))
    {
      y[0] = x[0] - 1.2;
      y[1] = x[1] - 1.0;
      y[2] = x[2] - 2.0;
    } else if(near(x[0],1.2) and near(x[1],1)) {
      y[0] = x[0] - 1.2;
      y[1] = x[1] - 1.0;
      y[2] = x[2];
    } else if(near(x[0],1.2) and near(x[2],2)) {
      y[0] = x[0] - 1.2;
      y[1] = x[1];
      y[2] = x[2] - 2.0;
    } else if(near(x[1],1) and near(x[2],2)) {
      y[0] = x[0];
      y[1] = x[1] - 1.0;
      y[2] = x[2] - 2.0;
    } else if(near(x[0],1.2)) {
      y[0] = x[0] - 1.2;
      y[1] = x[1];
      y[2] = x[2];
    } else if(near(x[1],1)) {
      y[0] = x[0];
      y[1] = x[1] - 1.0;
      y[2] = x[2];
    } else if(near(x[2],2)) {
      y[0] = x[0];
      y[1] = x[1];
      y[2] = x[2] - 2.0;
    } else {
      y[0] = -1000;
      y[1] = -1000;
      y[2] = -1000;
    }
  }
};

class Mu : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = 1000.0;
  }
};

int main()
{
  Point p0(0,0,0);
  Point p1(1.2,1,2);
  auto mesh = std::make_shared<BoxMesh>(p0, p1, 10,10,10);
  auto per_boundary = std::make_shared<PeriodicBoundary>();
  auto V = std::make_shared<PeriodicProblem::FunctionSpace>(mesh, per_boundary);

  auto mu = std::make_shared<Mu>();
  PeriodicProblem::BilinearForm a(V, V, mu);
  PeriodicProblem::LinearForm L(V);

  Function u(V);
  solve(a == L, u);

  plot(u);
  interactive();

  File file("u.pvd");
  file << u;

  return 0;
}​

But when I run the executable I get
*** Error: Unable to assemble form.
*** Reason: Coefficient number 0 ("mu") has not been set.
*** Where: This error was encountered inside AssemblerBase.cpp.
By looking at the PeriodicProblem.h I can clearly see that I have two ways setting up the Coefficient "mu". Giving it to a constructor or accessing the member directly by "a.mu". Neither of these approaches work.

Thank you for your time.
Community: FEniCS Project

2 Answers


3
12 months ago by
You have to set mu in your linear and bilinear forms.  Try something like this (see this example as reference):

auto mu = std::make_shared<Constant>(1.0); // or whatever what you want
a.mu = mu;
L.mu = mu;​
Thank you for your answer but this does not help. As I stated, neither option how to pass Coefficient mu to the BilinearForm does work:
auto mu = std::make_shared<Constant>(1000.0);
PeriodicProblem::BilinearForm a(V, V, mu);​
auto mu = std::make_shared<Constant>(1000.0);
PeriodicProblem::BilinearForm a(V, V);
a.mu = mu;

It seems that the problem must be in the UFL file. The only that differs from the C++ demos is the definition of the elasticity tensor and the index notation in the equation definition. But when I check the header file generated by FFC I see that the Coefficient mu gets recognized properly and is used in the definition of Form_a class.

written 12 months ago by Marek Netusil  
Sorry, but with the changes that I am proposing your code works fine for me. You just have to add (this lines doesn't appears in your main code)

a.mu = mu;
L.mu = mu;​​

 

written 12 months ago by Hernán Mella  
Thanks for the reply. I have realized where the problem is. The initialization of mu in the BilinearForm works fine, I just forgot to initialize it also for L, as the elasticity tensor appears also on the RHS. Code works fine now.
written 12 months ago by Marek Netusil  
0
12 months ago by
What about mu = Constant(tetrahedron)?
Please login to add an answer/comment or follow this question.

Similar posts:
Search »