Define different materials via the sub domain

Mirialex · June 13, 2024, 12:58pm

hi all,

im trying to set different material properties on my model and for some reason in the view of paraview it does not seems to be working, i dont see my model colored in 2 different colors
when i tried different subdomain its working as expected
for exp

import os
import fenics as fe
from dolfin import *
import matplotlib.pyplot as plt
from ffc.fiatinterface import create_quadrature as cquad
import numpy as np
import meshio
import dolfin

#rest of code
class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= 0.5 + DOLFIN_EPS

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

but when i try this subdomain config its not working

E1, E2, nu = 30e6, 3e6, 0.49

# Create mesh and define function space

length, width, height = 1.0, 0.2, 0.2  # Dimensions of the rod
num_elements = 40  # Number of elements along each dimension
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, 4, 4)

V = VectorFunctionSpace(mesh, "Lagrange", 1)

#creating subtomains and boundaty conditions
tol = 1E-14

class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return np.sqrt(((x[1] - 0.1)*(x[1] - 0.1)) + ((x[2] - 0.1)*(x[2] - 0.1))) <= 0.05

class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return np.sqrt(((x[1] - 0.1)*(x[1] - 0.1)) + ((x[2] - 0.1)*(x[2] - 0.1))) > 0.05

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

boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

left = AutoSubDomain(lambda x: near(x[0], 0))
right = AutoSubDomain(lambda x: near(x[0], 1))
up = AutoSubDomain(lambda x: near(x[2], 0.2))
down = AutoSubDomain(lambda x: near(x[2], 0))
front = AutoSubDomain(lambda x: near(x[1], 0))
back = AutoSubDomain(lambda x: near(x[1], 0.2))

left.mark(boundary_parts, 1)
right.mark(boundary_parts, 2)
up.mark(boundary_parts, 5)
down.mark(boundary_parts, 6)
front.mark(boundary_parts, 3)
back.mark(boundary_parts, 4)

bc1 = DirichletBC(V.sub(0), Constant(0), boundary_parts, 1)
bc3 = DirichletBC(V.sub(1), Constant(0), boundary_parts, 3)
bc5 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 5)
bc6 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 6)
tDirBC = Expression(('1.0*time_'), time_=0.0, degree=0)
tDirBC4 = Expression(('((1/(time_+1))-1)*0.2'), time_=0.0, degree=0)
bc2 = DirichletBC(V.sub(0), tDirBC, boundary_parts, 2)
bc4 = DirichletBC(V.sub(1), tDirBC4, boundary_parts, 4)

bcs = [bc1, bc2, bc3, bc4, bc5, bc6]
subdomain_ids = [2, 4]
ds = Measure("ds", subdomain_data=boundary_parts, subdomain_id=(2,4))
ds = ds(degree=4)
#creating material constants
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] == 1:
            values[0] = self.k_0
        else:
            values[0] = self.k_1

    def value_shape(self):
        return ()


E_global = K(materials, E1, E2, degree=0)

DG = FunctionSpace(mesh, "DG", 2)
materials_function = Function(DG)
E_values = [E1, E2]
materials_function = project(E_global, DG)

MU1, LAMBDA1 = Constant(E1 / (2 * (1 + nu))), Constant(E1 * nu / ((1 + nu) * (1 - 2 * nu)))
MU2, LAMBDA2 = Constant(E2 / (2 * (1 + nu))), Constant(E2 * nu / ((1 + nu) * (1 - 2 * nu)))


