Define two materials C tensor

1

I am trying to define two material properties but not getting correct answer. Is there anyway to define correctly? Please find the MWC.

Ex,Ey,nuxy,Gxy=100.,10.,0.3,5.0
S=as_matrix([[1./Ex,nuxy/Ex,0],[nuxy/Ex,1./Ey,0.0],[0.0,0.0,1./Gxy]])
C1=inv(S)
class K(UserExpression):def __init__(self,materials,d1,d2,**kwargs):
    super().__init__(**kwargs)
    self.materials=materials
    self.d1=d1
    self.d2=d2
    def value_shape(self):
          return (3,3)
    def eval_cell(self,values,x,cell):
          if self.materials[cell.index] == 1:
              values = self.d1 
          else:
              values = self.d2
              #print (values)
C = K(materials,C1,C1,degree=0)
2

Could you elaborate on what you imply by not Getting the right answer? As you have not supplied a fully minimal working code, it is hard to interpret your statement.

As far as I can tell, you are setting the same matrix regardless of your marking.

3 
from dolfin import *
from mshr import *

mesh = UnitSquareMesh(10,10)

tol = 1E-14

class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] <= 0.5 + tol

class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 0.5 - tol

materials = MeshFunction('size_t', mesh,mesh.topology().dim(), 0)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
materials.set_all(0)
subdomain_1.mark(materials, 1)

File("Materials.pvd")<< materials

Ex, Ey, nuxy, Gxy = 100., 10., 0.3, 5.
S = as_matrix([[1./Ex,nuxy/Ex,0.],[nuxy/Ex,1./Ey,0.],[0.,0.,1./Gxy]])
C1 = inv(S)
print (C1)
class K(UserExpression):
    def __init__(self,materials,d1,d2,**kwargs):
        super().__init__(**kwargs)
        self.materials=materials
        self.d1=d1
        self.d2=d2
    def value_shape(self):
        return (3,3)
    def eval_cell(self,values,x,cell):
        if self.materials[cell.index] == 1:
            values = self.d1 
        else:
            values = self.d2
        #print (values)
C = K(materials,C1,C1,degree=0)             
def eps(v):
    return sym(grad(v))
def strain2voigt(e):
    """e is a 2nd-order tensor, returns its Voigt vectorial representation"""
    return as_vector([e[0,0],e[1,1],2*e[0,1]])
def voigt2stress(s):
    """
    s is a stress-like vector (no 2 factor on last component)
    returns its tensorial representation
    """
    return as_tensor([[s[0], s[2]],
                     [s[2], s[1]]])
    def sigma(v):
        return voigt2stress(dot(C, strain2voigt(eps(v))))

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], L) and on_boundary
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0) and on_boundary

facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Top().mark(facets, 1)
Left().mark(facets, 2)
Bottom().mark(facets, 3)
ds = Measure('ds', subdomain_data=facets)

V = VectorFunctionSpace(mesh, 'Lagrange', 2)


du = TrialFunction(V)
u_ = TestFunction(V)
u = Function(V, name='Displacement')
a = inner(sigma(du), eps(u_))*dx


T = Constant((0, 1e-3))
l  = dot(T, u_)*ds(1)

bc = [DirichletBC(V.sub(0), Constant(0.), facets, 2),
      DirichletBC(V.sub(1), Constant(0.), facets, 3)]

solve(a == l, u, bc)

import matplotlib.pyplot as plt
p = plot(sigma(u)[1, 1]/T[1], mode='color')
plt.colorbar(p)
plt.title(r"$\sigma_{yy}$",fontsize=26)
plt.show()
4

First of all, encapsulate the code with ```
to format it properly.
Secondly, could you explain in words what is not working or going as expected.
Thirdly,
Can you reproduce the problem with a built in mesh? As you have not supplied us with a mesh, we cannot reproduce your error.

5

Mesh: I am not able to upload. I need to define two material (C1 and C2) but if I define a material (C1 and C1) using this way, It is unable to reproduce the results (for a single material, which we can get by using C1 in sigma definition).

6

I have defined mesh with builtin mesh. Please refer modified code.

1 Like
7

Note that the input to the eval_cell function in user-expression is an numpy array of shape (9,). Thus, you need to convert the input to this shape.
Additionally, you are not assigning values to values, but overwriting the variable. Below I’ve created a functioning version:

Ex, Ey, nuxy, Gxy = 100., 10., 0.3, 5.
mport numpy as np
C_numpy = np.linalg.inv([[1./Ex,nuxy/Ex,0.],[nuxy/Ex,1./Ey,0.],[0.,0.,1./Gxy]])
class K(UserExpression):
    def __init__(self,materials,d1,d2,**kwargs):
        super().__init__(**kwargs)
        self.materials=materials
        self.d1= d1.reshape(-1)
        self.d2=d2.reshape(-1)
    def value_shape(self):
        return (3,3)
    def eval_cell(self,values,x,cell):
        if self.materials[cell.index] == 1:
            values[:]= self.d1 
        else:
            values[:] = self.d2
    
C = K(materials,C_numpy,C_numpy,degree=0)
8

Thanks. It works :slight_smile: