imposing "measured data" to Dirichlet boundary conditions

Question

I'm relatively new to fenics and I just looked through all questions related to Dirichlet boundary conditions. I don't seem to find a well-described question or answer about what I'm about to ask.

I'm solving a PDE-constrained optimization problem and I have "measured data" from experiment defined on all nodes of my mesh. I want to impose part of the "measured data" to my Dirichlet boundary condition and impose 0 Neumann boundary condition to the nodes where the "measured data" is not specified.

After some reading, I think I can read in this "measured data" together with mesh as meshdata or mesh_function, since it is defined on mesh nodes. It seems possible, though I haven't figured out how to do this exactly.

But even with that "measured data" loaded in, I still don't know how to impose them on part of boundary nodes as Dirichlet boundary conditions and impose 0 Neumann bc for the rest of the boundary nodes. Applying them pointwisely seems to be very inefficient. Tagging subdomains doesn't seem to apply to this situation very well.

For example, in my case, I 'd like to have an array, measuredData that is of the same size as the coordinates x, so that I can do something like:

mesh = UnitSquareMesh(3,2)
V = VectorFunctionSpace(mesh,"CG",1)
measD0 = Expression("measuredData[0]")
measD1 = Expression("measuredData[1]")
bc0 = DirichletBC(V.sub(0), measD0, on_boundary)
bc1 = DirichletBC(V.sub(1), measD1, on_boundary)

But now I haven't figured out a way to overload Expression to make it see my measuredData.

Keen to hear experienced fenics' users' opinions! Thanks in advance!

Answer 1

Constant((0, 3, 2, 4)) is the vector \(\begin{pmatrix} 0 \\ 3 \\ 2 \\ 4 \end{pmatrix} \) and V is a scalar function space, so of course you cannot assign a vector if you only have degrees of freedom for a scalar. The error is pretty clear about this.

This, however, does not seem to match your original question.

Comments

Thanks for your heads-up. I see what you mean. It's true what I tried didn't apply to my original question. I intended to give the 4-component array to my 4 nodes in my mesh. Now I understand that can't be done like this. I think the Expression or Constant passed to DrichletBC is evaluated at nodal level and I may need to somehow overload the Expression operator or add some features into it.

---

There is a standard way to do this, which is

class Source : 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);
  }
};

as taken from https://fenicsproject.org/documentation/dolfin/1.6.0/cpp/demo/documented/poisson/cpp/documentation.html .

---

Thanks for bringing this up. I'm actually aware of this. But this, again, can only modify the coordinates in a constant way, meaning that the expression is constant for all nodes. Naively, I can iterate through all my points and do something like that, but that's simply too expensive and not really practically feasible.

In my case, I'd like to operate on "measured data" the same way as it operates on the coordinates x. Now the reason that fenics can operate on x is because x is available as an object in Expression. So I need to load my "measured data" also as an existing object in Expression, but I don't see this possible as far as I've read.

---

I am sorry, I am still not sure I get the question. Perhaps you are looking for this:

class ExprWithData : public Expression
{
  private:
  T* myData;

  public:
  ExprWithData(const T* _myData) : Expression(), myData(_myData);
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = myData->myMethod(...);
  }
};

If this is again not what you are asking, apologies. I'll read it again tomorrow and try to better interpret!

---

That's ok. Thanks for your quick reply. I think that your code is a user-defined JIT-compiled Expression, similar to this, right? In my opinion, this seems to be a different approach but doing the same thing comparing to the approach in your previous comment. In your code, myData can only be a scalar value and it still remains the same for all points, am I understanding your reply correctly?

In my case, I 'd like to have an array, measuredData that is of the same size as the coordinates x, so that I can do something like:

mesh = UnitSquareMesh(3,2)
V = VectorFunctionSpace(mesh,"CG",1)
measD0 = Expression("measuredData[0]")
measD1 = Expression("measuredData[1]")
bc0 = DirichletBC(V.sub(0), measD0, on_boundary)
bc1 = DirichletBC(V.sub(1), measD1, on_boundary)

I just deleted my "dumb code" in my original question and added the above code. Hopefully this can be clearer.

---

Yes, now it's much clearer, although my previous replies look stupid at this point :D
Also, I didn't realize you are using the Python interface, my code works fine with the C++ interface.
In any case, in the snippet you wrote it looks like the boundary condition is again a constant, since measuredData is an array [of floats, I presume]; I presume you instead want the content of Expression to be an array itself. In this case, it means you want to apply a CellFunction to the boundary condition, is this the correct statement of the question?

---

I just searched and found that CellFunction is equivalent to MeshFunction, right? You are right in a sense that I do think about importing my measured data as a MeshFunction but haven't found how. Even if I load the data as MeshFunction, I still don't see how to use it as what I suggested... Do you have any idea?

---

CellFunction is a template specialization of MeshFunction for cells.

I see a major problem here: Expression::eval works with the coordinates of a point, not with the single cells. What happens if your expression is evaluated at a vertex that is shared by two or more cells? Which cell value should be used? Expression is for analytically defined functions, not for "cellwise" functions, so perhaps it is not the right object.

You can define a DirichletBC using any GenericFunction, thus including Functions. You can try putting your data into a Function defined on a zeroth order Discontinuous Galerkin finite element space, which behaves essentially as a CellFunction, and use that to declare the DirichletBC.

---

Yes, it's true that Expression might not be the right object. Because I initially thought I could define functions not based on coordinates but based on my measured data, since it is also defined on points like coordinates.

Your suggestion about defining a Function seems interesting. Could you write a few lines of fenics code to expand a little bit? Thanks!

---

A DG function has a degree of freedom for every cell of a mesh. Since the basis functions of the DG space are 1 on the cell and 0 otherwise, you should assign your measured data for cell j to the j-th degree of freedom of the DG function.
I'm assuming that the indexings match this way and that everything goes fine, but you could investigate that.