How to get stress field of a subdomain created by MeshView?

I have created a new from from a deformed mesh using mesh view. I want to use the existing stress stress field for the deformed mesh. So I created a tensor functional space using a submesh and Tried to project the stress field onto the subspace. But I get error stating that the point are outside the domain. Ideals option is to allow extrapolation. But it doesn’t work either. Here is my code for reference.
I install fenics in ubuntu and I am running ubuntu in virtual box and dolfin verison ‘2019.2.0.dev0’.

sub_domains = df.MeshFunction(“size_t”, DeformedMesh, DeformedMesh.topology().dim(),0)

sub_domains.set_all(1)

domain1=LeftDom()

domain2=RightDom()

domain1.mark(sub_domains,0)

domain2.mark(sub_domains,1)

leftPart=df.MeshView.create(sub_domains,0)

sig_subspace=df.SubMesh(DeformedMesh,sub_domains,0)

V1 = df.VectorFunctionSpace(leftPart, “Lagrange”, degree=2)

V1sig = df.TensorFunctionSpace(sig_subspace, ‘DG’,degree=2)

sig_init_after.set_allow_extrapolation(True)

sig=df.project(sig_init_after,V1sig)

sig.set_allow_extrapolation(True)

Please consider the following link, particularly regarding construction of a MWE:

You need to have a complete code, this means that we can copy-paste the code in to an editor and run it, without modifications. Currently we are missing several key definitions to be able to reproduce the problem. As @nate stated, please follow the guidelines.

1 Like

I don’t know if I can post the entire code online, As it uses a lot of functions

You should try to reduce your problem to the smallest possible problem, with as few variables as possible, that reproduces your error. See for instance: How to define a Function in a submesh based on function values present in an adjoining submesh? - #4 by dokken which reduces the original example: How to define a Function in a submesh based on function values present in an adjoining submesh? - #3 by Sabarish_Gopi from over 100 lines of code to 48.

The key to get help is to identify what particular function causes the issue, remove all code that should be called after this error, and then remove one and one variable, and check that your code still produces the error.

