Finite element connectivities
Hand-made
Let’s build a very small connectivity consisting of 5 elements, 2 quadrangles Q1 and 3 triangles P1, as in the figure in comments:
const uint nElements = 5;
const uint nNodes = 8;
const uint dim = 2;
// (5)- - - - -(4)
// | \ 3 / |
// | \ / |
// | 4 \ / 2 |
// | \ / |
// (6)- -(7)- -(3)
// | | |
// | 0 | 1 |
// | | |
// (0)- -(1)- -(2)
// element number 0 1 2 3 4
auto nodesList = NaturalList { 0,1,7,6, 1,2,3,7, 3,4,7, 4,5,7, 5,6,7};
auto fePool = {FiniteElement(NV_TRIANGLE,"3P1"), FiniteElement(NV_QUADRANGLE,"4Q1")};
auto f_element_k = Natural1_FunctionN{1,1, 0,0,0}; // the first two elements are at index 1 in the pool, the next ones are at index 0
auto feFunc = Fe_FunctionN (fePool, f_element_k); // This is the list of finite elements, without nodes
evaluate("number_of_nodes_function",feFunc.number_of_nodes_function());
auto fec = FeConnectivity (feFunc, Numeric1_FunctionN__general (nodesList)); // now the connectivity is built
evaluate(fec);
number_of_nodes_function:
0 ( 4 )
1 ( 4 )
2 ( 3 )
3 ( 3 )
4 ( 3 )
reference_element_function:
0 ( 3 )
1 ( 3 )
2 ( 2 )
3 ( 2 )
4 ( 2 )
finite_element_function.pool_index_function:
0 ( 1 )
1 ( 1 )
2 ( 0 )
3 ( 0 )
4 ( 0 )
element_nodes_function:
0 ( 0 1 7 6 )
1 ( 1 2 3 7 )
2 ( 3 4 7 )
3 ( 4 5 7 )
4 ( 5 6 7 )
element_nodes_function.numericity():
0 ( 4 )
1 ( 4 )
2 ( 3 )
3 ( 3 )
4 ( 3 )
Neighbouring and reciprocal
To get the neighbours of each element, we cannot use the connectivity directly because the computation uses the vertex connectivity.
In this case, they are the same. Nethertheless, passing through a ReferenceDomain ensures this, because this object embeds
a vertex connectivity with no hole in the numbering, standing for a topology of reference elements with a given numbering of elements.
evaluate("neighbours of elements",neighbours(ReferenceDomain(fec)));
if (feFunc.number_of_nodes_function()==num(fec.element_nodes_function())) os << "correct\n";
neighbours of elements:
0 ( 1 4 )
1 ( 0 2 )
2 ( 1 3 )
3 ( 2 4 )
4 ( 0 3 )
correct
To get the reciprocal relation, that is, the elements owning a node, use the reciprocal function:
evaluate("node_elements function",reciprocal(fec.element_nodes_function()));
node_elements function:
0 ( 0 )
1 ( 0 1 )
2 ( 1 )
3 ( 1 2 )
4 ( 2 3 )
5 ( 3 4 )
6 ( 0 4 )
7 ( 0 1 2 3 4 )
To get all the nodes related to a node, use the related_values function:
evaluate("related_values function",related_values(fec.element_nodes_function()));
The node-to-node relation is given by this symmetric matrix:
: 1 2 3 4 5 6 7 8
-. - - - - - - - -.
1| x x x x|
2| x x x x x x|
3| x x x x|
4| x x x x x|
5| x x x x|
6| x x x x|
7| x x x x x|
8| x x x x x x x x|
-+ - - - - - - - -+
from which is the corresponding result:
related_values function:
0 ( 0 1 6 7 )
1 ( 0 1 2 3 6 7 )
2 ( 1 2 3 7 )
3 ( 1 2 3 4 7 )
4 ( 3 4 5 7 )
5 ( 4 5 6 7 )
6 ( 0 1 5 6 7 )
7 ( 0 1 2 3 4 5 6 7 )
Had we called the function with the upper flag set to true:
evaluate("related_values function",related_values(fec.element_nodes_function(),true));
that we would have got the upper matrix
: 1 2 3 4 5 6 7 8
-. - - - - - - - -.
1| x x x x|
2| x x x x x|
3| x x x|
4| x x x|
5| x x x|
6| x x x|
7| x x|
8| x|
-+ - - - - - - - -+
and the corresponding function:
related_values function:
0 ( 0 1 6 7 )
1 ( 1 2 3 6 7 )
2 ( 2 3 7 )
3 ( 3 4 7 )
4 ( 4 5 7 )
5 ( 5 6 7 )
6 ( 6 7 )
7 ( 7 )
We can get the list of elements sharing a node. An element being the index of the element_nodes_function, we shall write:
evaluate("related_values function",related_values(fec.element_nodes_function()));
and then obtain (yes, all the elements in our case share the central node #7):
: 1 2 3 4 5
-. - - - - -.
1| x x x x x|
2| x x x x x|
3| x x x x x|
4| x x x x x|
5| x x x x x|
-+ - - - - -+
related_indices function:
0 ( 0 1 2 3 4 )
1 ( 0 1 2 3 4 )
2 ( 0 1 2 3 4 )
3 ( 0 1 2 3 4 )
4 ( 0 1 2 3 4 )
or, with the upper matrix flag set:
evaluate("related_values function",related_values(fec.element_nodes_function(),true));
: 1 2 3 4 5
-. - - - - -.
1| x x x x x|
2| x x x x|
3| x x x|
4| x x|
5| x|
-+ - - - - -+
related_indices function:
0 ( 0 1 2 3 4 )
1 ( 1 2 3 4 )
2 ( 2 3 4 )
3 ( 3 4 )
4 ( 4 )
Transformation
From a connectivity, we can make another one with a different degree:
auto fec2 = FeConnectivity__lagrange (fec,2);
It is possible to be more precise in order to choose the degree 2 nodes: on edges, on sides, internal.
Sub-connectivity
auto fecQuad = fec( Natural1_FunctionN {0,1} );
auto fecTria = fec( Natural1_FunctionN {2,3,4} );
os << min(num(fecTria.element_nodes_function())) << "\n";
os << max(num(fecTria.element_nodes_function())) << "\n";
evaluate(fecTria);
evaluate("tria: element_nodes_function.numericity()",fecTria.element_nodes_function().numericity()());
NVI_TryProgram(num(fecQuad.element_nodes_function()).is_uniform());
NVI_TryProgram(num(fecTria.element_nodes_function()).has_only_value(3));
Reference sides
Specifying a side means giving an element index with a local side index: a side is the side of an element. We can make a list of them:
// give the reference elements on which the sides are defined, the element indices; the local side indices
auto sides1 = ReferenceSideN(fec.reference_element_function(),{1,2},{4,3}); // side #2 of element #1, side #3 of element #2
auto fecSides1 = fec(sides1);
evaluate(fecSides1)
reference_element_function - uniform in [0,1]:
( 1 )
finite_element_function.pool_index_function: - uniform in [0,1]:
( 0 )
element_nodes_function:
0 ( 6 0 )
1 ( 7 5 )
element_nodes_function.numericity() - uniform in [0,1]:
( 2 )
Let’s define the reference point being the middle of each side (the sides are segments), and then asking to sides1 for the corresponding reference point in the 2D domain.
auto f_Re = fecSides.reference_element_function();
// a point at the middle of each side of sides1; no need to provide the elements when all are taken once
auto f_sidepoint = ReferencePoint1_FunctionN(f_Re,Vector1_FunctionN(fecSides.domain_size(),Tensor::vector1({0.5})));
// the corresponding domain point
auto f_point = sides(f_sidepoint);
// list of the elements among the sides
evaluate("f_sidepoint.ie", f_sidepoint.element_index());
// list of the correponding elements of the domain
evaluate("f_point.ie", f_point.element_index());
// list of the local coordinates in the side elements
evaluate("f_sidepoint.ksi",f_sidepoint.coordinates());
// list of the corresponding local coordinates in the domain elements
evaluate("f_point.ksi",f_point.coordinates());
element_nodes_function.numericity() - uniform in [0,1]:
( 2 )
f_sidepoint.ie:
0 ( 0 )
1 ( 1 )
f_point.ie:
0 ( 0 )
1 ( 4 )
f_sidepoint.ksi - uniform in [0,1]:
[ 0.5]
f_point.ksi:
0
[ -1]
[-0.5]
1
[-0.5]
[ 1]
Rather than hand building, a more usual way is to search for some characteristical sides as:
auto domain = ReferenceDomain(fec);
auto bsides = boundary_reference_sides (domain);
auto isides1 = internal_reference_sides1(domain);
auto isides2 = internal_reference_sides2(domain);
A ReferenceSideN allows to map a side reference point to the global reference point:
auto rePointInDomain = sides(rePointInSideDomain,mirrorFlag);
The mirror flag specifies for sides appearing in two neighbour elements which element is to be taken. This way the same side reference point can be passed to the two elements and the local reference point can be computed in both elements, allowing mid value and jump computation between neighbours.
Interpolation and meshes
As we can see, the finite element connectivity does not use any node coordinate. But to be visualized as a mesh, we need to give the coordinates:
// node number: 0 1 2 3 4 5 6 7
const auto x_values = navlib::array<decimal_t> { 0,0, 1,0, 2,0, 2,1, 2,2, 0,2, 0,1, 1,1 };
// container for node coordinates, named X
auto dataX = DecimalDataN("X",nNodes*dim);
dataX.set_values(x_values);
auto X = vectorN(dataX,dim);// VectorN
auto mesh = Mesh(processor,{"MyMesh",fec,X}); // or Mesh(processor,{"MyMesh",fec,dataX}) : X will be built internally and obtained by mesh.X()
mesh.save(processor,"MyMesh.msh");
The resulting mesh displayed by Gmsh.
Let’s use the X container to make interpolations with respect to the finite element. We need as well some formal definitions of the element index \(e\) and of the list of local coordinates \((\xi)_i\), that is, the location where the interpolation occurs:
// local coordinates
auto dataKsi = DecimalDataN("ksi",dim);
auto ksi = vector1(dataKsi);
// element index
auto dataE = NaturalData1("E",{2,nElements-1});
auto ie = scalar1(dataE);/// Natural1
auto x = fec.interpolate_from_all_nodes_d<0>(X,ie,ksi); // Vector1
auto Dx = fec.interpolate_from_all_nodes_d<1>(X,ie,ksi); // Matrix1
dataE.set_value(3);
dataKsi.set_uniform_value(1/3.);
evaluate(x);
In the previous example, the reference point given to the interpolation was single. In the next example, we will use a multiple one given by the integration rule applied on the elements.
Loop on elements with integration points
Let’s see how we build a calculator to provide the measure of each element.
auto intRule = IntegrationRule(2);
auto f_Weights = integration_weights(fec,intRule); // ScalarN_FunctionN element -> list of weights
auto f_Points = integration_points (fec,intRule); // ReferencePointN_FunctionN element -> list of reference points (NB: probably, not necessarily, in the same element)
auto f_ie = f_Points.element_indices(); // ScalarN_FunctionN: i -> list of elements (elements are repeated by the number of integration points by element)
auto f_ksi = f_Points.coordinates(); // VectorN_FunctionN: i -> list of reference coordinates
auto w = f_Weights(); // ScalarN: list of weights at points (dummy loop index embedded)
auto ie = f_ie (); // ScalarN: list of elements (dummy loop index embedded)
auto ksi = f_ksi (); // VectorN: list of reference coordinates (dummy loop index embedded)
auto x = fec.interpolate_from_all_nodes_d<0>(X,ie,ksi); // VectorN x = x(ie,ksi)
auto G = fec.interpolate_from_all_nodes_d<1>(X,ie,ksi); // MatrixN G = dx/dksi
auto w = w/innersum(w); // normalisation (the length of the reference segment is 2; the measure of the reference quadrangle is 4, ... )
auto detG = det(G); // the jacobian of the real element at integration points
auto dx = w*det(G); // the *finite* integrand dx at integration point in the real elements
auto element_area = innersum(dx); // this is a numerical integration for ∫ dx on the element ie of the mesh defined by 'X' and the connectivity 'fec'
// A first use of a calculator
{
auto computer = computer(f_ie.arg_range(),element_area); // a calculator for the scalar 'element_area' embedding a dummy loop index
auto ptr_ix = computer.get<natural_t>(0); // a pointer on the memory containing the index of the loop
auto ptr_v = computer.get<decimal_t>(1); // a pointer on the memory containing the value of the first tensor (the unique one: element_area)
decimal_t total_area = 0.;
for (uint i : computer.range())
{
computer.update(i);
std::cout << i << " = " << *ptr_ix << " : " << *ptr_v << '\n';
total_area += *ptr_v;
}
}
// A more convenient way to have an automatic accumulation
{
auto scalarBuilder = ScalarBuilder{}; // a ScalarBuilder holds the formulations for a scalar
scalarBuilder.add(element_area);
auto smartScalar = SmartScalar__get("",scalarBuilder); // a SmartScalar holds the machinery for the computation specified in a ScalarBuilder
smartScalar.allocate(); // the allocation is not done before the user's request
smartScalar.update(); // the computation is done
auto total_area = smartScalar.value();
std::cout << "total mesh area: " << total_area << '\n';
}
If X depends on a variable as is often the case - simply consider the coordinates of the mesh being stored in a variable container -
then we can reuse the computer once the container values have changed and get the updated values.
...
while (< X was updated>)
{
auto total_area = smartScalar.updated_value();
std::cout << "total mesh area: " << total_area << '\n';
}
If we want to get the updated value of the area of each element, a SmartVector is more convenient.
...
{
auto vectorBuilder = VectorBuilder{}; // a VectorBuilder holds the formulations for a vector
vectorBuilder.add(ie,element_area);
auto smartVector = SmartVector__get("",vectorBuilder); // a SmartScalar holds the machinery for the computation specified in a ScalarBuilder
smartVector.allocate(); // the allocation is not done before the user's request
smartVector.update(); // the computation is done
auto values = navlib::array<decimal_t>(smartVector.size());
smartVector.save_to(values); // the SmartVector internal storage structure is not known
}