Create function living in subspace of product of spaces

Question

I am working with a mixed formulation problem, where in my .ufl file I defined

V = FiniteElement("Lagrange", triangle, 1)
W = V*V

In the .cpp file I obtain a solution $p$ living in the dual space W which components I can plot without concern

Function p(W);
solve(a == L, p, bcs);  

Function& p1 = p[0];
Function& p2 = p[1];

plot(p1,"Real part");
plot(p2,"Imaginary part");
interactive();

I would now like to plot the sources as well, so I tried doing

Function f1(W->sub(0));
Function f2(W->sub(1));

Source1 source1;
Source2 source2;

f1.interpolate(source1);
f2.interpolate(source2);

But I get the error

Error:   Unable to create function.
Reason:  Cannot be created from subspace. Consider collapsing the function space.

So, how do I access (or define) a subspace of my dual space?

---

As a small (not) working example consider

# Compile this form with FFC: ffc -l dolfin Complex.ufl.

V = FiniteElement("CG", triangle, 1)
W = V*V

(u_r, u_i)      = TrialFunction(W)
(v_r, v_i)      = TestFunction(W)
f_r                     = Coefficient(V)
f_i                     = Coefficient(V)

a_r = inner(grad(u_r),grad(v_r))*dx - inner(grad(u_i),grad(v_i))*dx
a_i = inner(grad(u_r),grad(v_i))*dx + inner(grad(u_i),grad(v_r))*dx

L_r = inner(f_r,v_r)*dx - inner(f_i,v_i)*dx
L_i = inner(f_r, v_i)*dx + inner(f_i,v_r)*dx

a = a_r + a_i
L = L_r + L_i

---

#include <dolfin.h>
#include "Complex.h"
#include <mshr.h>

using namespace dolfin;

// Source term
class SourceReal : public Expression
{
    public:
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = cos(2*DOLFIN_PI * x[0]);
  }
};

class SourceImag : public Expression
{
    public:
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = sin(2*DOLFIN_PI * x[0]);
  }
};

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

// Sub domains for Dirichlet boundary conditions
class DirichletBCBox : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return ((x[0] < -1 + sqrt(DOLFIN_EPS) or x[0] > 1 - sqrt(DOLFIN_EPS) or x[1] < 1 + sqrt(DOLFIN_EPS) or x[1] > 1 - sqrt(DOLFIN_EPS)) and on_boundary);
  }
};


int main()
{

    //File file("results/Round_no_cloak.pvd");

    // Mesh definition
    auto rectangle = std::make_shared<mshr::Rectangle>(dolfin::Point(-1.0, 1.0, 0.0), dolfin::Point(1.0, -1.0, 0.0));
    auto mesh = mshr::generate_mesh(*rectangle, 100);

    auto W = std::make_shared<Complex::FunctionSpace>(mesh);

    auto f_real = std::make_shared<SourceReal>();
    auto f_imag = std::make_shared<SourceImag>();

    auto boundary = std::make_shared<DirichletBCBox>();

    DirichletBC bc1(W->sub(0), std::make_shared<BoundaryValue>(), boundary);
    DirichletBC bc2(W->sub(1), std::make_shared<BoundaryValue>(), boundary);
    std::vector<const DirichletBC*> bcs;
    bcs.push_back(&bc1);
    bcs.push_back(&bc2);

    Complex::BilinearForm a(W, W);
    Complex::LinearForm L(W);

    L.f_r = f_real;
    L.f_i = f_imag;

  Function w(W);
  solve(a == L, w, bcs);

  // Extract sub functions (function views)
  Function& u_real = w[0];
  Function& u_imag = w[1];

    plot(u_real,"Real part");
  plot(u_imag,"Imaginary part");

  Function p_real(W->sub(0));
    Function p_imag(W->sub(1)); 

    SourceReal source1;
    SourceImag source2;

    p_real.interpolate(source1);
    p_imag.interpolate(source2);

    plot(p_real);
    interactive();

  return 0;
}

Comments

In the end I decided to create a Test.ufl file containing the definition of V plus some random terms a and L and call it in my .cpp file.

auto W = std::make_shared<Cloaking::FunctionSpace>(mesh);
auto V = std::make_shared<Test::FunctionSpace>(mesh);

So I can now plot my function easily

Function p_real(V);
Function p_imag(V);

Source1 source1;
Source2 source2;

p_real.interpolate(source1);
p_imag.interpolate(source2);

---

However, I now want to examine the sum of my previous function, i.e.

Function q_real(V);
Function q_imag(V);

q_real = u_real + p_real;
q_imag = u_imag + p_imag;

But this generates the error

*** Error:   Unable to Construct FunctionAXPY.
*** Reason:  Expected Functions to be in the same FunctionSpace.
*** Where:   This error was encountered inside FunctionAXPY.cpp.
*** Process: 0

since FEniCS does not seem to recognize my functions to be in the same space, although they ARE.

Answer 1

Reason:  Cannot be created from subspace. Consider collapsing the function space.

As the error says, you need to "collapse" the subspace. Try

Function f1(W->sub(0)->collapse());

Comments

Thank you for your answer!

I thought collapsing meant some obscure mathematical idea, I didn't even check if there was a method with that name.

Anyway, I tried implementing your suggestion and I can now create function in the subspaces of W. However, new problems arise when I try to manipulate these functions. For instance, I'm interested in plotting the sum of the solution and the source, i.e.

    Function w(W);
    solve(a == L, w, bcs);

    Function& w_r = w[0];
    Function& w_i = w[1];       

    SourceReal source1;
    SourceImag source2;

    Function f_r(W->sub(0)->collapse());
    Function f_i(W->sub(1)->collapse());
    Function sum_r(W->sub(0)->collapse());
    Function sum_i(W->sub(1)->collapse());
    f_r.interpolate(source1);
    f_i.interpolate(source2);

    sum_r = w_r + f_r;
    sum_i = w_i + f_i;

        plot(w_r, "Solution real");
        plot(w_i, "Solution imag");
        plot(f_r, "Source real");
        plot(f_i, "Source imag");
        plot(sum_r, "Sum real");
        plot(sum_i, "Sum imag");

    interactive();

The error I get seems to suggest FEniCS doesn't understand I'm dealing with the same space over and over :/

*** Error:   Unable to Construct FunctionAXPY.
*** Reason:  Expected Functions to be in the same FunctionSpace.
*** Where:   This error was encountered inside FunctionAXPY.cpp.
*** Process: 0

---

You can't just add two functions, but you can add their attached vectors via "axpy"

*sum_r.vector() = *w_r.vector();
*sum_r.vector() += *f_r.vector();
# same for sum_i

This should work. Maybe there are better solutions, I don't really use c++!

---

Again, thanks a lot.

Seems like a good idea, but for now I have another error

*** Error:   Unable to access vector of degrees of freedom.
*** Reason:  Cannot access a non-const vector from a subfunction.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 0

I'll look into it more in detail again tomorrow.

---

You were right, it works now!
I forgot a pointer somewhere which was invalidating the whole thing.

Thank you so much for your help!

For anyone who cares, here's complete code:

#include <dolfin.h>
#include "Complex.h"
#include <mshr.h>

using namespace dolfin;

// Source term
class SourceReal : public Expression
{
    public:
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = cos(2*DOLFIN_PI * x[0]);
  }
};

class SourceImag : public Expression
{
    public:
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = sin(2*DOLFIN_PI * x[0]);
  }
};

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

// Sub domains for Dirichlet boundary conditions
class DirichletBCBox : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return ((x[0] < -1 + sqrt(DOLFIN_EPS) or x[0] > 1 - sqrt(DOLFIN_EPS) or x[1] < 1 + sqrt(DOLFIN_EPS) or x[1] > 1 - sqrt(DOLFIN_EPS)) and on_boundary);
  }
};


int main()
{

    //File file("results/Round_no_cloak.pvd");

    // Mesh definition
    auto rectangle = std::make_shared<mshr::Rectangle>(dolfin::Point(-1.0, 1.0, 0.0), dolfin::Point(1.0, -1.0, 0.0));
    auto mesh = mshr::generate_mesh(*rectangle, 100);

    auto W = std::make_shared<Complex::FunctionSpace>(mesh);
    //auto V = W->sub(0);                                                                                           //This won't work 
    //auto V = std::make_shared<Test::FunctionSpace>(mesh);

    auto f_real = std::make_shared<SourceReal>();
    auto f_imag = std::make_shared<SourceImag>();

    auto boundary = std::make_shared<DirichletBCBox>();

    DirichletBC bc1(W->sub(0), std::make_shared<BoundaryValue>(), boundary);
    DirichletBC bc2(W->sub(1), std::make_shared<BoundaryValue>(), boundary);
    std::vector<const DirichletBC*> bcs;
    bcs.push_back(&bc1);
    bcs.push_back(&bc2);

    Complex::BilinearForm a(W, W);
    Complex::LinearForm L(W);

    L.f_r = f_real;
    L.f_i = f_imag;

        Function w(W);
        solve(a == L, w, bcs);

        Function w_r = w[0];
        Function w_i = w[1];        

    SourceReal source1;
    SourceImag source2;

    Function f_r(W->sub(0)->collapse());
    Function f_i(W->sub(1)->collapse());
    Function sum_r(W->sub(0)->collapse());
    Function sum_i(W->sub(1)->collapse());
    f_r.interpolate(source1);
    f_i.interpolate(source2);

        *sum_r.vector() = *w_r.vector();
        *sum_r.vector() += *f_r.vector();
        *sum_i.vector() = *w_i.vector();
        *sum_i.vector() += *f_i.vector();

        plot(w_r, "Solution real");
        plot(w_i, "Solution imag");
        plot(f_r, "Source real");
        plot(f_i, "Source imag");
        plot(sum_r, "Sum real");
        plot(sum_i, "Sum imag");

    interactive();

  return 0;
}

---

you're welcome, glad it worked!