Geometric domains#

Material Id#

The materialId property is managed by the functions:

• GeomDomain::setMaterialId that sets the material id of all elements of a domain,

• getMaterialId that gets the material id of the current element in loops of FE calculation.

Real rho(const Point&P, Parameters& pars=defaultParameters)
{ Number mat=getMaterialId();
  if(mat==1) return 1.;
  return 2.;}
...
Domain omega1=mesh2d.domain("omega_1"); // subdomain 1
Domain omega2=mesh2d.domain("omega_2"); // subdomain 2
omega1.setMaterialId(1);  omega2.setMaterialId(2);

Color#

The color property is used to bypass computations in FE loops if the colorFlag is on; more precisely when the color of an element is not 0 there is no FE computation on this element. This allows to compute integrals only on a subset of the elements that change during the calculation (e.g level set method). To change dynamically the color of elements of a domain, XLiFE++ use the setColor function with the following prototype:

void setColor(const GeomDomain& dom, const TermVector& val, ColoringRule cr)
void setColor(const GeomDomain& dom, const TermVector& val, VectorColoringRule cr)

with val any TermVector containing some values or vector values used by a ColoringRule or a VectorColoringRule that are C++ functions of the form:

Real coloringRule(const GeomElement& gelt, const vector<Real>& val);
Real vectorColoringRule(const GeomElement& gelt, const vector<Vector<Real>>& val);

where gelt will be the current geometric element and val values or vector values at element vertices.

A classical usage of color is the following:

Real myColoringRule(const GeomElement& gelt, const vector<Real>& val)
{ res=0;...; return res;}            // coloring gelt regarding val
...
Domain omega=mesh2d.domain("omg");   // shortcut to domain
omega.setColorFlag(true);            // activate color usage
TermVector tv=...;                   // compute a TermVector on omega
setColor(omega,tv,myColoringRule);   // set element colors using tv and myColoringRule

Note that the color management is not on by default. Use setColorFlag to switch on/off.

Sides domain#

Discontinuous Galerkin methods involve the set of sides (edges/faces) of elements of a plain domain. XLiFE++ deals with such sets: either the set of all element sides or the set of all internal element sides. These domains are either created by the Domain:: sidesDomain member function or by the sides function taking a domain as input argument. The new domain is allocated in the parent domain and the previous functions return it, if it has been already built.

Rectangle R(_origin=Point(0.,0.),_xlength=1., _ylength=1.,_nnodes=11,
            _domain_name="Omega",_side_names="Gamma");
Mesh mR(R, _shape=_triangle, _order=1, _generator=_structured);
Domain omega=mR.domain("Omega"), gamma=mR.domain("Gamma");
Domain somega=sides(omega);            // all side elements (edges) of omega
Domain isomega=internalSides(omega);   // all internal side elements of omega

Domain isomega=somega-gamma;

Sides domain may be used as the other domains for integrals involved in bilinear forms.

Extended domain#

The calculation of FE integrals on a side domain $\mathrm{\Gamma }$ of $\mathrm{\Omega }$ that involves non tangential derivatives has to take into account degrees of freedom that are inside $\mathrm{\Omega }$. To avoid to take into account all the elements and degrees of freedom of $\mathrm{\Omega }$, XLiFE++ is able to construct the domain containing the elements of $\mathrm{\Omega }$ having a side (edge or face) located on $\mathrm{\Gamma }$. This is the extended domain by sides:

Sometimes some FE data only defined on a side domain $\mathrm{\Gamma }$ of $\mathrm{\Omega }$ must be extended locally to the interior of $\mathrm{\Omega }$. In that case XLiFE++ builds the domain containing the elements of $\mathrm{\Omega }$ having a vertex located on $\mathrm{\Gamma }$. This is the extended domain by vertices:

To build an extended domain, use the Domain:: extendDomain with a boolean argument specifying the type of extension (true for vertex extension and false for side extension that is the default):

Rectangle R(_origin=Point(0.,0.),_xlength=1., _ylength=1.,_nnodes=6,
            _domain_name="Omega",_side_names={"Gamma","","Sigma",""});
Mesh mR(R, _shape=_triangle, _order=1);   saveToFile("mR", mR);
Domain Omega=mR.domain("Omega"), Gamma=mR.domain("Gamma"), Sigma=mR.domain("Sigma");
Domain sigma_exts=Sigma.extendDomain(false);   saveToFile("sigma_exts", sigma_exts);
Domain sigma_extv=Sigma.extendDomain(true);    saveToFile("sigma_extv", sigma_extv);

