Subject: [FEniCS-users] function spaces

From: David Beacham <[email-masked]>
Date: Mon, 5 Oct 2009 22:41:25 +0100

Hi all,

I'm playing around with using ALE methods for free surface fluids, and 
have solved for the initial free surface displacement, u, and pressure, 
p. I then want to displace the mesh using a linear elastic FEM soln, 
with the boundary displacement given by u on the boundary.

However, when I use the following code:

    linearElasticity::FunctionSpace E(mesh);
    linearElasticity::BilinearForm a1(E,E);
    linearElasticity::LinearForm L1(E);

    Function uE(E); uE = u;   //problem originates here?
    Constant lmbda(1,1.0);
    Constant mu(1,1.0);
    Constant f1(2,0.0);

    std::vector<BoundaryCondition*> bcs1;
    DirichletBC bc11(E, uE, extFacetDomains, 1);
    bcs1.push_back(&bc11);
    L1.f = f1;
    a1.mu = mu; a1.lmbda = lmbda;
    VariationalProblem pde(a1, L1, bcs1);

    Function e(E);
    pde.solve(e);

I get the error :
> terminate called after throwing an instance of 'std::runtime_error'
>   what():  *** Error: Unable to evaluate function at given point (not 
> inside domain).

However, if I try
    linearElasticity::FunctionSpace E(mesh);
    linearElasticity::BilinearForm a1(E,E);
    linearElasticity::LinearForm L1(E);

    Squish uE(E);
    Constant lmbda(1,1.0);
    Constant mu(1,1.0);
    Constant f1(2,0.0);

    std::vector<BoundaryCondition*> bcs1;
    DirichletBC bc11(E, uE, extFacetDomains, 1);
    bcs1.push_back(&bc11);
    L1.f = f1;
    a1.mu = mu; a1.lmbda = lmbda;
    VariationalProblem pde(a1, L1, bcs1);

    Function e(E);
    pde.solve(e);

where Squish is some given Function.

I assume the problem is the way I'm trying to move the function u (on 
NavierStokes::FunctionSpace V(mesh)) to uE (on 
linearElasticity::FunctionSpace E(mesh)), but I can't seem to fix it 
using interpolate() etc.

Thanks for any help.
David

---

From: Anders Logg <[email-masked]>
Date: Mon, 5 Oct 2009 23:53:14 +0200

Looks like the problem is that you are using functions defined on two
separate meshes. u is defined on the structure mesh and uE is defined
on the fluid mesh.

Try looking at the SubMesh class which can be used to extract two
matching meshes (fluid and structure) from a common global mesh. The
two meshes that you extract will contain mappings from between the
matching vertices in the two meshes so that you can transfer the
values from one mesh to the other on the boundary.

--
Anders

On Mon, Oct 05, 2009 at 10:41:25PM +0100, David Beacham wrote:
> _______________________________________________
> FEniCS-users mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/fenics-users

---

From: David Beacham <[email-masked]>
Date: Tue, 6 Oct 2009 14:05:52 +0100

Thanks for the response - still a bit confused about this though as I 
only have 1 mesh - I'm considering a freely oscillating droplet.

The free surface is defined as subdomain ds(1) in the ufl form for the 
Navier Stokes and passed through the MeshFunction extFacetDomains. The 
ale section requires the definition of an essential boundary condition u 
= uE on the boundary making use of the same MeshFunction. In this case 
would I still need to use SubMesh?

The condensed code (without using SubMesh) is below:

  Mesh mesh("../meshes/ellipse.xml");
  mesh.order();

  // Apply subdomains
  MeshFunction<unsigned int> cellDomains(mesh, mesh.topology().dim());
  cellDomains = 0;

  MeshFunction<unsigned int> intFacetDomains(mesh, mesh.topology().dim() 
- 1);
  intFacetDomains = 0;

  MeshFunction<unsigned int> extFacetDomains(mesh, mesh.topology().dim() 
- 1);
  extFacetDomains = 0; FreeSurfaceSD freeSurfaceSD; 
freeSurfaceSD.mark(extFacetDomains, 1);

  // Boundary conditions
  std::vector<BoundaryCondition*> bcs;

  // Navier-Stokes equations
  unsteadyNavierStokes::FunctionSpace V(mesh);

  unsteadyNavierStokes::BilinearForm a(V, V);
  unsteadyNavierStokes::LinearForm L(V);
  VariationalProblem pde(a, L, bcs, &cellDomains, &extFacetDomains, 
&intFacetDomains, true);

... //solve for u then consider mesh motion ...

  linearElasticity::FunctionSpace E(mesh);
  linearElasticity::BilinearForm ale_a(E,E);
  linearElasticity::LinearForm ale_L(E);

  Function uE(E); uE = u0;

  std::vector<BoundaryCondition*> ale_bcs;
  DirichletBC ale_bc1(E, uE, extFacetDomains, 1);
  bcs1.push_back(&ale_bc1);
  VariationalProblem ale_pde(ale_a, ale_L, ale_bcs);

which gives the 'not inside domain' error below.

In this instance, why would uE not be defined in the domain?

Thanks, David

---

From: Anders Logg <[email-masked]>
Date: Mon, 12 Oct 2009 19:13:21 +0200

> would I still need to use SubMesh?

No, I thought you had two domains but it looks like you just have one.

Install the latest DOLFIN (0.9.4) and try adding print some debug
information inside

  void Function::eval()

starting on line 243 in dolfin/function/Function.cpp and see if you
can figure out what goes wrong.

It could be a round-off error that prevents DOLFIN from finding the
correct point when you transfer the function from one mesh to the other.

--
Anders