Applying spatially varying material

1

Hello,

i am now solving a problem of a plate with three different kinds of materials, it’s a unit square mesh with a inner circle warped by a ring area as shown in the picture.
image

i want to apply the inner circle with material 1 and outside the outer circle area with second material. and inside the ring area is a functionally graded material towards outside the ring, and so the material properties is dependent on the distance from the center of the inner circle.

When i tried to realize it through fenics using the UserExpression class as shown below. i met some problem during the solving.

class InnerCircle1(SubDomain):
    def inside(self, x, on_boundary):
        return sqrt(pow((x[0]-0.8),2)+pow((x[1]-0.5),2)) <= 0.05 + tol

class OuterCircle1(SubDomain):
    def inside(self, x, on_boundary):
        return (sqrt(pow((x[0]-0.8),2)+pow((x[1]-0.5),2)) >= 0.05 - tol and\
               sqrt(pow((x[0]-0.8),2)+pow((x[1]-0.5),2)) <= 0.1 + tol)

materials = MeshFunction('size_t', mesh, 2)

subdomain_1 = InnerCircle1()
subdomain_2 = OuterCircle1()

materials.set_all(0)
subdomain_1.mark(materials, 1)
subdomain_2.mark(materials, 2)

plot(materials)
File('materials.xml.gz') << materials

class K(UserExpression):
    def __init__(self, materials, mesh, k_0, k_1, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        self.mesh = mesh
        self.x = SpatialCoordinate(self.mesh)
        self.k_0 = k_0
        self.k_1 = k_1
        
    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = self.k_0
        elif self.materials[cell.index] == 1:
            values[0] = self.k_1
        else:
            values[0] = (self.k_1 + (self.k_0-self.k_1)*\
                         ((sqrt(pow((self.x[0] - 0.8),2)+pow((self.x[1] - 0.5),2)) - 0.05) / 0.05))
    
E = K(materials, mesh, E1, E2, degree=0)

i encounterd the error shown as below:

ValueError: setting an array element with a sequence.

this error i believed is caused by the third kinds of material property assignments which is coordinate dependent. How can i solve this problem?

Thanks

2

I’ve attached a working code. I’ve removed your definition of SpatialCoordinate, as you send in the appropriate coordinates to eval_cell. In eval_cell only one coordinate is sent in at once, while spatial coordiante is an array.

from dolfin import *
mesh = UnitSquareMesh(100,100)
tol = 1e-6
E1, E2 = 0.5,0.3
class InnerCircle1(SubDomain):
    def inside(self, x, on_boundary):
        return sqrt(pow((x[0]-0.8),2)+pow((x[1]-0.5),2)) <= 0.05 + tol

class OuterCircle1(SubDomain):
    def inside(self, x, on_boundary):
        return (sqrt(pow((x[0]-0.8),2)+pow((x[1]-0.5),2)) >= 0.05 - tol and\
               sqrt(pow((x[0]-0.8),2)+pow((x[1]-0.5),2)) <= 0.1 + tol)

materials = MeshFunction('size_t', mesh, 2)

subdomain_1 = InnerCircle1()
subdomain_2 = OuterCircle1()

materials.set_all(0)
subdomain_1.mark(materials, 1)
subdomain_2.mark(materials, 2)

plot(materials)
File('materials.xml.gz') << materials

class K(UserExpression):
    def __init__(self, materials, mesh, k_0, k_1, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        self.mesh = mesh
        self.k_0 = k_0
        self.k_1 = k_1

    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = self.k_0
        elif self.materials[cell.index] == 1:
            values[0] = self.k_1
        else:
            values[0] = (self.k_1 + (self.k_0-self.k_1)*\
                         ((sqrt(pow((x[0] - 0.8),2)+pow((x[1] - 0.5),2)) - 0.05) / 0.05))

    def value_shape(self):
        return ()

E = K(materials, mesh, E1, E2, degree=0)
assemble(E*dx(domain=mesh))
1 Like
3

Thanks! The problem is solved now!