Regarding the extension type, the new extended domain is named side_domain_name_sextension or side_domain_name_vextension. Obviously, it can be renamed.

When extend an interface $\mathrm{\Gamma }$ between the domains ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$, it is possible to specify if the exetension concerns one of the domains ${\mathrm{\Omega }}_1$ or ${\mathrm{\Omega }}_2$ or both of them.

Rectangle R1(_origin=Point(0.,0.),_xlength=1., _ylength=1.,_nnodes=6, _domain_name="Omega1",
             _side_names={"","Sigma","",""});
Rectangle R2(_origin=Point(1.,0.),_xlength=1., _ylength=1.,_nnodes=6, _domain_name="Omega2",
            _side_names={"","","","Sigma"});
Mesh mR12(R1+R2,_shape=_triangle,_order=1,_generator=gmsh);
Domain omg1=mR12.domain("Omega1"), omg2=mR12.domain("Omega2"), sigma12=mR12.domain("Sigma");
sigma12.updateParentOfSideElements();
Domain int12_s=sigma12.extendDomain(false); saveToFile("int12_s",int12_s);
Domain int12_v=sigma12.extendDomain(true);  saveToFile("int12_v",int12_v);
Domain int12_s1=sigma12.extendDomain(false, omg1); saveToFile("int12_s1",int12_s1);
Domain int12_v1=sigma12.extendDomain(true, omg1);  saveToFile("int12_v1",int12_v1);

Fictitious domain#

A fictitious domain corresponds to the intersection of domains. More precisely, let $\mathrm{\Gamma }=\bigcup F_j$ and $\mathrm{\Omega }=\bigcup E_i$ two domains, fictitious domain is defined as follows:

The building process of a ficticious domain build an additionnal map that associates to an element of $\mathrm{\Gamma }$ the list of elements of $\mathrm{\Omega }$ intersecting it. Generally, $\mathrm{\Gamma }$ is of dimension n-1 when $\mathrm{\Omega }$ is of dimension n (curve in surface or surface in volume). This is useful to compute integrals involving unknown/test function on $\mathrm{\Gamma }$ and other unknown/test function on $\mathrm{\Omega }$. But, it also works with domains of same dimension; that allows XLiFE++ to deal with domain constructed on different meshes!

The name of the fictitious domain created is gammaName_F_omegaName.

The following code builds the fictitious domain of a circle domain immersed in a rectangle domain:

Rectangle R(_origin=Point(0.,0.),_xlength=1., _ylength=1.,_nnodes=11, _domain_name="R");
Disk D(_center=Point(0.5,0.5),_radius=0.25, _nnodes=10, _domain_name="D");
Mesh mR(R,_shape=_triangle,_order=1,_generator=gmsh);   // mesh of rectangle
Mesh mC(D,_shape=_segment,_order=1,_generator=gmsh);    // mesh of **circle**
Domain omega=mR.domain("R"),                 // 2d domain
       gamma=mC.domain("D");                 // 1d domain
Domain gammaF=gamma.ficticiousDomain(omega); // fictitious domain named "D_F_R"
saveToFile("omega",omega); saveToFile("gamma",gamma); saveToFile("gammaF",gammaF);

The next code (very similar) builds the fictitious domain of a disk domain immersed in a rectangle domain:

Rectangle R(_origin=Point(0.,0.),_xlength=1., _ylength=1.,_nnodes=11, _domain_name="R");
Disk Di(_center=Point(0.5,0.5),_radius=0.25, _nnodes=10, _domain_name="D");
Mesh mR(R,_shape=_triangle,_order=1,_generator=gmsh);    // mesh of rectangle
Mesh mD(Di,_shape=_triangle,_order=1,_generator=gmsh);   // ** mesh of disk **
Domain omega=mR.domain("R"),               // 2d domain
       D=mD.domain("D");                   // 2d domain
Domain Df=D.ficticiousDomain(omega);       // fictitious domain named "D_F_R"

If there exists a map $\phi :\mathrm{\Gamma }\to \mathrm{\Omega }$, that is a DomainMap object (see Defining map between two domains), it will be used to build the fictitious domain as follows:

Defining map between two domains#

DomainMap is a class handles a map from one domain (dom1) to an other one (dom2):

class DomainMap
{protected:
 const GeomDomain* dom1_p; // start domain
 const GeomDomain* dom2_p; // end domain
 Function map1to2_;        // map function from dom1 to dom2
 public:
  bool useNearest;          // use nearest method instead of locate method after mapping
  string_t name;            // a unique name
}

The Function objet map1to2_ must refer to a c++ function of the form:

Reals funMap(const Point& p, Parameters& pa = defaultParameters);

