Plotting a rod with two different materials (different modulues young) to paraview

Legacy DOLFIN
1

hi all,

im trying to create a laminated rode that is being defined by two subdomains, each domain for each material. those materials are defer by their moduoel young. this rode is being stretched on the left side of the rode. i want to represent this rode in paraview and see in different color each material. how can i achieve that part? what i specifically need in the part of the two color plotting , the rest i got

the relevant parts of the subdomains
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.5 + DOLFIN_EPS

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

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

class Boundaryleft(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, DOLFIN_EPS)

class Boundaryright(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1, DOLFIN_EPS)

ADDED

facets = MeshFunction(“size_t”, mesh, 2)
facets.set_all(0)
Boundaryleft().mark(facets, 1)
Boundaryright().mark(facets, 2)

class K(fe.UserExpression):
def init(self, materials, k_0, k_1, **kwargs):
super().init(**kwargs)
self.materials = materials
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
    else:
        values[0] = self.k_1

kappa = K(materials, 1.0e7, 2.0e4, degree=0)

File(“Materials.pvd”)<< materials
vtkfile_subdomains = File(‘subdomains.pvd’)
vtkfile_subdomains << materials

2

Please format your code with 3x` syntax, ie

```python
# add code here

```

Please also explain what is not working with your current code, as it is unclear what is not working for you.

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


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


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



# Function space
V = VectorFunctionSpace(mesh, 'Lagrange', 2)


# Boundary conditions
class Boundaryleft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0, DOLFIN_EPS)


class Boundaryright(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1, DOLFIN_EPS)


# ADDED
facets = MeshFunction("size_t", mesh, 2)
facets.set_all(0)
Boundaryleft().mark(facets, 1)
Boundaryright().mark(facets, 2)

'''
boundaryl = Boundaryleft()
boundaryr = Boundaryright()
bc1 = DirichletBC(V, Constant((0, 0, 0)), boundaryl)
bc2 = DirichletBC(V, Constant((1, 0, 0)), boundaryr)
'''
bc1 = DirichletBC(V.sub(0), Constant(0), facets, 1)
bc2 = DirichletBC(V.sub(0), Constant(1), facets, 2)

bc = [bc1, bc2]

class K(fe.UserExpression):
    def __init__(self, materials, k_0, k_1, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        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
        else:
            values[0] = self.k_1

kappa = K(materials, 1.0e7, 2.0e4, degree=0)

File("Materials.pvd")<< materials
vtkfile_subdomains = File('subdomains.pvd')
vtkfile_subdomains << materials
4

i want to see my rode in para view, each subdomain in different color, but im not able to doit, once i upload the subdomin file to paraview, im not able to find the material file in the coloring section
seems like there is no scalar defined or somrthing like that, so iwanted to know what im missing

image
image537×673 22.3 KB

5

It is looking fine to me. What are you expecting ?

image
image1920×1080 227 KB

6

hi, how did you do it? what did you set in paraview, im not able to see it like that?

Capture2
Capture21880×998 196 KB

this is what i see

7

image
As you can see that the max limit is no different than 0, therefore in order to change the max limit, go to Rescale to Custom Data Range and set the max limit to 1 and click rescale in order to see changes.
image

8

unfortunately still not working, any idea why? maybe its how you defined the mesh?

123
1231606×992 184 KB

9 
from dolfin import *
import numpy as np
import fenics as fe
import matplotlib.pyplot as plt

length, width, height = 1.0, 0.2, 0.2  # Dimensions of the rod
num_elements = 10  # Number of elements along each dimension
E1, E2 = 1E6, 2E6  # Young's moduli for the two halves of the rod

# Create a 3D mesh for the rod
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, num_elements, num_elements)

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


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


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



# Function space
V = VectorFunctionSpace(mesh, 'Lagrange', 2)


# Boundary conditions
class Boundaryleft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0, DOLFIN_EPS)


class Boundaryright(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1, DOLFIN_EPS)


# ADDED
facets = MeshFunction("size_t", mesh, 2)
facets.set_all(0)
Boundaryleft().mark(facets, 1)
Boundaryright().mark(facets, 2)

'''
boundaryl = Boundaryleft()
boundaryr = Boundaryright()
bc1 = DirichletBC(V, Constant((0, 0, 0)), boundaryl)
bc2 = DirichletBC(V, Constant((1, 0, 0)), boundaryr)
'''
bc1 = DirichletBC(V.sub(0), Constant(0), facets, 1)
bc2 = DirichletBC(V.sub(0), Constant(1), facets, 2)

bc = [bc1, bc2]

class K(fe.UserExpression):
    def __init__(self, materials, k_0, k_1, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        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
        else:
            values[0] = self.k_1

kappa = K(materials, 1.0e7, 2.0e4, degree=0)

File("Materials.pvd")<< materials
vtkfile_subdomains = File('subdomains.pvd')
vtkfile_subdomains << materials
10

ok i noticed that i had a mistake in the subdomain, now it working
thank you so much!

11

if i want the colors of the subdomain to also be reflected in the xdmf file, how can i do it? i want to see after applying the load from the boundary how each side deformed

12

You would have to store the MeshFunction and the Function to the same file, then you can apply Filters such as Extract Block, Warp by Vector etc.