How to update the initial mesh?

General
1

Hi

I am trying to generate a deformed 2D mesh by displacing the initial mesh using my solution (uh). However, I am getting an incorrect outcome. Could you please check if this operation is correct?
mesh_d.geometry.x[:,:mesh_d.geometry.dim] = uh.x.array.reshape((-1, mesh_d.geometry.dim))

I would also want to keep both the initial mesh and deformed mesh to check the difference. Could you please help me how to keep the initial after updating it to deformed mesh?
How can I check the difference, since the operation I tested is not valid for mesh?
L2_error = fem.form(ufl.inner(mesh - mesh_d, mesh - mesh_d) * ufl.dx)

Thanks!

Here is the code:


from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, locate_dofs_geometrical,
                         locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary, GhostMode
from ufl import Identity, Measure, TestFunction, TrialFunction, VectorElement, dot, dx, inner, grad, nabla_div, sym
from mpi4py import MPI
from petsc4py.PETSc import ScalarType

import numpy as np
import pyvista
from dolfinx import mesh



lambda_ = 1.25
mu = 1

H=1
L=1



mesh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (H, L)), n=(10, 10),
                            cell_type=CellType.triangle,
                            ghost_mode=GhostMode.shared_facet)

element = VectorElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, element)

from dolfinx.mesh import locate_entities, meshtags


def right(x):
    return  np.isclose(x[1], 0)
boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim-1, right)
boundary_dofs_x = locate_dofs_topological(V.sub(1), mesh.topology.dim-1, boundary_facets)

bcx = dirichletbc(ScalarType(0), boundary_dofs_x, V.sub(1))
bcs = [bcx]



tol=0.000000001
boundaries_n = [(1, lambda x: np.isclose(x[1], 0,tol)),
                (2, lambda x: np.isclose(x[1], H,tol)),
              (3, lambda x: np.isclose(x[0], 0,tol)),
              (4, lambda x: np.isclose(x[0], L,tol))]

g2=Constant(mesh, ScalarType((0,0.1)))
g3=Constant(mesh, ScalarType((0,0)))
g4=Constant(mesh, ScalarType((0,0)))


facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries_n:
    facets = locate_entities(mesh, fdim, locator)
    print(facets)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])


from ufl import Measure
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

def epsilon(u):
    return sym(grad(u)) 
def sigma(u):
    return lambda_ * nabla_div(u) * Identity(len(u)) + 2*mu*epsilon(u)

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(mesh, ScalarType((0, 0)))
a = inner(sigma(u), epsilon(v)) * dx
L = dot(f, v) * dx +inner(g2, v) * ds(2)+inner(g3, v) * ds(3)+inner(g4, v) * ds(4)

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

mesh_d=mesh


# calculating the deformed mesh by adding the displacement (uh) to the initial mesh
mesh_d.geometry.x[:,:mesh_d.geometry.dim] = uh.x.array.reshape((-1, mesh_d.geometry.dim))

#Difference between the mesh
from dolfinx import fem
import ufl
L2_error = fem.form(ufl.inner(mesh - mesh_d, mesh - mesh_d) * ufl.dx)
2

You are not adding, but reassigning the geometry.
Try:
mesh_d.geometry.x[:,:mesh_d.geometry.dim] += uh.x.array.reshape((-1, mesh_d.geometry.dim))

The best way to look at the difference is to save the mesh to file and compare the two in say paraview (or with pyvista)

2 Likes
3

Yes, that worked. Thanks!

But I still have a problem with evaluating the difference between two meshes and changing the mesh using a vector.

I used the following method to find the two mesh coordinates. But after finding the difference R, I can not update the mesh with the R. Could you please help me?

 x_ref= V_ref.tabulate_dof_coordinates()
 x_def= V_def.tabulate_dof_coordinates()

R= x_def-x_ref
update_mesh=mesh_def-R
4

Why do you want to subtract the two meshes. The displacement vector uh is the difference between the unperturbed and perturbed geometry.

5

I aim to achieve a target mesh. To do this, I calculate the difference (residual R) between the current mesh and the target mesh at each iteration, and then subtract it from the current mesh for next iteration.

I read in one of your post that V.tabulate_dof_coordinates() does not give the coordinates of the nodes, but the coordinates of the degrees of freedom. I did not understand it as I get the size of out come as (number of nodes, 3) same as my mesh coordinates.

But, the size of my displacement uh is as I expected (number of degree of freedom)

Then when I do the following it is not correct and I can not have the correct subtract.

x_ref= V_ref.tabulate_dof_coordinates()
 x_def= V_def.tabulate_dof_coordinates()

R= x_def-x_ref
update_mesh.geometry.x[:,:] -= R

Could you please advice me on this?

Thanks!

6

This is because in your special case, the mesh is described with the same Finite element (linear polynomials) as your function space.

You have not supplied enough code for anyone to understand what the difference between V_ref and V_def is. Please provide a minimal reproducible example

1 Like
8

You haven’t included the perturbation file, so there is no way to reproduce your behavior