Please format your code using ``` such that the code is properly formatted (and has the correct indentation).

First of all, why do you have both SubMesh and MeshView in the same code?
You should only use MeshView.

But I did not use it.

even if I remove It won’t change anything.

The error comes only if I deform the mesh in line 97-100.

As I have said in the previous posts, please remove all code that is not needed to reproduce the error message. It makes it alot easier for others to help you if only the relevant code is shown.

Following is a minimal example that reproduces your error message. Here I have removed all code not relevant for reproducing the error. For instance, I have removed your solver, and simply added an expression into the appropriate space.

import numpy as np
import dolfin as df
E = df.Constant(110000.)
nu = df.Constant(0.3)
df.set_log_level(1)
mu = E / (2 * (1 + nu))
lmbda = (E * nu) / ((1 + nu) * (1 - 2 * nu))

def eps(v):
    return df.sym(df.grad(v))
def botLeftCorner(x, on_boundary):
    tol = 10**(-5)
    return abs(x[0]) < tol and abs(x[1]) <tol and on_boundary
def topleftSide(x, on_boundary):
    tol = 10**(-5)
    return abs(x[0]) < tol and on_boundary
mesh = df.UnitSquareMesh(36,36)
V = df.VectorFunctionSpace(mesh, "Lagrange",2,dim=2)
u=df.project(df.Expression(("x[0]", "x[1]"), degree= 1),df. VectorFunctionSpace(mesh, "CG", 1))
class LeftDom(df.SubDomain):
    def inside(self,x,on_boundary):
        return x[0] <= 0.5 + df.DOLFIN_EPS
class RightDom(df.SubDomain):
    def inside(self,x,on_boundary):
        return x[0] > 0.5 - df.DOLFIN_EPS
sub_domains = df.MeshFunction("size_t", mesh, mesh.topology().dim(),0)
sub_domains.set_all(1)
domain1=LeftDom()
domain2=RightDom()
domain1.mark(sub_domains,0)
DeformedMesh = df.Mesh(mesh)
X = DeformedMesh.coordinates()
X+=np.vstack(map(u,X))
submesh1 = df.MeshView.create(sub_domains, 1)
submesh2 = df.MeshView.create(sub_domains, 0)
V1 = df.VectorFunctionSpace(DeformedMesh, "Lagrange", degree=2)
V1sig = df.TensorFunctionSpace(DeformedMesh, 'Lagrange',degree=2)
V1sig_1= df.TensorFunctionSpace(submesh1,"DG",degree=1)

def sigmaE(v, lmbda, mu):
    sigE = lmbda * df.tr(eps(v)) * df.Identity(2) + 2.0 * mu * eps(v)
    return sigE
sigma=sigmaE(u,lmbda,mu)

#u.set_allow_extrapolation(True)
sig = df.project(sigma,V1sig)

Note that adding

u.set_allow_extrapolation(True)

before the projection removes the error message.

sorry, I made a mistake in pasting correct version of the code, Please check the following code:

   import dolfin as df
   E = df.Constant(110000.)
   nu = df.Constant(0.3)
   df.set_log_level(1)
   mu = E / (2 * (1 + nu))
   lmbda = (E * nu) / ((1 + nu) * (1 - 2 * nu))

   def eps(v):
       return df.sym(df.grad(v))
   def botLeftCorner(x, on_boundary):
       tol = 10**(-5)
       return abs(x[0]) < tol and abs(x[1]) <tol and on_boundary
   def topleftSide(x, on_boundary):
       tol = 10**(-5)
       return abs(x[0]) < tol and on_boundary
   mesh = df.UnitSquareMesh(36,36)
   V = df.VectorFunctionSpace(mesh, "Lagrange",2,dim=2)
   u=df.project(df.Expression(("x[0]", "x[1]"), degree= 1),df. VectorFunctionSpace(mesh, "CG", 1))
   class LeftDom(df.SubDomain):
       def inside(self,x,on_boundary):
           return x[0] <= 0.5 + df.DOLFIN_EPS
   class RightDom(df.SubDomain):
       def inside(self,x,on_boundary):
           return x[0] > 0.5 - df.DOLFIN_EPS

   DeformedMesh = df.Mesh(mesh)
   X = DeformedMesh.coordinates()
   X+=np.vstack(map(u,X))
   sub_domains = df.MeshFunction("size_t", DeformedMesh, DeformedMesh.topology().dim(),0)
   sub_domains.set_all(1)
   domain1=LeftDom()
   domain2=RightDom()
   domain1.mark(sub_domains,0)
   submesh1 = df.MeshView.create(sub_domains, 1)
   submesh2 = df.MeshView.create(sub_domains, 0)
   V1 = df.VectorFunctionSpace(DeformedMesh, "Lagrange", degree=2)
   V1sig = df.TensorFunctionSpace(DeformedMesh, 'Lagrange',degree=2)
   V1sig_1= df.TensorFunctionSpace(submesh1,"DG",degree=1)

   def sigmaE(v, lmbda, mu):
       sigE = lmbda * df.tr(eps(v)) * df.Identity(2) + 2.0 * mu * eps(v)
       return sigE
   sigma=sigmaE(u,lmbda,mu)

   u.set_allow_extrapolation(True)
   sig = df.project(sigma,V1sig_1)```

The error message you get informs you that dolfin is unable to evaluate the function at a point because the point is not inside the domain, suggesting you to use u.set_allow_extrapolation(True) before the projection as mentioned by @dokken. And you noticed that the error occurs only when you deform the mesh.

Indeed, the function u which is used in the variable sigma=sigmaE(u,lmbda,mu) that you want to project is defined over the initial mesh (and not the deformed mesh). But the submesh (submesh1) on which you want to do the projection is extracted from the deformed mesh.
So there are points in submesh1 (subset of DeformedMesh points) that are not in mesh, and u cannot be evaluated at these points. (Note : That is why the error occurs only when deforming the mesh.)

If you redefine u on the DeformedMesh, then it works without the need of extrapolation :

u=df.project(df.Expression(("x[0]", "x[1]"), degree= 1),df.VectorFunctionSpace(DeformedMesh, "CG", 1))
sigma=sigmaE(u,lmbda,mu)