Domain maps can be used in different context: periodicity condition, transport of a function from one domain to an other domain (spectral basis for instance, see ?). Sometimes, the transporting point has to be localized in the arrival domain using localization function that proposed two methods: an exact localization or an approximate localization using the nearest point or element. This is the purpose of the DomainMapuseNearest boolean variable.

The DomainMap holds the list of all domain maps so it’s easy to find a map if it’s been defined. To define a domain map, the simplest way consists in using the defineMap function:

Reals transxm1(const Point& p, Parameters& pa = defaultParameters)
{ Reals q={p(1)-1,p(2)}; return q;}  // translation (x,y) -> (x-1,y))
...
Rectangle R(_origin=Point(0.,0.),_xlength=1., _ylength=1.,_nnodes=6,
            _domain_name="Omega",_side_names={"","x=1","","x=0"});
Mesh mR(R,_shape=triangle);
Domain omega=mR.domain("Omega"), sigma=mR.domain("x=1"), gamma=mR.domain("x=0");
defineMap(sigma,gamma,transxm1);  // define DomainMap with useNearest=false

Other posibilities are

defineMap(sigma,gamma,transxm1,true);           // define with useNearest=true
defineMap(sigma,gamma,transxm1,"sigma->gamma"); // define with a name,useNearest=false
defineMap(sigma,gamma,transxm1,"sigma->gamma",true);// define with a name,useNearest=true

When no name passed, the name of the domain map is “dom1_name–>dom2_name”. Name has to be unique!

Some processes (e.g. processing periodic conditions between two boundaries) requires the map related one boundary (dom1) to the other one (dom2). XLiFE++ uses the findMap function, that’s either find a map related dom1 to dom2 if it has been defined or try to build one (only translation).

Merging domains#

XLiFE++ can merge domains coming from the same mesh and having the same dimension. The process creates a new DomainMap object by

• merging the list of geometric elements (that are supposed to non overlap)

• merging the list of nodes (duplicate nodes are not reduced)

To merge domains use the merge functions:

Domain merge(Domain dom1, Domain dom2, Domain dom3, ..., String name);

The number of domains can be up 20.

The following code illustrates the merging of two 2d domains:

Rectangle R(_origin=Point(0.,0.),_xlength=1., _ylength=1.,_nnodes=5,
           _domain_name="Omega0");
Disk Dp(_center=Point(1.,0.5),_radius=0.5, _angle1=-pi_/2,_angle2=pi_/2,_nnodes=3,
        _domain_name="Omega+");
Disk Dm(_center=Point(0.,0.5),_radius=0.5, _angle1=pi_/2,_angle2=3*pi_/2,_nnodes=3,
        _domain_name="Omega-");
Mesh mRD(R+Dm+Dp,_shape=triangle);
Domain omg0=mRD.domain("Omega0"), omgp=mRD.domain("Omega+"), omgm=mRD.domain("Omega-"),
omega=mRD.domain("#Omega");      // full mesh domain
Domain omg0p=merge(omg0,omgp,"omg0p");
saveToFile("omg0",omg0); saveToFile("omgp",omgp); saveToFile("omgm",omgm);
saveToFile("omg0p",omg0p); saveToFile("omega",omega);

Note that nodes on interface are duplicated in the original mesh done by Gmsh and remain duplicated in the merging process!

In an alternative way, the operator + may be used to merge two domains:

Domain omg0p = omg0 + omgp;   // named "0mega0 + Omega+"

Domains difference#

The difference between two domains, say ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ of same dimension, is the domain defined as

There are two ways to get the new domain $\mathrm{\Delta }$:

Domain delta = omega1.removeElts(omega2, [name]);
Domain delta =  omega1-omega2;  // named "omega1 - omega2"

Using domains defined in merging example, the code

Domain omg2=omega-omgp; saveToFile("omg2",omg2);

gives the same domains as omg0 + omgm:

Updating domains#

In some numerical methods, geometric domains can changed and have to be redefined during an iterative process (think about level set methods for instance). XLiFE++ provides a general tool to do this : the rebuild function with the following prototypes, rebuiding 2 domains and optionnally 1 or 2 side domains:

void rebuild(dom1, comparisonFun1, dom2, comparisonFun2);
void rebuild(dom1, comparisonFun1, dom2, comparisonFun2, sdom);
void rebuild(dom1, comparisonFun1, dom2, comparisonFun2, sdom1, sdom2);

The comparisonFun are objects of the ComparisonFunction class which handles boolean criteria involving either _color or _index element property, for instance:

ComparisonFunction cr1 = (_color == 0) && (_index>0),
                   cr2 = (_color>0);

