Transformations, extrusions and other operations#

Geometrical transformations#

XLiFE++ allows you to apply geometrical transformations on Mesh, Geometry and Geometry children objects. The main type is Transformation. It can be a canonical transformation or a composition of transformations.

Canonical transformations#

In the following, we will consider straight lines and planes. A straight line is fully defined by a point and a direction. The latter is a vector of components (2 or 3). This is a reason why we will write a straight line as follows : \(\left(\Omega, \vec{d}\right)\). A plane is fully defined by a point and a normal vector. This is a reason why we will write a plane as follows : \(\left[\Omega, \vec{n}\right]\).

All canonical transformations are related to classes inheriting from the Transformation class and the way to define them is based on a key-value system.

Translations#

Point P is the image of point M by a translation of vector \(\vec{u}\) if and only if

\[\overrightarrow{MP}=\vec{u}\]

A translation can be defined by a vector (size 1, 2 or 3) or by its components, thanks to the key _direction:

Reals u;
Real ux, uy, uz;
Translation t1(_direction=u), t2(_direction={ux, uy, uz});

Important

Default value for _direction is the 3D zero vector. In the second syntax, uy and uz can be omitted too. If so, their default value is 0.

2D rotations#

Point P is the image of point M by the 2D rotation of center \(\Omega\) and of angle \(\theta\) if and only if

\[\begin{split}\overrightarrow{\Omega P} = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} \overrightarrow{\Omega M}\end{split}\]

A 2D rotation is defined by a point and an angle (in radians), thanks to keys _center and _angle:

Point omega;
Real angle;
Rotation2d r(_center=omega, _angle=angle);

Important

Default value for _center is the 2D zero vector. Default value for _angle is 0.

3D rotations#

Point P is the image of point M by the 3D rotation of axis \(\left(\Omega, \vec{d}\right)\) and of angle \(\theta\) (in radians) if and only if

\[\overrightarrow{\Omega P} = \cos\theta \; \overrightarrow{\Omega M} + \left( 1 -\cos\theta\right) \; \overrightarrow{\Omega M} \cdot \vec{n} + \sin\theta \; \vec{n} \wedge \overrightarrow{\Omega M} \quad \text{(Rodrigues rotation formulae)}\]

where \(\displaystyle \vec{n}=\dfrac{\vec{d}}{||\vec{d}||}\) (the unitary direction).

A 3D rotation is defined by a point, a direction vector and an angle (in radians), thanks to keys _center, _axis and _angle:

Point omega;
Reals d;
Real dx, dy, dz;
Real theta;
Rotation3D r1(_center=omega, _axis=d, _angle=theta), r2(_center=omega, _axis={dx, dy, dz}, _angle=theta);

Important

Default value for _center and _axis is the 3D zero vector. Default value for _angle is 0. In the second syntax, dz can be omitted. If so, its default value is 0.

Homotheties#

Point P is the image of point M by the homothety of center \(\Omega\) and of factor \(k\) if and only if

\[\overrightarrow{\Omega P} = k \; \overrightarrow{\Omega M}\]

A homothety is defined by a point and a scale factor, thanks to keys _center and _scale:

Point omega(1.,2.,3.);
Real k=2.;
Homothety h(_center=omega, _scale=k);

Important

Default value for _scale is 1. Default value for _center is the 3D zero vector.

Point reflections#

Point P is the image of point M by the point reflection of center \(\Omega\) if and only if

\[\overrightarrow{\Omega P} = - \overrightarrow{\Omega M}\]

It is a homothety of factor -1 and same center.

Point omega(1.,2.,3.);
PointReflection h(_center=omega); // omega can still be omitted, as for homothety

2D reflections#

Point P is the image of point M by the 2D reflection of axis \(\left(\Omega,\vec{d}\right)\) if and only if

\[\overrightarrow{MP} = 2 \overrightarrow{MH}\]

where \(H\) is the orthogonal projection of \(M\) on \(\left(\Omega,\vec{d}\right)\)

A 2D reflection is then defined by a Point and a vector defining an axis, thanks to keys _center and _direction:

