How to define a subdomain using material data (from data file)?

General
1

I want to define a subdomain for 3 materials within my main domain. However, my material is randomly distributed within the domain. And I have a .dat file which contains the information about the material at each point within the domain. My .dat files has the information as follows:
(x , y, z, material_identity)
0,0,0,1
1,0,0,1
2,0,0,1
3,0,0,2
4,0,0,2
5,0,0,3
6,0,0,3
.
.
.
99,99,99,1

I have also mapped this data on a function and it looks something like this:
material

However, I am not able to define the subdomains based on this material data. I finally aim to solve mechanics for this domain. Can someone please help me to solve this issue?

thanks a lot.

1 Like
2

Without providing a minimal example, it is really hard to give you any guidance.
For instance, are you using legacy DOLFIN or DOLFINx?

Secondly, does your “point” coincide with the vertices of the cells? If so, how do you want a cell with more than one distinct value at is vertices to be defined?

3

Hi Dokken,

I am currently using DOLFIN.
I mistakenly wrote “cells” as “points” in my original question. The data that I have is for each cell within the domain. The first 3 columns (x,y,z) in my data file correspond to the cell index in the respective directions, and the 4th column denotes the material identity (1-> material 1, 2-> material 2, etc.).

My code goes as follows:

Step1: Define rectangular mesh (here I am trying to do a 2D mesh)
Step2: Input material data file (this contains x, y, material identity. x and y denotes cell indices in corresponding direction as explained above)
Step3: Define distinct subdomains for each material in the above-defined rectangular mesh

Please let me know if you need any more details regarding the query.

4

First of all, does your mesh consist of triangles or quadrilaterals?
Secondly, is your mesh uniform and structured (i.e. is it a RectangleMesh with 100 x 100 cells?)

I would suggest something like:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(100, 100)

V = FunctionSpace(mesh, "DG", 0)
u = Function(V)

dx = 1/100

values = np.arange(100*100)
local_values = np.zeros_like(u.vector().get_local())
for cell in cells(mesh):
    midpoint = cell.midpoint().array()
    i = int(midpoint[0] // dx)
    j = int(midpoint[1] // dx)
    pos = i + (j+1)
    local_values[cell.index()] = values[pos]
    print(midpoint, i, j , pos)

u.vector().set_local(local_values)

File("u.pvd") << u
5

Dear Dokken,

Thank you so much for your help! I have used your piece of code to assign material information to the cells in the domain and it has worked perfectly. However, I am still not able to define subdomains based on the material information. I am attaching my code below. Could you please help me with this step? Thank you so much!

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
import csv

NX_points = 500
domain_length = 200e-6

# mesh = UnitSquareMesh(100, 100)
mesh = RectangleMesh(Point(0, 0), Point(domain_length, domain_length), int(NX_points), int(NX_points))

V = FunctionSpace(mesh, "DG", 0)
u = Function(V)

dx = domain_length/500

material_array = list()

#---------------Import material file----------------#

with open ('uStr_30AM_50SE.dat','r') as f:          
    csv_reader = csv.reader (f, delimiter = "\t")


    for line in csv_reader:
        x_temp = float(line[0])
        y_temp = float(line[1])
        z_temp = float(line[2])

        c_temp = line[3]

        material_array.append(c_temp)

material_array = np.array(material_array) # Material array
print(len(material_array))
#-------------------------------------------------------#

values = np.arange(500*500)

# print(values)
local_values = np.zeros_like(u.vector().get_local())
local_values_material = np.zeros_like(u.vector().get_local())
for cell in cells(mesh):
    midpoint = cell.midpoint().array()
    i = int(midpoint[0] // dx)
    j = int(midpoint[1] // dx)
    pos = i + 500*j
    local_values[cell.index()] = values[pos]
    local_values_material[cell.index()] = material_array[pos]
    print(midpoint, i, j , pos)

# u.vector().set_local(local_values)
u.vector().set_local(local_values_material)

File("u.pvd") << u

#----------------- How to set subdomains based on the material here on? ------------------------#
6

Just add the setting values to a MeshFunction inside this loop.

7

I am very new to fenics and I apologize if this sounds silly, but I am not able to understand your suggestion. If possible, could you please explain it with a coding example?

Sincerely,
ABS

8

I was able to define subdomains as per your suggestion:

f = MeshFunction('double', mesh, 2)
local_values = np.zeros_like(u.vector().get_local())
local_values_material = np.zeros_like(u.vector().get_local())
for cell in cells(mesh):
    midpoint = cell.midpoint().array()
    i = int(midpoint[0] // dx)
    j = int(midpoint[1] // dx)
    pos = i + 500*j
    local_values[cell.index()] = values[pos]
    local_values_material[cell.index()] = material_array[pos]
    f[cell] = local_values_material[cell.index()]
    print(midpoint, i, j , pos)

Thank you so much for your help.

9

Though the above piece of code works for defining the subdomain, I am getting the following error when I solve for the linear elastic mechanics:

Traceback (most recent call last):
  File "Chemo_mechanics.py", line 280, in <module>
    solve(a == l, u, bc)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 268, in _solve_varproblem
    problem = LinearVariationalProblem(eq.lhs, eq.rhs, u, bcs,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/problem.py", line 58, in __init__
    L = Form(L, form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 82, in __init__
    self.set_cell_domains(subdomains)
TypeError: set_cell_domains(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfin.cpp.fem.Form, arg0: dolfin.cpp.mesh.MeshFunctionSizet) -> None

Invoked with: <dolfin.fem.form.Form object at 0x7fbba1084720>, <dolfin.cpp.mesh.MeshFunctionDouble object at 0x7fbb98516070>

I am using the MeshFunction of type ‘double’ in the following line of the code:

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

# print(values)
local_values = np.zeros_like(u.vector().get_local())
local_values_material = np.zeros_like(u.vector().get_local())
for cell in cells(mesh):
    midpoint = cell.midpoint().array()
    i = int(midpoint[0] // dx)
    j = int(midpoint[1] // dx)
    pos = i + 500*j
    local_values[cell.index()] = values[pos]
    local_values_material[cell.index()] = material_array[pos]
    materials[cell] = local_values_material[cell.index()]
    # f.array()[:] = local_values_material[cell.index()]
    print(midpoint, i, j , pos)

I think the error lies in using the MeshFunction type ‘double’ instead of ‘size_t’. When I use the type ‘size_t’ in the above code, it is giving me the following error:

Traceback (most recent call last):
  File "Chemo_mechanics.py", line 53, in <module>
    materials[cell] = local_values_material[cell.index()]
TypeError: __setitem__(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfin.cpp.mesh.MeshFunctionSizet, arg0: int, arg1: int) -> None
    2. (self: dolfin.cpp.mesh.MeshFunctionSizet, arg0: dolfin.cpp.mesh.MeshEntity, arg1: int) -> None

Invoked with: <dolfin.cpp.mesh.MeshFunctionSizet object at 0x7f8b582511f0>, <dolfin.cpp.mesh.Cell object at 0x7f8b58c4c870>, 1.0

Could you please help me with resolving this issue. Thank you so much for your help.

Sincerely,
ABS

10

Use an integer instead of the actual physical floating point value on the rhs here.

11

Thank you so much for your suggestion. The code is working now.

12

Dear Dokken,

Would you please guide me on how to define the internal boundaries between any 2 materials within such domain? I have successfully defined the external boundaries (boundaries of the rectangular domain), but I am not able to define the internal boundaries. I have also asked this question in my recent post: here is the link
Please let me know if you need any more information which might help you to solve this query.

Thanks a lot for your help.