Obviously, the rebuild function works in conjunction with the setColor function. How does it work?

• the dom1 is rebuilt from elements for which the comparisonFun1 is true (e.g. elements with their color property set to 0),

• the dom2 is rebuilt from elements for which the comparisonFun2 is true (e.g. elements with their color property positive),

• when side domains are specified, they are updated if dom1 or dom2 have changed.

There exists also a version of the rebuild function that works with lists of domains, of comparison functions and of side domains:

void rebuild(vector<GeomDomain*>& doms, vector<ComparisonFunction<> >& crs,
             set<GeomDomain*>& sidedoms);

The next example shows how to rebuild two domains ${\mathrm{\Omega }}_1$, ${\mathrm{\Omega }}_2$, using the coloring rule on an element ($x_i$ abcissa of the i-th vertex of element):

Real colorRule(const GeomElement& gelt, const std::vector<Real>& val)
{ Real mx=*std::max_element(val.begin(),val.end());
  if (mx>=1.5) return 1; else return 0;}
...
Rectangle R1(_origin={0.,0.},_xlength=1.,_ylength=1.,_hsteps=0.1,_domain_name="Omega_1");
Rectangle R2(_origin={1.,0.},_xlength=1.,_ylength=1.,_hsteps=0.1,_domain_name="Omega_2");
Mesh mR12(R1+R2,_shape=triangle);
Domain omega=mR12.domain("Omega"),omega1=mR12.domain("Omega_1"),
                                  omega2=mR12.domain("Omega_2");
Space V(_domain=omega,_interpolation=P1,_name="V"); Unknown u(V,_name="u");
TermVector tv(u,omega,_x);      // vector on mesh vertices containing x_i
setColor(omega,tv,colorRule);   // compute color according to the colorRule
ComparisonFunction<> cr1= _color==0., cr2= _color>0.; // rebuild criteria
rebuild(omega1,cr1,omega2,cr2); // rebuild omega1 (elements with color = 0)
                                // and     omega2 (elements with color > 0)

for(auto gp : omega.meshDomain()->geomElements)  // loop on elements of omega
   {gp->color= ...}

for(Number k=1; k<=omega.numberOfElements(); k++)  // loop on elements of omega
   {omega.element(k)->color= ...}

Partitioning a domain#

The partition of a domain is a set of non-overlapping sub-domains covering the original domain $\mathrm{\Omega }$:

Here $\mathrm{\Omega }$, ${\mathrm{\Omega }}_i$ are considered as sets of geometric elements of same dimension.

XLiFE++ provides Domain::partition member function to create a partition of a domain:

dom.partition(par1,[par2], [par3],...);

where par1, par2, … (up to 10) are partitioning parameters passed using key/value:

• _nbParts : number of partitions (2 by default)

• _partmesh: type of partition: dual (default) or nodal

• _ncommon : number of common nodes that two elements must have in order to put an edge between them in the dual graph (default 1).

• _ptype : partitionning method: rec (recursive bisectioning) or kway (default)

• _ctype : matching scheme to be used during coarsening: rm (random matching) or shem (sorted heavy-edge matching, default)

• _iptype : scheme to compute the initial partitioning (only when _ptype= rec): grow (greedy strategy, default) or iprandom (random)

• _objtype : partitioning’s objective function (only when _ptype= kway): cut (edge cut minimization, default) or vol (total communication volume minimization)

• _rtype : refinement algorithm: fm (FM-based cut refinement, default), greedy (greedy strategy), sep2sided or sep1sided

When executing the Domain::partition function on a domain, say omega object , a PartitionData object is allocated in omega (partitionData_p pointer) and Domain::partitionData function gives an access to this object. This object provides many informations on the partition:

The mesh sub-domains of the partition are available using Mesh::domain function with the name: omegaName_i. The interface between sub-domains are also created as Domain objects and they are named: D_omegaName_i_omegaName_j (i<j).

The following example illustrates how Domain::partition function works (5 sub-domains and default patitioning parameters):

Rectangle Rp(_origin=Point(0.,0.),_xlength=1., _ylength=1.,_hsteps=0.1, _domain_name="Omega");
Disk Di(_center=Point(0.5,0.5),_radius=0.25,_hsteps=0.1);
Mesh mRpDi(Rp-Di,_shape=triangle);
Domain omega=mRpDi.domain("Omega");
omega.partition(_nbParts=5);
saveToFile("mRpDi.msh",mRpDi);  // to visualize partitions with gmsh

The same example but using _ncommon=2, that makes interfaces more regular:

omega.partition(_nbParts=5, _ncommon=2);