Point omega(1.,2.,3.);
Vector<Real> d(1.,0.,0.);
Real dx=1., dy=0.;
Reflection2d r1(_center=omega, _direction=d), r2(_center=omega, _direction={dx, dy});

Important

Default value for _direction is the 2D zero vector. Default value for _center is the 2D zero point.

3D reflections#

Point P is the image of point M by the 3D reflection of plane \(\left[\Omega,\vec{n}\right]\) if and only if

\[\overrightarrow{MP} = 2 \overrightarrow{MH}\]

where \(H\) is the orthogonal projection of \(M\) on plane \(\left[\Omega,\vec{n}\right]\).

A 3D reflection is then defined by a Point and a normal vector defining a plane, thanks to keys _center and _normal:

Point omega(1.,2.,3.);
Vector<Real> n;
Real nx, ny, nz;
Reflection3d r1(_center=omega, _normal=n), r2(_center=omega, _normal={nx, ny, nz});

Important

Default value for _normal is the 3D zero vector. Default value for _center is the 3D zero point.

Linear transformations#

Any linear transformation can be set from its matrix-vector (\(x\rightarrow \mathbb{A}x+b\)) using the following Transformation constructors:

Transformation(const Matrix<Real>& A, const Vector<Real>& b = Vector<Real>(3,0.))
Transformation(Real a11, Real a12, Real a13, Real a21, Real a22, Real a23,
               Real a31, Real a32, Real a33, Real b1=0, Real b2=0, Real b3=0);

Note

The matrix representation of a transformation can be accessed using the member functions mat() and vec(). Matrix representations are also available for canonical geometries.

Composition of transformations#

To define a composition of transformations, use the operator * between canonical transformations, an is the following example:

Rotation2d r1(_center=Point(0.,0.), _angle=120.);
Reflection2d r2(_center=Point(1.,-1.), _axis={1.,2.5, -3.});
Translation t1(_direction={-1.,4.});
Homothety h (_center=Point(-1.,0.), _scale=-3.2);
Transformation t = r1*h*r2*t1;

Composition * has to be understood as usual composition operator \(\circ\) : \(t(P)=r1(h(r2(t1(P))))\).

Note

The matrix representation (mat() and vec() member functions) is available for composition of transformations.

Transformations on Point, geometries and Mesh#

How to apply a transformation ?#

On Point and geometries#

Transformations can be applied on every canonical geometry, composite geometry or loop geometry. The only exceptions are Rotation2d and Reflection2d for 3D geometries, as they are meaning less in theses cases.

In the following, to illustrate application of transformations, we will use Segment.

The general case is to use transform() :

  • either as a member function, to modify the object on which the transformation is applied

    Segment s;
s.transform(Translation(_direction={0.,0.,1.}));

    Warning

    This is only available for geometries, not for Point

  • or as an external function, to build a new object next to the one on which the transformation is applied

    Segment s;
Segment s2 = transform(s, Translation(_direction={0.,0.,1.}));

Shortcuts are naturally available in both contexts (only the external functions for Point) with routines translate(), rotate2d(), rotate3d(), homothetize(), pointReflect(), reflect2d() and reflect3d().

These routines also work with a key-value system, the same as the corresponding transformation object:

Segment s;
Segment s2 = translate(s, _direction={0.,0.,1.});
s.translate(_direction={0.,0.,1.});

On meshes#

On Mesh objects, transformations apply in the same way as geometries and points (See above).

Mesh m;
Mesh m2 = translate(m, _direction={0.,0.,1.});
m.translate(_direction={0.,0.,1.});

There is a little difference concerning the external shortcuts, in so are as additional keys can be used: _name to set the Mesh name (as in every Mesh constructor), and _suffix to add a suffix to every domain name. In the following code sample, plain domain in Mesh m2 will be S_prime and side domains will be Sigma_prime and Gamma_prime:

Mesh m(Segment(_v1=Point(0.,0.), _v2=Point(1.,2.), _domain_name="S"), _nnodes=10, _side_names={"Sigma, "Gamma"});
Mesh m2 = translate(m, _direction={0.,1.}, _name="Translation of m", _suffix="prime");

What does a transformation really do ?#

Applying a transformation to an object means computing the image of each point defining the object. It also calls the transformation routine on some underlying objects, such as the Geometry pointer in a Mesh, or the BoundingBox and the MinimalBox in a Geometry.

Extrusions of geometries and meshes#

This is another way to define geometries: by extrusion of geometries of lesser dimension. Extruded geometries can be surfaces or volumes, defined by a geometry (the section of the extruded geometry) and a geometrical transformation. This feature can be used to generate meshes with the GMSH interface with some restrictions about the transformation : only translations or rotations are authorized. There is also another parameter : the number of layers. Let’s see the following figures:

On the left, extrusion of a disk by a translation, with 3 layers. On the right, extrusion of a circular arc by rotation, with 4 layers

Figure made with TikZ

Figure made with TikZ

Extrusion of geometries#

It is possible to extrude a 1D or a 2D geometry, applying a Transformation (Translation, Rotation, …) using the routine extrude(). More precisely, you can extrude:

  • a canonical geometry (1D or 2D). Here, a CircArc:

    Point b(1.5,-4.,0.);
Point c(2.,-3.5,0.);
CircArc g(_center=Point(1.5,-3.5,0.), _v1=b, _v2=c, _nnodes=5, _domain_name="BC");
Geometry e2D=extrude(g, Translation(_direction={0.,0.,4.}), _layers=5, _domain_name="Omega");
    extrusion1DByTranslation (png)

    Fig. 87 Extrusion of a circular arc by translation, with 5 layers.#

  • a loop geometry (1D or 2D). Here, a rounded rectangle defines as in Definition of a geometry from its boundary:

    Point a(-1.5,-4.,0.); Point b(1.5,-4.,0.); Point c(2.,-3.5,0.); Point d(2.,3.5,0.);
Point e(1.5,4.,0.); Point f(-1.5,4.,0.); Point g(-2.,3.5,0.); Point h(-2.,-3.5,0.);
Segment s1(_v1=a, _v2=b, _nnodes=21, _domain_name="AB");
CircArc c1(_center=Point(3.5,0.5,0.), _v1=b, _v2=c, _nnodes=5, _domain_name="BC");
Segment s2(_v1=c, _v2=d, _nnodes=11,  _domain_name="CD");
CircArc c2(_center=Point(3.5,1.5,0.), _v1=d, _v2=e, _nnodes=5, _domain_name="DE");
Segment s3(_v1=e, _v2=f, _nnodes=21,  _domain_name="EF");
CircArc c3(_center=Point(0.5,1.5,0.), _v1=f, _v2=g, _nnodes=5, _domain_name="FG");
Segment s4(_v1=g, _v2=h, _nnodes=11,  _domain_name="GH");
CircArc c4(_center=Point(0.5,0.5,0.), _v1=h, _v2=a, _nnodes=5, _domain_name="HA");
Geometry g=planeSurfaceFrom(s1+c1+s2+c2+s3+c3+s4+c4,"Omega");
Geometry e3d=extrude(g, Translation(_direction={0.,0.,4.}), _layers=10, _domain_name="Omega");
    extrusion2DloopByTranslation (png)

    Fig. 88 Extrusion of a rounded rectangle (loop geometry) by a rotation, with 10 layers.#

  • Every composite geometry composed exclusively of a geometry and its holes (1D or 2D). That is to say only operator- or operator-= is used to define the geometry:

    Ellipse e1(_center=Point(0.,0.,0.), _v1=Point(4,0.,0.), _v2=Point(0.,5.,0.), _nnodes=12, _domain_name="Omega1", _side_names=Strings("Gamma_1","Gamma_2","Gamma_3","Gamma_4"));
Ellipse e2(_center=Point(1.,2.,0.), _v1=Point(1.5,2.,0.), _v2=Point(1.,3.,0.), _nnodes=12, _domain_name="Omega3", _side_names=Strings("Gamma_9","Gamma_10","Gamma_11","Gamma_12"));
Geometry e3d2=extrude(e1-e2, Rotation3d(_center=Point(5.,0.,0.), _axis={0.,5.,0.}, _angle=pi_/2.), _layers10, _domain_name="Omega", _side_names="Gamma");
    extrusion2DCompositeByRotation (png)

    Fig. 89 Extrusion of an ellipse with an elliptical hole by rotation, with 10 layers.#

How to define names of lateral domains of an extrusion ?#

Instead of giving the same name to every lateral surface of an extrusion, it is possible to name each of them, but what about sides numbering ?

First example, let’s take the extrusion of a CircArc :

Point b(1.5,-4.,0.);
Point c(2.,-3.5,0.);
CircArc g(_center=Point(1.5,-3.5,0.), _v1=b, _v2=c, _nnodes=5, _domain_name="BC");
Geometry e2D=extrude(g, Translation(_direction={0.,0.,4.}), _layers=5, _domain_name="Omega", _side_names=Strings("Gamma1", "Gamma2"));

Point b is the front below left corner and point c is the front top right corner of the circular arc g. As g is defined from b to c, the first lateral side, corresponding to domain Gamma1, will be the edge below, starting from b. If g was defined from c to b, Gamma1 would have corresponded to the edge above, starting from c.

Second example, let’s take the extrusion of an ellipse with an elliptical hole:

Ellipse e1(_center=Point(0.,0.,0.), _v1=Point(4,0.,0.), _v2=Point(0.,5.,0.), _nnodes=12, _domain_name="Omega1");
Ellipse e2(_center=Point(1.,2.,0.), _v1=Point(1.5,2.,0.), _v2=Point(1.,3.,0.), _nnodes=12, _domain_name="Omega3");
Geometry e3d3=extrude(e1-e2, Rotation3d(_center=Point(5.,0.,0.), _axis={0.,5.,0.}, _angle=pi_/2.), _layers=10, _domain_nae="Omega", _side_names={"Gamma1", "Gamma2", "Gamma3", "Gamma4", "Gamma5", "Gamma6", "Gamma7", "Gamma8});

This time, lateral surfaces are ordered as follows:

  • Lateral surfaces from the outer ellipse are ordered the same way as borders of the ellipse. They will become Gamma1, Gamma2, Gamma3 and Gamma4.

  • Lateral surfaces from the inner ellipse (and every hole in general) are ordered in the reverse order of borders of the ellipse. They will become Gamma5, Gamma6, Gamma7 and Gamma8.

Tip

Contrary to Gmsh, extrusion of a geometry by rotation of angle greater than \(\pi\) is available, by splitting extrusion in 2 half extrusions when angle is not \(2\pi\) or in 4 quarter extrusions when angle is \(2\pi\). As a result, the number of lateral surfaces is multiplied by 2 or 4.

Extrusion of a mesh#

As an alternative to mesh an extruded geometry, it is possible to extrude a 1D or a 2D mesh, applying a Transformation (Translation, Rotation, …) using the routine extrude():

Mesh mdi(Disk(_center=Point(0.,0.),_radius=1.,_nnodes=8),_shape=_triangle);
Mesh m=extrude(mdi, Rotation3d(_center=Point(3.,0.,0.), _axis={0.,1.,0.}, _angle=2*pi_/40), _layers=40);

Here a 3D rotation of angle \(\frac{\pi}{20}\) is applied to each section and _layers key-value specifies the number of layers (40). Note that in that case, a full turn is done and the last section coincides with the first section, the nodes being not duplicated.

mesh_extrusion_torus

Fig. 90 Extrusion using a rotation leading to a torus#

The transformation may be a basic transformation (translation, rotation, …) or a combination of transformations (see Geometrical transformations for more details):

Mesh mr(Rectangle(_xmin=0, _xmax=1, _ymin=1, _ymax=3, _nnodes=Numbers(6, 10)), _shape=_triangle, _order=1, _generator=_structured);
Mesh m=extrude(mr, Translation(0.,0.1,0.)*Rotation3d(Point(3.,0.,0.), Reals(0.,1.,0.), 2.5*pi_/100), _layers=100);
mesh_extrusion_helix

Fig. 91 Extrusion using a combination of transformations#

The extruded mesh is a mesh of:

  • quadrangles if the section mesh is 1D,

  • prisms if the section mesh is made of triangles,

  • hexahedra if the section mesh is made of quadrangles.

Instead of specifying a transformation, it is possible to specify the list of transformations to apply to the sections.

Mesh mr(Rectangle(_xmin=0, _xmax=1, _ymin=1, _ymax=3, _nnodes=Numbers(6, 10)), _shape=_triangle, _order=1, _generator=_structured);
Translation tran(_direction={0.2,0.,0.});
Rotation3d rot1(_center=Point(3.,0.,0.), _axis={0.,1.,0.}, _angle=pi_/40);
Rotation3d rot2(_center=Point(5.,0.,6.), _axis={0.,1.,0.}, _angle=-pi_/40);
std::vector<Transformation*> trs(50);
for (number_t i=0;i<20;i++) { trs[i]=&rot1; trs[i+30]=&rot2; }
for (number_t i=0;i<10;i++) { trs[i+20]=&tran; }
Mesh m=extrude(mesh2DP1, trs, _naming_domain=1);
mesh_extrusion_S

Fig. 92 Extrusion using a list of Transformations#

To deal with more general extrusions, a C++ function giving the parametrization of the extrusion path: \(t\in[0,1]\rightarrow \gamma(t)\in \mathbb R^d\) can be given. Be caution, the parametrization has to statisfy \(\gamma(0)=0\) and \(\gamma'(0)=n_s\) where \(n_s\) is the normal vector to the section to be extruded. The extrusion works as following:

  • the nodes \(M_0\) of section \(S_0\) are translated: \(M_k=M_0+\gamma(\frac{k}{n})\)

  • then, a rotation defined from the center \(C_k=C_0+\gamma(\frac{k}{n})\) and the two vectors \(n_0,\ \gamma'(\frac{k}{n})\) is applied to the nodes \(M_k\)

By default, the rotation center \(C_0\) is the barycenter of all nodes of the section \(S_0\), but it can be changed using the key-value _center. Note that this process guarantees that all sections are orthogonal to the extrusion curve.

Reals gamma(const Point& P, Parameters& par, DiffOpType dop)
{
  real_t t=P(1);
  if (dop==_dt) return Reals{16*pi_*std::sin(4*pi_*t), 0., 40.};
  return Reals{2*(1-std::cos(4*pi_*t)), 0., 40*t};
}

Mesh mdi(Disk(_center=Point(0.,0.), _radius=1., _nnodes=8), _shape=_triangle, _generator=subdiv);
Mesh m=extrude(mdi, gamma, _layers=100);
mesh_extrusion_snake

Fig. 93 Extrusion using a curve parametrization#

Tip

The orientation of normal vector to the section to be extruded is not predictable. But using the routine setNormalOrientation() it can be set:

Mesh mdi(Disk(_center=Point(0.,0.),_radius=1.,_nnodes=8),_shape=_triangle);
mdi.domain(0).setNormalOrientation(_towardsInfinite,Point(0.,0.,-1.));

where the given point (here (0.,0.,-1.)) specifies an ‘interior’ point.

Danger

If curvature radius is too small compared to the section size, the extrusion process may lead to incoherent domain!

extrude() manages the following key-values:

  • _layers to give the number of layers (mandatory)

  • _center to give the center of rotation when using a curve parametrization (computed by default)

  • _init_transformation to specify a transformation to apply to the reference section before extruding it (none by default)

  • _naming_domain to control the naming of domains (0,1 or 2)

  • _naming_section to control the naming of sections (0,1 or 2)

  • _naming_side to control the naming of boundaries (0,1 or 2)

The meaning of the naming level is

  • 0: domains/sections/boundaries are not created

  • 1: one extruded domain/section/boundary is created for any domain/section/boundary of the original section

  • 2: for each layer, one extruded domain/section/boundary is created for any domain/section/boundary of the original section

Be cautious, the side domains of extruded domain are created from the side domains of the given section. Thus, if the given section has no side domains, the extruded domains will have no side domain! The naming convention is the following:

  • Domains and side domains keep their name with the extension “_e” or “_e1”, “_e2”, …, “_en”

  • Section domains have the name of original domains with the extension “_0”, …, “_n”

  • When _naming_domain=0, the full domain is always created and named “Omega”.

The following figure illustrates the naming rules of domains.

mesh_extrusion

For instance, to mesh a tubular domain using the mesh extrusion of a crown:

Disk dext(_center=Point(0.,0.), _radius=1., _nnodes=20, _domain_name="Omega", _side_names="Sigma");
Disk dint(_center=Point(0.,0.), _radius=0.5, _nnodes=10, _side_names="Gamma");
Mesh crown(dext-dint, _generator=gmsh);
Mesh tube(crown, Translation(_direction={0,0,.1}), _layers=10, _naming_domain=1, _naming_section=1, _naming_side=1);
crown_extrude

Fig. 94 Tube prismatic mesh from extrusion of a crown mesh#

Example: definition of a conesphere#

To define geometries based on cones, extrusions have to be used. It is the case for the conesphere:

Real rb=1., hc=3.;
Real hs=rb*rb/hc;etrahedra
Real rs=sqrt(rb*rb + hs*hs);

Point origin(0.,0.,0.), apex(0.,0.,hc), p1(rb,0.,0.), p2(0.,0.,-hs-rs);

Segment s1(_v1=p1, _v2=apex, _hsteps=0.05);
Segment s2(_v1=apex, _v2=origin, _hsteps=0.05);
Segment s3(_v1=origin, _v2=p2, _hsteps=0.05);
CircArc c1(_center=Point(0.,0.,-cssphereheight), _v1=p2, _v2=p1, _hsteps=0.05);
Disk d1(_center=0.5*p1, _v1=0.5*p1+Point(0.2*rb,0.,0.), _v2=0.5*p1+Point(0.,0.,0.2*rb), _domain_name="Sigma", _hsteps=0.05);

Geometry base=planeSurfaceFrom(s2+s3+c1+s1, "Gamma");
Geometry g=extrude(base, Rotation3d(_center=Point(0.,0.,0.), _axis={0.,0.,1.}, _angle=2.*pi_), _domain_name="Omega1", _side_names={"Gamma1", "Gamma2", "Gamma3", "Gamma4", "Gamma5", "Gamma6", "Gamma7", "Gamma8"});
conesphere (png)

Fig. 95 Mesh of a conesphere.#

Other Mesh operations#

Mesh refinement#

It is possible to refine an existing mesh of order 1: a new mesh is created using the subdivision generator mentioned above. It is a uniform refinement, say, all the elements are subdivided in a same way. The corresponding constructor is defined as follows:

Disk d(_center=Point(0.,0.),_radius=1.,_nnodes=5);
Mesh md(d,_generator=gmsh,_shape=triangle);
Mesh md1(md,_refinement_depth=1);
Mesh md2(md,_refinement_depth=2);

This builds the meshes md2 and md2 which are obtained by subdividing, respectively, once and twice the mesh md:

meshRefinement

Fig. 96 Refine once and twice the P1 mesh of a disk.#

The optional keys are:

_refinement_depth: the refinement depth of the subdivision generator. 0 means no refinement (the default value).

_order: by using it, it possible to increase the order of the refined mesh (defaut value is 1).

Note

Up to now, there is no adaptative refinement tool!

Splitting mesh element#

Sometimes it may be useful to split elements into elements of another type, for instance to produce mesh of pyramids that are not provided by standard meshing software products. To do this, the routine split() is available:

Cuboid c(_xmin=0,_xmax=1,_ymin=1,_ymax=3,_zmin=1,_zmax=5,_nnodes={3, 5, 9});
Mesh meshQ1(c, _shape=hexahedron, _generator=structured, _name="Q1 mesh");
Mesh meshPyramid=split(meshQ1, pyramid, "Py1 mesh");
mesh3dQ1
mesh3DPyramidSplit

Fig. 97 Hexahedron mesh split to pyramid mesh.#

Note

Up to now, only one splitting process is available: hexahedron of order 1 into six pyramids of order 1, based on the hexahedron faces and the centroid of the hexahedron as tip.

Merging meshes#

It is possible to merge some meshes of same dimension. They can share common parts or common boundaries but with same elements/nodes. If elements or nodes on common parts do not coincide, the resulting mesh would not be conform to do finite element approximations. The merging process combines the elements and nodes of the meshes, eliminating duplicate elements or duplicate nodes. Domains with same name are also merged, names being preserved. The new whole domain is created with default name #Omega or one given and shared boundary domains are merged on demand (default).

Mesh m1(...);
Mesh m2(...);
Mesh m3(...);
m1.merge(m2,[mergeSharedBoundary=true],[name="#Omega"]);//merging m2 with m1
m1.merge(m3,[mergeSharedBoundary=true],[name="#Omega"]);//merging m2, m3 with m1
Mesh m = merge(m1,m2,m3,[mergeSharedBoundary = true],[name = "#Omega"]);
// alternate syntax merging m1,m2 and m3 in a new mesh m (up to 5 meshes)

As an exemple:

Rectangle R1(_xmin=0, _xmax=1, _ymin=1, _ymax=3, _nnodes={6,10});
Rectangle R2(_xmin=1, _xmax=2, _ymin=1, _ymax=3, _nnodes={3, 5});
Mesh mP1(R1, _shape=triangle, _generator=structured, _name="P1 mesh");
Mesh mQ1(R2, _shape=quadrangle, _generator=structured, _name="Q1 mesh");
Mesh mP1Q1 = merge(mP1,mQ1);

Rectangle R3(_xmin=1, _xmax=2, _ymin=1, _ymax=3, _nnodes={6,10});
Mesh mQ1a(R3, _shape=quadrangle, _generator=structured, _name="Q1 mesh a");
mP1Q1 = merge(mP1,mQ1a);
mergeP1Q1

Fig. 98 Non-conforming and conforming merging of a triangular mesh with a quadrangular mesh.#

In merging process, the following rules are applied:

  • all named domains of original meshes are preserved, with the same names,

  • if two different domains have the same name, say domname, they are merged together

  • if two side domains define an interface and do not have the same name, say domname1 and domname2:

    • if mergeSharedBoundary is true, a new domain named domname1_domname2 is created and the original ones are deleted

    • if mergeSharedBoundary is false, the original domains domname1 and domname2 are retained

  • if two side domains define an interface and have the same name, say domname, only one domain remains, again named domname and regarded as an interface.

Attention

To be merged, two side domains must have exactly the same nodes.

The following example shows the merging of two rectangles with different side names (mergeSharedBoundary set to true by default):

Rectangle R1(_xmin=0, _xmax=1, _ymin=0, _ymax=1, _nnodes={3,3}, _side_names={"g1","g2","g3","g4"});
Rectangle R2(_xmin=1, _xmax=2, _ymin=0, _ymax=1, _nnodes={3,3}, _side_names={"s1","s2","s3","s4"});
Mesh mR1P1(R1, _shape=_triangle, _order=1, _generator=_structured);
Mesh mR2Q1(R2, _shape=_quadrangle, _order=1, _generator=_structured);
Mesh mR1R2=merge(mR1P1,mR2Q1);

producing the domains:

Domain 'g1' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 13: side 3 of element 1
 geometric side element 14: side 3 of element 3
Domain 'g3' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 17: side 3 of element 6
 geometric side element 18: side 3 of element 8
Domain 'g4' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 19: side 2 of element 1
 geometric side element 20: side 2 of element 5
Domain 's1' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 21: side 1 of element 9
 geometric side element 22: side 1 of element 10
Domain 's2' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 23: side 2 of element 10
 geometric side element 24: side 2 of element 12
Domain 's3' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 25: side 3 of element 11
 geometric side element 26: side 3 of element 12
Domain 'g2_s4' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 29: side 2 of element 4, side 4 of element 9
 geometric side element 30: side 2 of element 8, side 4 of element 11
Domain '#Omega' of dimension 2 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 12 elements
 geometric element 1 (material 0, color 0, twin 0): triangle_Lagrange_1
 geometric element 2 (material 0, color 0, twin 0): triangle_Lagrange_1
 ...
 geometric element 12 (material 0, color 0, twin 0): quadrangle_Lagrange_1

This is the same example but with mergeSharedBoundary set to false:

Mesh mR1R2=merge(mR1P1,mR2Q1, false);

producing the domains:

Domain 'g1' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 13: side 3 of element 1
 geometric side element 14: side 3 of element 3
Domain 'g2' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 15: side 2 of element 4
 geometric side element 16: side 2 of element 8
Domain 'g3' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 17: side 3 of element 6
 geometric side element 18: side 3 of element 8
Domain 'g4' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 19: side 2 of element 1
 geometric side element 20: side 2 of element 5
Domain 's1' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 21: side 1 of element 9
 geometric side element 22: side 1 of element 10
Domain 's2' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 23: side 2 of element 10
 geometric side element 24: side 2 of element 12
Domain 's3' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 25: side 3 of element 11
 geometric side element 26: side 3 of element 12
Domain 's4' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 27: side 4 of element 9
 geometric side element 28: side 4 of element 11
Domain '#Omega' of dimension 2 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 12 elements
 geometric element 1 (material 0, color 0, twin 0): triangle_Lagrange_1
 geometric element 2 (material 0, color 0, twin 0): triangle_Lagrange_1
 ...
 geometric element 12 (material 0, color 0, twin 0): quadrangle_Lagrange_1

Finally, the case when the side domains defining the interface have the same name (here g2):

Rectangle R1(_xmin=0, _xmax=1, _ymin=0, _ymax=1, _nnodes={3,3}, _side_names={"g1","g2","g3","g4"});
Rectangle R2(_xmin=1, _xmax=2, _ymin=0, _ymax=1, _nnodes={3,3}, _side_names={"s1","s2","s3","g2"});
Mesh mR1P1(R1, _shape=_triangle, _order=1, _generator=_structured);
Mesh mR2Q1(R2, _shape=_quadrangle, _order=1, _generator=_structured);
Mesh mR1R2=merge(mR1P1,mR2Q1);

provides the following domains:

Domain 'g1' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 13: side 3 of element 1
 geometric side element 14: side 3 of element 3
Domain 'g3' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 17: side 3 of element 6
 geometric side element 18: side 3 of element 8
Domain 'g4' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 19: side 2 of element 1
 geometric side element 20: side 2 of element 5
Domain 's1' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 21: side 1 of element 9
 geometric side element 22: side 1 of element 10
Domain 's2' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 23: side 2 of element 10
 geometric side element 24: side 2 of element 12
Domain 's3' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 25: side 3 of element 11
 geometric side element 26: side 3 of element 12
Domain 'g2' of dimension 1 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 2 elements
 geometric side element 29: side 2 of element 4, side 4 of element 9
 geometric side element 30: side 2 of element 8, side 4 of element 11
Domain '#Omega' of dimension 2 from mesh 'P1 mesh of [0,1]x[0,1] + Q1 mesh of [1,2]x[0,1]', orientation not computed, 12 elements
 geometric element 1 (material 0, color 0, twin 0): triangle_Lagrange_1
 geometric element 2 (material 0, color 0, twin 0): triangle_Lagrange_1
 ...
 geometric element 12 (material 0, color 0, twin 0): quadrangle_Lagrange_1, orientation +0
mergeMesh

Fig. 99 Naming rules of side domains when names are distinct (left: mergeSharedBoundary=true, center: mergeSharedBoundary=false), and when there are the same on interface (right)#

Note

As it can be seen in previous examples, an element of an interface domain handles two parent elements whereas an element of a side domain handles only one parent element. For continuous-FE computations, both types of domain can be used but sometimes it may be important to work with the interface domain. Note that the member function updateParentOfSideElements() transforms a side domain into an interface domain.

Cracking mesh#

Cracking a mesh means defining an interface between two domains as a crack; the nodes located on the crack are duplicated, one belonging to a domain and the second to the other domain (see ?).

On this page