Add support for dual triangle/tetra-regions

[adtzlr]

No description provided.

[adtzlr]

Example 1:

import felupe as fem
import numpy as np

line = fem.mesh.Line(a=15, b=30, n=6)
mesh = fem.mesh.revolve(line, phi=180, axis=2, n=13)
bush = fem.mesh.triangulate(mesh)
mesh = fem.mesh.add_midpoints_edges(bush)

# fix midpoints
pts = fem.mesh.revolve(line, phi=180, axis=2, n=25).points
for p in pts:
    x = np.linalg.norm(mesh.points, axis=1)
    dx = np.linalg.norm(mesh.points - p, axis=1)
    i = np.argmin(dx)
    if dx[i] / x[i] < 0.01:
        mesh.points[i] = p
        
region = fem.RegionQuadraticTriangle(mesh)
field = fem.FieldsMixed(region, n=3, planestrain=True)
umat = fem.ThreeFieldVariation(fem.NeoHooke(mu=1, bulk=5000))

solid = fem.SolidBody(umat, field)
inner = np.isclose(np.linalg.norm(mesh.points, axis=1), 15, rtol=0.01)
outer = np.isclose(np.linalg.norm(mesh.points, axis=1), 30, rtol=0.01)
bounds = dict(
    sym=fem.Boundary(field[0], fy=0, skip=(1, 0)),
    inner=fem.Boundary(field[0], mask=inner, skip=(1, 0)),
    outer=fem.Boundary(field[0], mask=outer),
    move=fem.Boundary(field[0], mask=inner, skip=(0, 1)),
)
step = fem.Step(
    items=[solid], 
    ramp={bounds["move"]: fem.math.linsteps([0, 5], 50)}, 
    boundaries=bounds
)
job = fem.CharacteristicCurve(steps=[step], boundary=bounds["move"])
job.evaluate(filename="bushing.xdmf")
job.plot(yscale=2)

(unstable)

[codecov]

Codecov Report

Base: 97.36% // Head: 97.36% // Increases project coverage by +0.00% 🎉

Coverage data is based on head (02b5f23) compared to base (2bfe392).
Patch coverage: 100.00% of modified lines in pull request are covered.

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #373   +/-   ##
=======================================
  Coverage   97.36%   97.36%           
=======================================
  Files          81       81           
  Lines        3568     3574    +6     
=======================================
+ Hits         3474     3480    +6     
  Misses         94       94           

Impacted Files	Coverage Δ
felupe/region/_templates.py	99.03% <100.00%> (+0.05%)	⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[adtzlr]

and the 2D plane-strain MINI-element:

import felupe as fem
import matadi as mat
import numpy as np

line = fem.mesh.Line(a=15, b=30, n=21)
mesh = fem.mesh.revolve(line, phi=180, axis=2, n=21)
bush = fem.mesh.triangulate(mesh)
mesh = fem.mesh.add_midpoints_faces(bush)
mesh.points[mesh.npoints - mesh.ncells:] = 0
mesh.cell_type = "triangle"

# fix midpoints
pts = fem.mesh.revolve(line, phi=180, axis=2, n=41).points
for p in pts:
    x = np.linalg.norm(mesh.points, axis=1)
    dx = np.linalg.norm(mesh.points - p, axis=1)
    i = np.argmin(dx)
    if dx[i] / x[i] < 0.01:
        mesh.points[i] = p
        
region = fem.RegionTriangleMINI(mesh)
field = fem.FieldsMixed(region, n=3, planestrain=True)
umat = fem.ThreeFieldVariation(fem.NeoHooke(mu=1, bulk=5000))

solid = fem.SolidBody(umat, field)
inner = np.isclose(np.linalg.norm(mesh.points, axis=1), 15, rtol=0.01)
outer = np.isclose(np.linalg.norm(mesh.points, axis=1), 30, rtol=0.01)
bounds = dict(
    sym=fem.Boundary(field[0], fy=0, skip=(1, 0)),
    inner=fem.Boundary(field[0], mask=inner, skip=(1, 0)),
    outer=fem.Boundary(field[0], mask=outer),
    move=fem.Boundary(field[0], mask=inner, skip=(0, 1)),
)
step = fem.Step(
    items=[solid], 
    ramp={bounds["move"]: fem.math.linsteps([0, 9], 10)}, 
    boundaries=bounds
)
job = fem.CharacteristicCurve(steps=[step], boundary=bounds["move"])
job.evaluate(filename="bushing.xdmf", mesh=bush.as_meshio())
job.plot(yscale=2)

(quite bad large strain prediction)