class MUMU(fe.UserExpression):
    def __init__(self, materials, MU_0, MU_1, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        self.MU_0 = MU_0
        self.MU_1 = MU_1

    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 1:
            values[0] = self.MU_0
        else:
            values[0] = self.MU_1

    def value_shape(self):
        return ()


MU = MUMU(materials, MU1, MU2, degree=0)


class lamabada(fe.UserExpression):
    def __init__(self, materials, lam_0, lam_1, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        self.lam_0 = lam_0
        self.lam_1 = lam_1

    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 1:
            values[0] = self.lam_0
        else:
            values[0] = self.lam_1

    def value_shape(self):
        return ()

LAMBDA = lamabada(materials, LAMBDA1, LAMBDA2, degree=0)

# rest of the code
File("Materials.pvd") << materials
vtkfile_subdomains = File('subdomains.pvd')
vtkfile_subdomains << materials

#writing into files in order to check it in paraview
materials_function.rename("materials", "")
u.rename("Displacement Vector", "")
beam_deflection_file = fe.XDMFFile("beam_deflection.xdmf")
beam_deflection_file.parameters["flush_output"] = True
beam_deflection_file.parameters["functions_share_mesh"] = True
beam_deflection_file.write(u, 0.0)
beam_deflection_file.write(materials_function, 0.0)

would appreciate to know if you have any idea why

this is what i see

francesco-ballarin · June 13, 2024, 4:16pm

My understanding is that you are trying to define a circle/sphere around (x,y) = (0.1, 0.1). You may need to cut the visualization in paraview with a clip to see that.

I’ve already told you in the past to make sure to include all imports in your snippet, so that people who answer you can easily do a copy and paste. Last time I fixed that for you, this time I won’t, hence I won’t be trying to run your code (because it is not runnable). Please try to remember to do this: it’s a simple task, but you must put the others in the best position possible to help you.

Mirialex · June 13, 2024, 4:29pm

added the imports, sorry for that

i tried the cute and as you can see still no change

francesco-ballarin · June 14, 2024, 11:28am

It seems to me that you are trying to identify a circle in the cross section of the beam, but the mesh of each cross section is so coarse that you’ll never identify the circle.

The following code

from dolfin import *
import numpy as np

length, width, height = 1.0, 0.2, 0.2  # Dimensions of the rod
num_elements = 40  # Number of elements along each dimension
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, num_elements, num_elements)

class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return np.sqrt(((x[1] - 0.1)*(x[1] - 0.1)) + ((x[2] - 0.1)*(x[2] - 0.1))) <= 0.05

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

# rest of the code
File("Materials.pvd") << materials

does produce two subdomains. Note that:

I had to increase the number of elements in the cross-section
there is no need to have a SubDomain for \Omega_0, and then negate that expression to identify \Omega_1. Simply set the default value to 2 to say that all cells initially belong to \Omega_1, and then identify \Omega_0.

Mirialex · June 14, 2024, 3:54pm

hi thanks for you reply and time checking it
but when run it this is the error i got

UMFPACK V5.7.9 (Oct 20, 2019): ERROR: out of memory

Traceback (most recent call last):
  File "/home/mirialex/PycharmProjects/fibers/13-06-24-take2.py", line 390, in <module>
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     https://fenicsproject.discourse.group/
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside ./dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.13.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------


Process finished with exit code 1

im using a solver once without time dependency and once with
without the time dependency im able to get resolution but once i getting to the loop of the time dependency it fails with the error above
i tried to lower the load steps in the loop or the number of the elements but nothing helped.

if you have any other piece of advice, would be more then happy to read

import fenics as fe
from dolfin import *
import matplotlib.pyplot as plt
from ffc.fiatinterface import create_quadrature as cquad
import numpy as np

# rest of the code

du = TrialFunction(V)  # Incremental displacement
v = TestFunction(V)  # Test function
u = Function(V)  # Displacement from previous iteration
B = Constant((0.0, 0.0, 0.0))  # Body force per unit volume
T = Constant((0.0, 0.0, 0.0))  # Traction force on the boundary

# Project or interpolate the initial guess onto u
initial_guess = Expression(("x[0]", "x[1]", "x[2]"), degree=1)
# Project or interpolate the initial guess onto u
u.interpolate(initial_guess)
# u.interpolate(Constant((0,0,0)))

# Kinematics
d = len(u)
I = Identity(d)  # Identity tensor
F = variable(I + grad(u))  # Deformation gradient
C = F.T * F  # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J = det(F)

# Stored strain energy density (compressible neo-Hookean model)

W = MU / 2 * (tr(C) - 3 - 2 * ln(J)) + LAMBDA / 2 * pow((J - 1), 2)

n = FacetNormal(mesh)
P = diff(W, F)
# Total potential energy
L = inner(P, grad(v)) * dx - dot(B, v) * dx - dot(T, v) * ds

Je = derivative(L, u, du)

problem = NonlinearVariationalProblem(L, u, bcs=bcs,J=Je)

solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'lu'

#rest of the code
dt = 0.0001
t, T = 0.0, 1 * dt
TT = Traction()

while t <= T:
    print('time: ', t)

    # Increase traction
    TT.t = t
    tDirBC.time_ = t
    tDirBC4.time_ = t
    #rest of the code
    solver.solve()
    t += float(dt)

Mirialex · June 14, 2024, 4:25pm

when i lower it back to number of elemnts of the number 4
its working but i still dont see the 2 materials

francesco-ballarin · June 15, 2024, 4:20pm

Try using mumps instead of umfpack. If you don’t know how to do that, it would be appreciated that you made a search (useful keywords could be NonlinearVariationalSolver and mumps) before asking us how to do that.

its working but i still dont see the 2 materials

As I told you, 4 elements in the cross section is too coarse to obtain a decent representation of the circle.
If you want to use a very coarse mesh, the suggestion would be to have it generated by gmsh (you know how to do that, since we discussed it last week), so that the subdomains are build during the mesh generation process, rather than afterwards.

Topic		Replies	Views
Subdomains with different material properties	25	4884	April 22, 2024
Defining and plotting subdomains for different materials	3	667	May 31, 2020
Problem when defining different materials General	1	200	November 8, 2023
Plotting a rod with two different materials (different modulues young) to paraview Legacy DOLFIN visualization	11	223	December 29, 2023
Hyperelastic model with two different materials General	1	215	February 4, 2024

Define different materials via the sub domain

Related topics