How to obtain the element-wise stiffness matrix?

1

Hi all, how to obtain the element / local / cell matrix of the following example in dolfinx? And is there a function that can directly calculate the coordinates of the center of each cell?

import dolfinx
import numpy as np
from mpi4py import MPI
from dolfinx.cpp.mesh import CellType
import ufl

mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [2, 1], CellType.quadrilateral)
V = dolfinx.VectorFunctionSpace(mesh, ("CG", 1))
lambda_ = 1
mu = 1
T = dolfinx.Constant(mesh, (0, -1)) # traction
f = dolfinx.Constant(mesh, (0, 0)) # body force
def epsilon(u):
    return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu*epsilon(u)
ds = ufl.Measure("ds", domain=mesh)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
2

As the first part of the question is a duplicate of:

I’ll just supply you with a reference on how to compute the midpoint of all cells local to a process:

1 Like
3

Hi Dokken, thanks for your reply. Do we have a way to calculate the area/volume of each cell? Or at least could we obtain the coordinate list of each cell?

4

If you are using an affine mesh, this can be easily done using the following code to obtain the coordinates per cell:

import dolfinx
from mpi4py import MPI
import numpy

mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 2, 2)

num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
g_indices = dolfinx.cpp.mesh.entities_to_geometry(
    mesh, mesh.topology.dim, numpy.arange(num_cells, dtype=numpy.int32), False)
for c in range(num_cells):
    print(c, mesh.geometry.x[g_indices[c]])

For non-affine cases, it becomes more complicated, and a simple way of doing this for small problems is:

mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 2, 2)

num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
mts = [dolfinx.MeshTags(mesh, mesh.topology.dim, numpy.array([c], dtype=numpy.int32),
                        numpy.array([c], dtype=numpy.int32)) for c in range(num_cells)]
forms = [dolfinx.Constant(mesh, 1) * ufl.dx(domain=mesh,
                                            subdomain_data=mts[c],
                                            subdomain_id=c)
         for c in range(num_cells)]
volumes = [dolfinx.fem.assemble_scalar(forms[c]) for c in range(num_cells)]
print(volumes)

However, note that the last approach is not very efficient for large number of cells, and should be implemented in C++ or using numba.

2 Likes
5

Hi, I am trying to apply a technique called the decimation technique to my stiffness matrix in a 2D plane strain problem. Let’s say the domain is a simple rectangle that has a homogeneous Dirichlet BC on its left edge and homogeneous Neumann BC on other edges. This method needs to consider the whole domain (here the rectangle) as connected sub-elements (small rectangles connected to each other). More specifically, if we consider that the whole rectangle is made of 20 thin rectangles, to apply the decimation technique I need to know what are those parts of the stiffness matrix corresponding to each of the thin rectangles. Certainly, in addition to those blocks in the stiffness matrix which are related to the thin rectangles, there should be some blocks which define the connection between each pair of rectangles which are next to each other. Is there any way I can do it in Fenics?