### C++ interface - Coefficient has not been set.

182
views
0
7 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.

Community: FEniCS Project

3
7 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 7 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 7 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 7 months ago by Marek Netusil
0
7 months ago by
What about mu = Constant(tetrahedron)?