Field transfer is an essential tool in the framework of the finite elements method.
Field transfer operations are used in many cases, principally to transfer a field from one discretization (mesh) to another.
But not only, distance computation is another use of field transfer algorithms.
It enables us to (by increasing complexity):
Compute the values of shape functions at an arbitrary point in the space.
Compute the value of a (finite element) field at an arbitrary point on the space (even outside the field domain).
Compute the operator to transfer the information from a (finite element) field to a set of points.
Construct a complex transfer operator from one mesh to another (e.i. Mortar [1] operator).
Fig. 8 Example of a field transfer: Mesh and some points where the value of a field is required (by interpolation or other mean).¶
As mentioned before, the algorithm for computing the values of the shape functions of an approximation space at an arbitrary point has great value.
This basic (but technical) procedure enables us to construct more sophisticated (high-level) functions.
Here we present a selection of use cases where the field transfer has an important role.
In weakly-coupled computations, (sometimes multiphysics), the coupling loading of a problem comes from the output of a previous problem.
In many cases, the domain discretizations (mesh) used for the two problems are not compatible (even the discretized geometries can be different).
The field transfer operator is needed to transfer the information from one computation to another.
Often the solution from the simulation must be compared with real data.
The position of the sensor rarely coincides with a node of the mesh.
To be able to compare the results, an evaluation at the exact sensor location is required.
Also, the extraction of a quantity over a predefined path demands the evaluation of the field at the points of the path.
In the context of multiphysics problems several scenarios can be encountered.
In the general case, the information from one physic must be available (in some form) to the other physic to correctly take into account the coupling.
We can point out two different configurations:
When two physics are defined on the same domain (same \(\Omega\)), but with different discretization (for example to capture unknowns gradients), and the coupling is performed over the full domain.
Fig. 9 One domain (\(\Omega\)) with two discretizations.¶
When two physics are defined on different domains (\(\Omega_1\) and \(\Omega_2\)).
In this, case the information transits over the common boundary (\(\Gamma := \Omega_1 \cap \Omega_1\))
This case applies also when we have only one physic but the domain is partitioned with different discretizations.
This is the case when discretizations of different parts are imposed but a simulation of the assembly is performed.
Fig. 10 One domain (\(\Omega\)) splitted into two sub-domains (\(\Omega_1\) and \(\Omega_2\)).¶
Let \(d^\prime = 2,3\) be the dimension of the ambient space, \(d = 2,3 \leq {d^\prime}\) and let \(\Omega \subset \R^d\) be an open bounded domain in \(\R^{d^\prime}\).
Let \(\mathcal{M}\) be a conforming mesh representing the domain \(\Omega\), this can be achieved typically by employing simplicial mesh cells.
\(\mathcal{M}\) can be decomposed in \(n_e\) non-overlapping elements \(\mathcal{M_e}\).
The use of simplicial meshes is not required but makes the transfer algorithm much simpler (mapping between the parametric and the geometrical space is constant).
In the general case, conformal mixed elements meshes are often used (mix of triangles, quadrangles in 2D; tetrahedrons, hexahedrons, pyramids, wedges in 3D).
Let \(\{T_1, \ldots, T_n \}\) be a collection of \(n\) points in \(\R^{d^\prime}\) called target points (or simply \({T}_j\) for \(j\) in \([1,\dots,n]\)).
Let \(u \in V_h(\R^d)\) be the finite element field (over \(\mathcal{M}\)) to evaluate.
if \(N = \text{dim }V_h(\R^d)\), we consider a basis \((\phi_i)_{1\leq i \leq N}\) of \(V_h(\R^d)\).
The decomposition of \(u\) in the basis \(V_h(\R^d)\) is given by \(u = \sum^{N}_{i=1} u_i\phi_i\).
Note
In the case of an iso-parametric field, the approximation space used for the geometry is also used for the interpolation of the field.
In many cases, the approximation space \(V_h(\R^d)\) of the field (i.e. P0, discontinues-P1 or P2 over a P1 geometry) is not the same as the geometry one.
In this case, \(V_h(\R^d) \neq {^{geo}V_h}(\R^d)\) where \(^{geo}V_h(\R^d)\) is the approximation space of the geometry.
This makes the transfer algorithm a little bit more complex for the general case.
From Fig. 8, we can see that two following cases can be identified:
The case where the target point\(T_j\)lies inside the domain\(\mathcal{M}\)
Classically, for a given element \(e\), the shape functions are defined as a function of the barycentric coordinates \(\eta^e_k\) with \(k = [1,..,d].\)
So the first step is to find the element \(e\) containing the point \(T_j\).
This means, for a given element, the equation (1) is written in terms of the element barycentric coordinates \(\eta^e_k\),
Here \({N}^e\), \(u^e_i\) and \(\eta^e\) denote the number of shape functions in \(V_h(\R^d)\) of the element \(e\), the elementary and shape functions values extracted from the globally defined quantities, respectively.
where the inverse mapping \(\eta^e = {^{geo}\mathbf{Q}^{e}}(x)\) cannot be given explicitly in general. However, \(\eta^e\) can be computed very efficiently for any given \(\mathbf{x}\) by the means of the Newton-Raphson method.
The case where the target points lies outside of the domain\(\mathcal{M}\)
This transfer can’t be defined using the equation (3) because no element contains the target point \(T_j\).
Different strategies can be used:
Assign the value of the field at the closest node of the discretization (Nearest).
\[\widetilde{u}_j = u(x_i) = u_i\]
Here \(x_{i}\) is the position of the closest node (\(i\)) of the mesh.
Assign the value of the field at the closest position of the discretization\(\mathcal{M}\) (Clamp)
The position does not necessarily coincide with a node on the mesh.
Here \(e\) is the closest element to the target point \(T_j\).
\(\widetilde{\eta}^e\) corresponds to the barycentric coordinate inside the element \(e\), which minimizes the distance between the target point and the coordinate inside the element.
This translates to:
where \(d(a,b)\) denotes the distance between points \(a\) and \(b\).
Note
This definition is not unique in the case of discontinues spaces (e.g., P0, or discontinues-P1).
In the case where the closest point coincides with a node shared between 2 elements, only the final user can decide if this ambiguity is acceptable for the final application.
Assign the value using the shape functions of the closest element but evaluated outside the element (Extrapolation)
This is dangerous for different reasons:
The extrapolated value does not necessarily preserve field properties (i.e. a damage field is generally defined on the interval \([0,1]\) and the extrapolated values can be beyond the limits).
The notion of the closest element is not well defined (i.e. if the closest point is a node attached to multiple elements).
In this case, we search for the closest element and solve the problem (4) but without the constrain \(x(\eta^e) \in \Omega_e\)
In the case the point is outside \(\mathcal{M}\) then the value of the field is zero.
The complexity of the algorithm is linear with respect to \(N\) (the number of target points).
For each target point, we need to find the candidate elements.
This is done using a geometry search algorithm.
\(k\)-d tree [4] and R-tree [5] data structures are good choices to reduce the complexity of the search step at the expense of the tree construction.
Finally with all the shape functions values the transfer operator \(P_{N,n}\) is given by
In the context of the finite element method, this operator is sparse and can be built independently line-by-line (one line per target point).
This makes possible a parallel implementation capable of taking advantage of multi-processor computers (CPUs/GPUs).
In this section, we describe in detail the C++ Muscat implementation.
In Muscat, different user options are available to control the kind of algorithm to apply in the case where the target point \(T_j\) lies inside or outside the mesh \(M\).
Currently Muscat supports the following options:
“Interp/Nearest” : Interpolation inside, and nearest node if \(T_j\) lies outside \(M\)
“Nearest/Nearest” : Nearest node value for every \(T_j\)
“Interp/Clamp” : Interpolation inside, and the value of the orthogonal projection if \(T_j\) lies outside \(M\)
“Interp/Extrap” : Interpolation inside, and extrapolation (use of the shape function evaluated outside) if \(T_j\) lies outside \(M\)
“Interp/ZeroFill” : Interpolation inside, and no transfer (empty coefficient in \(P\)) if \(T_j\) lies outside \(M\)
A point \(T\) is in the element \(e\) if \(T \in M_e\).
Projection
A point \(T_j\) has a projection into the element \(e\) if we can find a point \(\widetilde{T}_j^e\) inside the element \(e\) such that the vector
\(\overline{T_j\widetilde{T}_j^e}\) form an orthogonal vector to the element \(e\).
Fig. 11 Two elements and four points: points \(T_2\) & \(T_3\) have projected points (\(\widetilde{T}_2^e\) & \(\widetilde{T}_3^e\)) while points \(T_1\) & \(T_4\) don’t.¶
For 3D elements, points outside of the element don’t have a projection.
For 0D element (point elements), all points have a projection.
In the next sections we present, in detail, the implementation of the transfer field operation in Muscat.
Input Field: of type Muscat.FE.Fields.FEField.FEField.
This object represent the input field.
It contains the mesh (\(M\)), the finite element space, and the associated numbering \(V_h(\R^d)\), and the values of the field (\(u_i\)) for each shape function.
Target Points: of type Muscat.Types.ArrayLike.
This object contains the target points \(T_j\), which is a simple numpy array of size (n,[2-3]).
Method: of type str.
This argument selects the method depending on the position of the point.
Element Filter: of type Muscat.Containers.Filters.FilterObjects.ElementFilter.
In many cases the user want to transfer only a partial part of the field, this filter select the elements (part of \(M\)) to be considered as the source for the transfer.
options: of type dict.
Extra options are available to select the different strategies for the potential element search and to activate the computation of the derivation operator.
Available options are:
usePointSearch (bool): to activate the search of potential elements by points.
Find the closest node of the mesh for each target point and add all the elements touching this node to the list of potential elements.
useElementSearch (bool): to activate the search of potential elements by element center.
Find the closest element center for each target point and add all the elements touching the nodes of the closest element center.
useElementSearchFast (bool): to activate the search of potential elements by element center (variant).
To make the list of potential elements smaller, the user can activate this option to only add the closest element (in the sense of the element center) and not all its neighbors.
useEdgeSearch (bool): to activate the search of potential elements by an edge search.
Find the closest edge in the mesh and add all the elements touching the two extremities to the list of potential elements.
differentialOperator (int): Activate the derivation computation.
By default, the algorithm only calculates the transfer operator (default -1).
The user can compute the derivation operator \(\frac{\partial u}{\partial x_c}\) for \(c \in [0,1,2]\).
And some less relevant, arguments:
verbose: of type (bool).
To print the progress of the computation (for big problems).
nbThreads: of type (int).
To set the number of threads to use for the transfer (by default the of threads is selected automatically).
The algorithm produces 3 outputs:
The transfer operator (6) to compute the values at each target point.
The status, which is an int numpy vector of size \(n\). For each target point, the status represents the algorithm used for the computation.
Options are:
0: “Nearest”: Information for the \(T_j\) is extracted from the position of the closest node.
1: “Interp”: Information for the \(T_j\) is extracted by interpolation (this means the point is inside \(M\)).
2: “Extrap”: Information for the \(T_j\) is extracted by extrapolation (this means the point is outside \(M\)).
3: “Clamp”: Information for the \(T_j\) is extracted from the orthogonal projection to the mesh (this means the point is outside \(M\)).
4: “ZeroFill”: No information is computed (the point is outside \(M\)).
5: “MultipleExtrapPoint”: for some points outside the mesh, the concept of closest element is not well defined. Extrapolation is done using one element.
The entities. For every target point, the entities vector holds the id of the entity (node or element) used to compute the transfer operator.
If status = 0, then entity is a node id. If status \(\in [1,2,3,5]\), the entity is an element id.
The first step of the algorithm is to build all the structures for the fast geometrical search of geometrical primitives (points or edges).
We use the boost::kdtree [4] for this job.
The three instances of kd-tree to be populated are:
nodeRTree: instance for the fast geometrical search for points.
centerRTree: instance for the fast geometrical search for the element center.
segmentRTree: instance for the fast geometrical search for edges (using the bg::model::segment concept).
We also need to compute the inverse mapping from point-to-element to correctly retrieve element ids attached to a particular node id.
Now that all the preliminary work is done (kd-trees building), we loop over the target points to compute the transfer operator.
In the case of the user chooses the “Nearest/Nearest” approach, a simplified algorithm is used because we don’t need to find any elements.
We use the nodeRTree instance to find, for every target point, the closest node on the mesh.
The construction of the output (transfer operators, status, entities) is trivial.
For all other cases we have:
for each target point \(T_j\):
Generate a list of elements that potentially contain the \(T_j\)
if option usePointSearch is active: find the closest node and add all the elements touching the node to the list of the potential elements.
if option useElementSearch is active: find the closest element center
if option useElementSearchFast is active: add the closest element to the potential elements.
if option useElementSearchFast is inactive: add all the elements touching the closest element nodes to the potential elements.
if useEdgeSearch is active: find the closes edge (bar) and add all the elements touching the nodes of the bar to the list of potential elements.
Clean repeated element ids for the potential elements list.
Sort potential elements by the distance to the center (closest first).
For each candidate element \(\mathbf{e}\):
Compute the barycentric \(\mathbf{\eta^e}\) coordinates for the point \(T_j\) (solve (3) for \(\mathbf{\eta^e}\)).
If the barycentric coordinates define a point inside the element \(e\):
Assemble line \(j\) of the transfer operator, and set correct values for the status[j] = 1, and entities[j]= e
Go to the next target point
If the barycentric coordinates define a point outside the element \(e\):
Select level 1 faces (child elements) having the normal vector pointing in the same direction as the vector from the center of the face to the target point (positive dot product).
Check if the target point has a projection into the face.
This is done using the same algorithm to compute the barycentric coordinates on an element (solve (3) for \(\mathbf{\eta^e_{L1_k}}\)).
The only difference is that \(\mathbf{\eta^e_{L1_k}}\) are the barycentric coordinates of the \(k\) level 1 child element.
if no level 1 face has a projection restart the projection algorithm with the level 2 faces.
if no level 2 face has a projection restart the projection algorithm with the level 3 faces.
With the projected point on the border (surface, line or points), we calculate the barycentric coordinate on the parent element \(e\).
Update the distance from the target point to the element \(e\).
If no element containing the target points is found:
- Use the shape functions values (for extrapolation) or the clamped shape functions (shape function from the first child element with a projection), to construct the transfer operator, and set the status and the entities outputs accordingly.
In the case of meshes with distorted elements, the basic strategy of finding the potential elements using the simple strategy of the closest node does not work properly.
Fig. 12 Examples of a 2D mesh with distorted elements¶
Fig. 12, shows a case where the strategy of the closest point fails because the closest node is not attached to the element that contains the target point.
