Conformal Unstructured Mesh Adaptation for Linear Elastic Tensial Loading

# Adaptive MFEM+PUMI Workflow

MFEM with PUMI Conformal Unstructured Mesh Adaptation

## At A Glance

Questions | Objectives | Key Points |

How can I define a simulation on a complex geometric model? | Demonstrate model defeaturing and mesh generation pre-processing tools | Complex models require robust pre-processing tools! |

Demonstrate the power of geometric model classification | Automated simulation tools designed for scientists and engineers working with real CAD models problem definition information must be specified on the geometric model - not the mesh. | |

How can I reliably and efficiently execute simulations? | Demonstrate an adaptive MFEM+PUMI elastic analysis | Adaptation is critical to automated, robust, and efficient simulation of simulations with transient behavior in which a static mesh defined a-priori will fail. |

Before you begin, first Open the Answers Form in a separate browser tab/window. We will be entering responses to questions here that are placed throughout the lesson.

The MFEM+PUMI examples are prepared for execution on Cooley at ALCF. Except where noted, all commands shown below should be executed in a terminal logged into Cooley.

Each step can be performed independently, but we encourage you to go
through them in order. If there are issues running the examples, or you
wish to skip ahead, you will find the execution output in the `expectedResults`

directory.

## Geometric Model Defeaturing

Geometric models are often provided by design engineers and include many features that are not required for simulation. The left side of Figure 1 depicts the geometric model of a Tokamak fusion reactor that has multiple bolts, nuts, and brackets that have no influence on the frequency analysis. These features are removed from the Parasolid CAD model by the Simmetrix SimModSuite tools by calling Parasolid kernel APIs to delete the faces and ‘heal’ the remaining hole via extensions of the bounding edges and faces. The right side of Figure 1 depicts the defeatured Tokamak model. In a model with this level of complexity that engineers and scientists running the simulation often must coordinate to create the `as-simulated’ model.

*Figure 1. Initial (left) and defeatured (right) Tokamak antenna geometric model*

For this hands-on we removed a few small faces on the RPI Formula Hybrid suspension upright shown in Figure 2 to illustrate a few key characteristics of the mesh generation procedures.

*Figure 2. RPI Formula Hybrid suspension upright*

*Figure 3. RPI Formula Hybrid suspension upright with small faces removed*

In the Mesh Generation step we will generate and compare meshes of the initial and defeatured geometric models.

### Optional - Run SimModeler to defeature the upright.

Run SimModeler within a terminal inside the VNC remote desktop session:

```
cd ~/track-5-numerical/mfem-pumi-lesson/defeaturing
./simmodeler
```

The following mouse commands, with the `Ctrl`

button held down, will manipulate
the view:

- Ctrl+Left Mouse Button: Rotate
- Ctrl+Right Mouse: Zoom
- Ctrl+Middle Mouse: Pan
- Ctrl+Wheel: Zoom

Remove a small face as shown in Figure 4.

- click ‘File->Import Geometry’ and select the
`upright.x_t`

model - click the ‘Modeling’ tab
- click ‘Delete Features’ (0)
- click a small geometric model face (1)
- click the ‘+’ button below ‘Max Num. Faces’ (2)
- click ‘Apply’ (3)

*Figure 4. Face removal example.*

Try removing the ring inside the hub by selecting all three of its faces, clicking the ‘+’ symbol to add them to the list, then clicking ‘Apply’.

## Problem Definition Introduction

The boundary conditions, initial conditions, material parameters, etc. are specified on the geometric model - not the mesh. In the Problem Definition Section we will describe the mechanism that supports this and the specific boundary conditions for this problem.

## Mesh Generation

Now that we have prepared the model and defined the boundary conditions we can proceed with mesh generation using the Simmetrix SimModSuite library APIs. The mesh generation procedures are driven by size controls defined by the user on geometric model entities or defined within geometric primitives that intersect the model (e.g., cubes, spheres, cylinders, etc.). In this example we generate a coarse mesh using a single relative size control applied to the entire model. This control guides the generation procedures to create edges whose lengths are proportional, using our prescribed value, to the length of local geometric model features.

The mesh generator is automated, which in the domain of mesh generation means that it will always produce a valid mesh while doing its best to respect the user specified controls. Note, the features of the geometric model highly dictate what mesh sizes are possible. Specifically:

- the mesh topology must exactly respect the geometric model topology, i.e., each geometric model edge has a set of meshes edges classified on it that cover the entire model edge.
- in the case of problems on curved domains, the geometry of the typical interpolating point mesh (even using higher-order curved elements) is only an approximation to the original geometric model.

Once the mesh is generated we convert it, in-memory (i.e., using PUMI and SimModSuite APIs), to a PUMI mesh and then write the PUMI mesh to file.

a

```
cd ~/track-5-numerical/mfem-pumi-lesson/meshGeneration
export SIM_LICENSE_FILE=.
# generate the mesh on the defeatured model and create paraview vtu files
./generate upright_defeatured.smd --native-model=upright_defeatured_nat.x_t case1
./render upright_defeatured.smd case1/ case1_vtu
```

*Figure 5. Mesh of initial upright model*

*Figure 6. Mesh of defeatured upright model*

### Optional - Visualize the Initial Meshes

Download the `.*vtu`

files from
Here, or
use a file transfer utility (e.g., `scp`

, `rsync`

, `putty`

, etc.)
to copy the ones you generated to your local machine.

## Partitioning

Efficient parallel execution requires equally distributing the units of work among the processing resources while keeping communication costs low. The multi-level graph and recursive sectioning geometric methods are among the most commonly available and used methods for partitioning unstructured meshes. In this lesson we will exercise the multi-level graph partitioner (referred to as ‘graph’ below) provided by the Zoltan interface and implemented by ParMETIS. We will then run recursive coordinate bisection (‘RCB’) and recursive inertial bisection (‘RIB’) to produce meshes of with the same part count (number of sub-domains), and compare the results. All three partitioners are targeting an imbalance of 5%; defined as the max element count on any part divided by the average element count across all parts.

```
cd ~/track-5-numerical/mfem-pumi-lesson/partition
export SIM_LICENSE_FILE=.
isLocal=0
render=1
runParma=0
mdl=upright_defeatured.smd
mesh=case1/
mpiexec -np 4 ./ptnParma $mdl $mesh parmetis/ 4 pmetis kway $isLocal $runParma $render &> graph.out
mpiexec -np 4 ./ptnParma $mdl $mesh rcb/ 4 rcb ptn $isLocal $runParma $render &> rcb.out
mpiexec -np 4 ./ptnParma $mdl $mesh rib/ 4 rib ptn $isLocal $runParma $render &> rib.out
```

The entity imbalance information is listed on the lines containing ```
entity
imbalance <v e f r>:
```

.
The `grep`

command to extract these lines from the outputs is:

```
grep 'afterSplit entity imbalance' *.out
```

The imbalance values follow the definition given above (`max count/avg count`

);
i.e., a 5% imbalance will be listed as `1.05`

.

The number mesh entities that are shared with other processes is indirectly
indicated on the lines containing `weighted TYPE <tot max min avg>`

where `TYPE`

is `vtx`

, `edge`

, or `face`

.
The `grep`

command to extract these lines from the outputs is:

```
grep 'afterSplit weighted vtx' *.out
```

If the total (`tot`

) or average (`avg`

) number of entities of a given TYPE is
larger relative to another partition then that partition has more replicated
entities.

Given the values of the imbalance and shared entity count for each partition
(graph, RCB, RIB) above and the content of the `[graph|rcb|rib].out`

files,

c

a

c

The partitioning decisions of RIB and RCB are based on the centroid of mesh elements while the multi-level graph partitioner is using mesh adjacency information (graph vertex = mesh element, graph edge = mesh face shared by two mesh elements).

All of the above.

*Figure 7. Partitions created with multi-level graph (left), RCB (middle), and RIB (right)*

### Optional - Visualize the Partitioned Meshes

Download the `.*vtu`

files from
Here, or
use a file transfer utility (e.g., `scp`

, `rsync`

, `putty`

, etc.)
to copy the ones you generated to your local machine.

## Problem Definition

We will define a tensile loading on the upright by applying a uniform force on one end (min Y face) and fixing the displacements on the other end (max Y face); as depicted in Figure 8. All other geometric model faces will be unconstrained. No body forces are applied.

*Figure 8. Boundary conditions*

Using geometric classification of the mesh we can define the boundary conditions on the geometric model without having any knowledge/consideration of the mesh. Geometric classification defines the equal, or greater, order association of mesh entities to geometric model entities. Specifically, a mesh region can only be classified on a geometric model region, a mesh face on a geometric face or region, and so on. Geometric classification is used throughout the set of unstructured mesh tools. For example, in mesh adaptation, classification (along with knowledge of the geometries parametric shape information) is critical for improving the meshes approximation to the geometric model when moving, adding, or removing mesh entities. Further reading on the use of classification is available in Beall et al. 1999 and Ibanez 2016.

A very simple text based interface was defined for this demonstration example to
define the two boundary conditions.
To specify the geometric model faces with
the load applied we add one line with `Load`

followed by another line with the
number of faces with that load applied. The geometric model face ids, one per
line, are then listed. Using the same convention we define the fixed faces using
the string `Dirichlet`

, the number of listed faces, then the list of face ids.
The boundary condition specification for this example is listed below.
Using the geometric model classification mechanism supports creation of more
feature-rich interfaces (command line, file, GUI).

```
$ cat upright.def
Dirichlet
2
43
838
Load
1
344
```

All of the above.

## Adaptive Simulation

This problem solves a linear elasticity model on a cantilever beam. Specifically, we approximate the weak form of \[-div(\sigma(u))=0\]

where \[\sigma(u)=\lambda*div(u)*I+\mu*(\nabla*u+u*\nabla)\]

is the stress tensor corresponding to displacement field \(u\), and \(\lambda=1\) and \(\mu=1\) are the material Lame constants.

Figure 8 depicts the applied boundary conditions; a fixed displacement \(u=0\) on the max Y face (#43), and the small circular hole in the mounting bracket marked with the red arrow (#838), and pull force \(f=1.0e-1\) on the min Z face. The specification of these boundary conditions using geometric model entity ids is described in the Problem Definition section.

To control the error without a-priori knowledge of solution we will run this example using mesh adaptation. The adaptive cycle:

- solves the problem on an initial uniform coarse mesh,
- estimates the discretization error,
- defines an isotropic size field from the error estimate,
- adapts the mesh using the size field, and
- solves the problem on the adapted mesh.

Multiple methods are available for error estimation (see the Further Reading section) and the subsequent size field computation. In this example a patch recovery approach is used.

Given the structure of the geometric model and the specified boundary conditions, we expect that the upper portion (positive Z) of the main hub experiences higher stresses as the lower part is much stiffer. Given this intuition we would expect higher error there and therefore a finer mesh is expected as the result of the meshAdapt.

Run the adaptive simulation:

```
cd ~/track-5-numerical/mfem-pumi-lesson/analysis
source setup.sh #set shared library paths
mpiexec -np 2 ./pumi_upright_ex2p -p upright_defeatured_nat.x_t -bf upright.def -m 2p_10k/
```

*Figure 9. Initial mesh (left) and final adapted mesh with stress field (right).

The requested isotropic (same in all directions) edge length is specified with a scalar value at each mesh vertex. In Section 2.1 of Li et al. edge lengths that vary with direction, anisotropic, are specified.

Using a mesh metric tensor that can be decomposed into eigen values and eigen vectors that specify the length of the edges in the principle directions.

### Optional - Visualize the analysis results

Download the `.*vtu`

files from
Here, or
use a file transfer utility (e.g., `scp`

, `rsync`

, `putty`

, etc.)
to copy the ones you generated to your local machine.

## Out-Brief

We have exercised a workflow using PUMI and MFEM software components that is based on a high level definition of the problem domain; the geometric model. The PUMI and MFEM software components exchange information on the geometric model, the mesh, and fields on the mesh to provide a robust and extensible set of functionalities that support the needs of multiple applications.

### Further Reading

- SCOREC tools
- C. W. Smith, M. Rasquin, D. Ibanez, K. E. Jansen, and M. S. Shephard
**Improving Unstructured Mesh Partitions for Multiple Criteria Using Mesh Adjacencies**, SIAM Journal on Scientific Computing, 2018. DOI: 10.1137/15M1027814 - D.A. Ibanez, E.S. Seol, C.W. Smith and M.S. Shephard,
**PUMI: Parallel Unstructured Mesh Infrastructure**, ACM Transactions on Mathematical Software, 42(3): Article 17 (28 pages), 2016. DOI: 10.1145/2814935. - M. Zhou, O. Sahni, T. Xie, M.S. Shephard and K.E. Jansen,
**Unstructured Mesh Partition Improvement for Implicit Finite Element at Extreme Scale**, Journal of Supercomputing, 59(3): 1218-1228, 2012. DOI 10.1007s11227-010-0521-0 - M. Zhou, T. Xie, S. Seol, M.S. Shephard, O. Sahni and K.E. Jansen,
**Tools to Support Mesh Adaptation on Massively Parallel Computers**, Engineering with Computers, 28(3):287-301, 2012. DOI: 10.1007s00366-011-0218-x - M. Zhou, O. Sahni, M.S. Shephard, K.D. Devine and K.E. Jansen,
**Controlling unstructured mesh partitions for massively parallel simulations**, SIAM J. Sci. Comp., 32(6):3201-3227, 2010. DOI: 10.1137090777323 - M. Zhou, O. Sahni, H.J. Kim, C.A. Figueroa, C.A. Taylor, M.S. Shephard, and K.E. Jansen,
**Cardiovascular Flow Simulation at Extreme Scale**, Computational Mechanics, 46:71-82, 2010. DOI: 10.1007s00466-009-0450-z

- C. W. Smith, M. Rasquin, D. Ibanez, K. E. Jansen, and M. S. Shephard
- Partitioning
- K. Devine, E. Boman, R. Heaphy, B. Hendrickson, and C. Vaughan,
**Zoltan data management services for parallel dynamic applications**, Computing in Science & Engineering, Volume: 4, Issue: 2, 2002, DOI: 10.1109/5992.988653. - G. Karypis and V. Kumar,
**Parallel Multilevel k-way Partitioning Scheme for Irregular Graphs**, SIAM Review, Vol. 41, No. 2, pp. 278 - 300, 1999, DOI: 10.1137/S0036144598334138.

- K. Devine, E. Boman, R. Heaphy, B. Hendrickson, and C. Vaughan,
- Mesh data and geometry interactions
- Seol, E.S. and Shephard, M.S.,
**Efficient distributed mesh data structure for parallel automated adaptive analysis**, Engineering with Computers, 22(3-4):197-213, 2006, DOI: 10.1007s00366-006-0048-4 - Beall, M.W., Walsh, J. and Shephard, M.S,
**A comparison of techniques for geometry access related to mesh generation**, Engineering with Computers, 20(3):210-221, 2004, DOI: 10.1007s00366-004-0289-z. - Beall, M.W.. and Shephard, M.S.,
**An Object-Oriented Framework for Reliable Numerical Simulations**, Engineering with Computers, 15(1):61-72, 1999, DOI: 10.1007/s003660050005. - Beall, M.W. and Shephard, M.S.,
**A general topology-based mesh data structure**, Int. J. Numer. Meth. Engng., 40(9):1573-1596, 1997, DOI: 10.1002(SICI)1097-0207(19970515)40:9<1573::AID-NME128>3.0.CO;2-9.

- Seol, E.S. and Shephard, M.S.,
- Adaptivity
- “Morteza H. Siboni, and Mark S. Shephard,
**Adaptive Workflow for Simulation of RF Heaters**, 2021. Submitted to Computer Physics Communications.”:https://www.scorec.rpi.edu/REPORTS/2021-3.pdf - “Aleksandr Ovcharenko,
**Parallel Anisotropic Mesh Adaptation with Boundary Layers**, Ph.D. Dissertation, RPI, 2012”:http://www.scorec.rpi.edu/REPORTS/2012-20.pdf - Q. Lu, M.S. Shephard, S. Tendulkar and M.W. Beall,
**Parallel Curved Mesh Adaptation for Large Scale High-Order Finite Element Simulations**, Proc. 21 Roundtable, Springer, NY, pp. 419-436, 2012, DOI 10.1007978-3-642-33573-0. - A. Ovcharenko, K. Chitale, O. Sahni, K.E. Jansen and M.S. Shephard, S. Tendulkar and M.W. Beall,
**Parallel Adaptive Boundary Layer Meshing for CFD Analysis**, Proc. 21st International Meshing Roundtable, Springer, NY, pp. 437-455, 2012, DOI 10.1007978-3-642-33573-0 - X.-J. Luo, M.S. Shephard, L.-Q. Lee and C. Ng,
**Moving Curved Mesh Adaption for Higher Order Finite Element Simulations**, Engineering with Computers, 27(1):41-50, 2011. DOI: 10.1007/s00366-010-0179-5 - O. Sahni, X.J. Luo, K.E. Jansen, M.S. Shephard,
**Curved Boundary Layer Meshing for Adaptive Viscous Flow Simulations**, Finite Elements in Analysis and Design, 46:132-139, 2010. DOI: 10.1007/s00366-008-0095-0 - Alauzet, F., Li, X., Seol, E.S. and Shephard, M.S.,
**Parallel Anisotropic 3D Mesh Adaptation by Mesh Modification**, Engineering with Computers, 21(3):247-258, 2006, DOI: 10.1007s00366-005-0009-3 - Li, X., Shephard, M.S. and Beall, M.W.,
**3-D Anisotropic Mesh Adaptation by Mesh Modifications**, Comp. Meth. Appl. Mech. Engng., 194(48-49):4915-4950, 2005, doi:10.1016/j.cma.2004.11.019 - Li, X., Shephard, M.S. and Beall, M.W.,
**Accounting for curved domains in mesh adaptation**, International Journal for Numerical Methods in Engineering, 58:246-276, 2003, DOI: 10.1002/nme.772

- “Morteza H. Siboni, and Mark S. Shephard,