GCC Code Coverage Report

Directory: ./
File: DegreeOfFreedom/supportDofMesh.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 447 776 57.6%
Branches: 342 1108 30.9%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: J. Foulon and others
13 //
14
15 /*!
16 \file supportDofMesh->cpp
17 \author J. Foulon V. Martin and others
18 \date 07/10/2010
19 \brief File contains SupportDofMesh implementation
20 */
21
22 // System includes
23
24 // External includes
25
26 // Project includes
27 #include "DegreeOfFreedom/supportDofMesh.hpp"
28 #include "FiniteElement/refElement.hpp"
29 #include "FiniteElement/curvilinearFiniteElement.hpp"
30 #include "FiniteElement/currentFiniteElement.hpp"
31
32 namespace felisce
33 {
34 SupportDofMesh::StringToReferenceElement SupportDofMesh::m_eltRefNameToRefEle;
35
36 769 void SupportDofMesh::m_initMap()
37 {
38
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementNode" ] = & refElementNode;
39
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementSegmentP1" ] = & refElementSegmentP1;
40
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementSegmentP1b" ] = & refElementSegmentP1b;
41
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementSegmentP2" ] = & refElementSegmentP2;
42
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementTriangleP1" ] = & refElementTriangleP1;
43
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementTriangleP1b" ] = & refElementTriangleP1b;
44
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementTriangleP2" ] = & refElementTriangleP2;
45
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementQuadrangleQ1" ] = & refElementQuadrangleQ1;
46
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementQuadrangleP1xP2" ] = & refElementQuadrangleP1xP2;
47
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementQuadrangleQ1b" ] = & refElementQuadrangleQ1b;
48
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementQuadrangleQ2" ] = & refElementQuadrangleQ2;
49
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementQuadrangleQ2c" ] = & refElementQuadrangleQ2c;
50
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementTetrahedronP1" ] = & refElementTetrahedronP1;
51
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementTetrahedronP1b" ] = & refElementTetrahedronP1b;
52
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementTetrahedronP2" ] = & refElementTetrahedronP2;
53
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementHexahedronQ1" ] = & refElementHexahedronQ1;
54
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementHexahedronQ1b" ] = & refElementHexahedronQ1b;
55
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementHexahedronQ2" ] = & refElementHexahedronQ2;
56
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementHexahedronQ2c" ] = & refElementHexahedronQ2c;
57
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementPrismR1" ] = & refElementPrismR1;
58
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementPrismP1xP2" ] = & refElementPrismP1xP2;
59
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementPrismR2" ] = & refElementPrismR2;
60
2/4
✓ Branch 2 taken 769 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 769 times.
✗ Branch 6 not taken.
 769 m_eltRefNameToRefEle[ "refElementTetrahedronRT0" ] = & refElementTetrahedronRT0;
61 769 }
62
63 /*!
64 \brief Constructor for the global SupportDofMesh
65 \param[in] variable Variable associated to this SupportDofMesh
66 \param[in] mesh Reference to the mesh to etablish a connexion with finite element
67 geometric and finite element reference
68 */
69 769 SupportDofMesh::SupportDofMesh(const Variable& variable, const GeometricMeshRegion::Pointer& mesh):
70 769 m_mesh(mesh),
71
1/2
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
 769 m_variable(variable),
72 769 m_copiedMeshPoint(false),
73 769 m_numDofSupportedByEdge(0),
74 769 m_resizedEdgeNode(false),
75 769 m_resizedFaceNode(false),
76 769 m_resizedVolumeNode(false),
77 769 m_createdEmbeddedInterfaceMaps(false),
78 769 m_createdCrackInterfaceMaps(false)
79 {
80
1/2
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
 769 m_initMap();
81 //=============================
82 // initialize the mesh/supportDof : copy points, resize m_listNode for internal edges and faces if necessary
83 769 const std::vector<ElementType>& bagElementTypeDomain = m_mesh->bagElementTypeDomain();
84
2/2
✓ Branch 1 taken 769 times.
✓ Branch 2 taken 769 times.
 1538 for (std::size_t ieltype = 0; ieltype < bagElementTypeDomain.size(); ++ieltype) {
85 769 const ElementType eltType = bagElementTypeDomain[ieltype];
86
87 //! definition of the finite element type ( geoEle, RefEle )
88 //! allow to pass from geoEleP1 to refEleP2
89 769 const GeoElement& geoEle = *GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
90
1/2
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
 769 const RefElement& refEle = m_chooseRefEleFromFEType(geoEle, eltType);
91
92 769 const int numDOFNodeVertex = refEle.numDOFNodeVertex();
93 769 const int numDOFEdge = refEle.numDOFEdge();
94 769 const int numDOFFace = refEle.numDOFFace();
95 769 const int numDOFVolume = refEle.numDOFVolume();
96
97 //! if at least one element is vertex based: copy the nodes:
98
1/2
✓ Branch 0 taken 769 times.
✗ Branch 1 not taken.
 769 if ( numDOFNodeVertex > 0 )
99
1/2
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
 769 m_copyNodes();
100
101
2/2
✓ Branch 0 taken 168 times.
✓ Branch 1 taken 601 times.
 769 if ( numDOFEdge > 0) {
102
2/2
✓ Branch 2 taken 100 times.
✓ Branch 3 taken 68 times.
 168 if ( refEle.numDOF() == geoEle.numPoint() ) m_mesh->statusEdges() = true;
103
104
2/2
✓ Branch 2 taken 68 times.
✓ Branch 3 taken 100 times.
 168 if ( m_mesh->statusEdges() == false ) {
105
1/2
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
 68 std::stringstream msg;
106
2/4
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 68 times.
✗ Branch 5 not taken.
 68 msg << "-->building internal edges" << std::endl;
107
1/2
✓ Branch 2 taken 68 times.
✗ Branch 3 not taken.
 68 m_mesh->buildEdges();
108
2/4
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 68 times.
✗ Branch 5 not taken.
 68 msg << "<--edges built" << std::endl;
109
2/4
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 68 times.
✗ Branch 6 not taken.
 68 PetscPrintf(PETSC_COMM_WORLD,"%s",msg.str().c_str());
110 68 }
111
2/2
✓ Branch 1 taken 108 times.
✓ Branch 2 taken 60 times.
 168 if (refEle.edgesCanHaveDifferentNumberOfDof())
112
1/2
✓ Branch 3 taken 108 times.
✗ Branch 4 not taken.
 108 m_numDofSupportedByEdge = m_mesh->listEdges().countAndComputeIdEdgesSupportingDofs(refEle);
113 else
114
1/2
✓ Branch 3 taken 60 times.
✗ Branch 4 not taken.
 60 m_numDofSupportedByEdge = m_mesh->numEdges() * refEle.numDOFPerEdge();
115
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
 168 m_resizeEdgeNodes();
116 }
117
118
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 769 times.
 769 if ( numDOFFace > 0 ) {
119 ✗ if ( refEle.numDOF() == geoEle.numPoint() ) m_mesh->statusFaces() = true;
120
121 ✗ if ( m_mesh->statusFaces() == false ) {
122 ✗ std::cout << "-->building internal faces" << std::endl;
123 ✗ m_mesh->buildFaces();
124 ✗ std::cout << "<--faces built" << std::endl;
125 }
126 ✗ m_resizeFaceNodes(refEle);
127 }
128
129 //! if at least one element is volume based: build the volumes
130
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 769 times.
 769 if( numDOFVolume > 0)
131 ✗ m_resizeVolumeNodes(refEle);
132 }
133
134 769 felInt cptElt = 0;
135 felInt countEltPerType[ GeometricMeshRegion::m_numTypesOfElement ];
136
2/2
✓ Branch 0 taken 19225 times.
✓ Branch 1 taken 769 times.
 19994 for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ )
137 19225 countEltPerType[ityp] = 0;
138
139 //! build the edge and face nodes
140 //! Domain
141 769 bool boundary_flag = false;
142
1/2
✓ Branch 3 taken 769 times.
✗ Branch 4 not taken.
 769 m_buildNodesOfEdgeFaceVolumePerBag(m_mesh->bagElementTypeDomain(), countEltPerType, boundary_flag);
143
144 //! Boundary (it should not be useful to build the edge/face nodes: already done via the domain!!...?)
145 769 boundary_flag = true;
146
1/2
✓ Branch 3 taken 769 times.
✗ Branch 4 not taken.
 769 m_buildNodesOfEdgeFaceVolumePerBag(m_mesh->bagElementTypeDomainBoundary(), countEltPerType, boundary_flag);
147 // Do it for all bags (3d, 2d, 1d??)
148
149 // Set m_numSupportDof
150 769 m_numSupportDof=m_listNode.size();
151
152
3/4
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 20 times.
✓ Branch 4 taken 749 times.
 769 if(FelisceParam::instance().FusionDof) { //For periodic boundary conditions
153
1/2
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 20 m_listPerm.resize(m_listNode.size(), 0);
154
2/2
✓ Branch 1 taken 15360 times.
✓ Branch 2 taken 20 times.
 15380 for (std::size_t lp = 0; lp<m_listNode.size(); lp++ ) {
155 15360 m_listPerm[lp] = lp ;//this is just the list of node if variable not fusioned
156 }
157
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 fusionSupportDof(variable);
158 }
159
160 //! Compute the size of support vectors:
161 //! m_iSupportDof (sum of the number of Dof in each elements)
162 //! m_iEle (Total number of elements + 1)
163 //! m_vectorIdElementSupport (Total number of elements)
164
1/2
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
 769 m_resizeSupportVectors();
165
166 // reset the counter
167
2/2
✓ Branch 0 taken 19225 times.
✓ Branch 1 taken 769 times.
 19994 for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ )
168 19225 countEltPerType[ityp] = 0;
169
170 //! compute support vectors in CSR format: m_iSupportDof, m_iEle, and compute m_vectorIdElementSupport
171
1/2
✓ Branch 3 taken 769 times.
✗ Branch 4 not taken.
 769 m_buildSupportCSR(m_mesh->bagElementTypeDomain(), cptElt, countEltPerType);
172
1/2
✓ Branch 3 taken 769 times.
✗ Branch 4 not taken.
 769 m_buildSupportCSR(m_mesh->bagElementTypeDomainBoundary(), cptElt, countEltPerType);
173
174 // Embedded Interface
175
2/4
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 769 times.
 769 if(FelisceParam::instance().EmbeddedInterface)
176 ✗ matchEmbeddedLabelPairs(variable);
177
178 // Cracks
179
2/4
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 769 times.
 769 if(FelisceParam::instance().hasCrack){
180 ✗ matchCrackLabelPairs(variable);
181 }
182 769 }
183
184 /*!
185 \brief Constructor of local SupportDofMesh (partitionning of SupportDofMesh in n locals SupportDofMesh).
186 \param[in] variable Variable associate to this SupportDofmesh
187 \param[in] meshLocal Reference to the local mesh.
188 \param[in,out] loc2GlobElemSupport Local to global mapping of element/element support ids.
189 \param[in] globSupportDofMesh Global SupportDofmesh
190 \param[in] loc2GlobSupport Mapping on global and local numbering of Support.
191 \param[in] rank Rank of the working process.
192 */
193 762 SupportDofMesh::SupportDofMesh(const Variable& variable, const GeometricMeshRegion::Pointer& meshLocal, std::vector<felInt>& loc2GlobElemSupport, SupportDofMesh& globSupportDofMesh):
194 762 m_mesh(meshLocal),
195
1/2
✓ Branch 1 taken 762 times.
✗ Branch 2 not taken.
 762 m_variable(variable),
196 762 m_numDofSupportedByEdge(0),
197 762 m_createdEmbeddedInterfaceMaps(false),
198 762 m_createdCrackInterfaceMaps(false)
199 {
200 // Copy all the support dof from the global mesh
201
1/2
✓ Branch 2 taken 762 times.
✗ Branch 3 not taken.
 762 m_listNode = globSupportDofMesh.listNode();
202
203 // Initialize variables
204 ElementType eltType; // Element type
205 762 std::size_t cptElt = 0; // Number of local support element
206 762 felInt cptSupportElem = 0; // Count the number of support element
207 762 felInt cptSupport = 0; // Number of local support dof
208 762 felInt numSupportDofElem = 0; // Number of support dof on an element
209 762 felInt idGlobalElementSupport = 0; // Global id of a support element
210
211 // Local bags of element type
212 762 const std::vector<ElementType>& bagElementTypeDomain = m_mesh->bagElementTypeDomain();
213 762 const std::vector<ElementType>& bagElementTypeDomainBoundary = m_mesh->bagElementTypeDomainBoundary();
214
215 // Total number of local elements
216 762 cptElt = loc2GlobElemSupport.size();
217
218 // m_link between mesh element and element support
219
1/2
✓ Branch 1 taken 762 times.
✗ Branch 2 not taken.
 762 m_vectorIdElementSupport.resize(cptElt);
220
221 // Total number of support dof per local elements
222
2/2
✓ Branch 0 taken 1452334 times.
✓ Branch 1 taken 762 times.
 1453096 for(std::size_t iel=0; iel<cptElt; iel++) {
223
1/2
✓ Branch 3 taken 1452334 times.
✗ Branch 4 not taken.
 1452334 globSupportDofMesh.getIdElementSupport(loc2GlobElemSupport[iel], m_vectorIdElementSupport[iel]);
224
2/2
✓ Branch 2 taken 1453304 times.
✓ Branch 3 taken 1452334 times.
 2905638 for(std::size_t iSup=0; iSup < m_vectorIdElementSupport[iel].size(); ++iSup) {
225 1453304 cptSupport += globSupportDofMesh.getNumSupportDof(m_vectorIdElementSupport[iel][iSup]);
226 1453304 ++cptSupportElem;
227 }
228 }
229
230 // Allocate the arrays
231
1/2
✓ Branch 1 taken 762 times.
✗ Branch 2 not taken.
 762 m_iEle.resize(cptSupportElem+1, 0); // m_iEle array of the CSR format: the support elements
232
1/2
✓ Branch 1 taken 762 times.
✗ Branch 2 not taken.
 762 m_iSupportDof.resize(cptSupport, 0); // m_iSupportDof of the CSR format: the support dof
233
1/2
✓ Branch 1 taken 762 times.
✗ Branch 2 not taken.
 762 loc2GlobElemSupport.resize(cptSupportElem, 0);
234
235 // Filling of m_iEle and m_iSupportDof and m_vectorIdElementSupport arrays
236 762 m_iEle[0] = 0;
237 762 cptElt = 0;
238 762 cptSupportElem = 0;
239
240 // for each element type in the domain
241
2/2
✓ Branch 1 taken 741 times.
✓ Branch 2 taken 762 times.
 1503 for (std::size_t i=0; i < bagElementTypeDomain.size(); ++i) {
242 741 eltType = bagElementTypeDomain[i];
243
244 // for each element of the current element type
245
2/2
✓ Branch 2 taken 1310439 times.
✓ Branch 3 taken 741 times.
 1311180 for ( felInt iel = 0; iel < m_mesh->numElements(eltType); iel++ ) {
246
2/2
✓ Branch 2 taken 1311357 times.
✓ Branch 3 taken 1310439 times.
 2621796 for( std::size_t iSup=0; iSup < m_vectorIdElementSupport[cptElt].size(); ++iSup) {
247 1311357 idGlobalElementSupport = m_vectorIdElementSupport[cptElt][iSup];
248 1311357 numSupportDofElem = globSupportDofMesh.getNumSupportDof(idGlobalElementSupport);
249
250 // fill m_iEle for this element (CSR format)
251 1311357 m_iEle[cptSupportElem+1] = m_iEle[cptSupportElem] + numSupportDofElem;
252
253 // fill m_iSupportDof (add the support dof of this element) (CSR format)
254
2/2
✓ Branch 0 taken 5136168 times.
✓ Branch 1 taken 1311357 times.
 6447525 for ( int iSupport = 0; iSupport < numSupportDofElem; iSupport++)
255 5136168 m_iSupportDof[m_iEle[cptSupportElem] + iSupport] = globSupportDofMesh.iSupportDof()[globSupportDofMesh.iEle()[idGlobalElementSupport] + iSupport];
256
257 // fill m_vectorIdElementSupport
258 1311357 m_vectorIdElementSupport[cptElt][iSup] = cptSupportElem;
259
260 // fill loc2GlobElemSupport
261 1311357 loc2GlobElemSupport[cptSupportElem] = idGlobalElementSupport;
262
263 // increment the count of local support elements
264 1311357 ++cptSupportElem;
265 }
266
267 // increment the count of local elements
268 1310439 ++cptElt;
269 }
270 }
271
272 // same for the boundary domain elements
273
2/2
✓ Branch 1 taken 798 times.
✓ Branch 2 taken 762 times.
 1560 for (std::size_t i=0; i < bagElementTypeDomainBoundary.size(); ++i) {
274 798 eltType = bagElementTypeDomainBoundary[i];
275
276 // for each element of the current element type
277
2/2
✓ Branch 2 taken 141895 times.
✓ Branch 3 taken 798 times.
 142693 for (felInt iel = 0; iel < m_mesh->numElements(eltType); iel++) {
278
2/2
✓ Branch 2 taken 141947 times.
✓ Branch 3 taken 141895 times.
 283842 for(std::size_t iSup=0; iSup < m_vectorIdElementSupport[cptElt].size(); ++iSup) {
279 141947 idGlobalElementSupport = m_vectorIdElementSupport[cptElt][iSup];
280 141947 numSupportDofElem = globSupportDofMesh.getNumSupportDof(idGlobalElementSupport);
281
282 // fill m_iEle for this element (CSR format)
283 141947 m_iEle[cptSupportElem+1] = m_iEle[cptSupportElem] + numSupportDofElem;
284
285 // fill m_iSupportDof (add the support dof of this element) (CSR format)
286
2/2
✓ Branch 0 taken 422826 times.
✓ Branch 1 taken 141947 times.
 564773 for ( int iSupport = 0; iSupport < numSupportDofElem; iSupport++)
287 422826 m_iSupportDof[m_iEle[cptSupportElem] + iSupport] = globSupportDofMesh.iSupportDof()[globSupportDofMesh.iEle()[idGlobalElementSupport] + iSupport];
288
289 // fill m_vectorIdElementSupport
290 141947 m_vectorIdElementSupport[cptElt][iSup] = cptSupportElem;
291
292 // fill loc2GlobElemSupport
293 141947 loc2GlobElemSupport[cptSupportElem] = idGlobalElementSupport;
294
295 // increment the count of local support elements
296 141947 ++cptSupportElem;
297 }
298
299 // increment the count of local elements
300 141895 ++cptElt;
301 }
302 }
303
304 // Set m_numSupportDof
305
1/2
✓ Branch 3 taken 762 times.
✗ Branch 4 not taken.
 762 m_numSupportDof = std::set<felInt>(m_iSupportDof.begin(), m_iSupportDof.end()).size();
306 762 }
307
308 /*!
309 \brief Resize the following vectors :
310 m_iSupportDof (sum of the number of Dof in each elements),
311 m_iEle (Total number of elements + 1)
312 */
313 769 void SupportDofMesh::m_resizeSupportVectors()
314 {
315 769 felInt numDofTotal = 0;
316 769 felInt numElemTotal = 0;
317 ElementType eltType;
318 769 const std::vector<ElementType>& bagEleDomain = m_mesh->bagElementTypeDomain();
319 769 const std::vector<ElementType>& bagEleDomainBoundary = m_mesh->bagElementTypeDomainBoundary();
320
321
2/2
✓ Branch 1 taken 769 times.
✓ Branch 2 taken 769 times.
 1538 for (std::size_t ieltype=0; ieltype < bagEleDomain.size(); ++ieltype) {
322 769 eltType = bagEleDomain[ieltype];
323 769 numElemTotal += m_mesh->numElements(eltType);
324
325 //! definition of the finite element type ( geoEle, RefEle )
326 769 const GeoElement& geoEle = *GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
327
1/2
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
 769 const RefElement& refEle = m_chooseRefEleFromFEType(geoEle, eltType);
328
329 769 numDofTotal += refEle.numDOF()*m_mesh->numElements(eltType);
330 }
331
332
2/2
✓ Branch 1 taken 839 times.
✓ Branch 2 taken 769 times.
 1608 for (std::size_t ieltype=0; ieltype < bagEleDomainBoundary.size(); ++ieltype) {
333 839 eltType = bagEleDomainBoundary[ieltype];
334 839 numElemTotal += m_mesh->numElements(eltType);
335
336 839 const GeoElement& geoEle = *GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
337
1/2
✓ Branch 1 taken 839 times.
✗ Branch 2 not taken.
 839 const RefElement& refEle = m_chooseRefEleFromFEType(geoEle, eltType);
338
339 839 numDofTotal += refEle.numDOF()*m_mesh->numElements(eltType);
340 }
341
342
1/2
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
 769 m_iSupportDof.resize(numDofTotal, 0);
343
1/2
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
 769 m_iEle.resize(numElemTotal+1, 0);
344
1/2
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
 769 m_vectorIdElementSupport.resize(numElemTotal);
345 769 }
346
347 /*!
348 \brief make the link between the nodes of the (reference) finite element, and
349 the nodes of the geometric element of the mesh.
350 For instance, identify nodes between linear geometric element and quadratic treatment in the problem.
351 \param[in] fe finite element reference
352 \param[out] feIdLocalDof return the ID of the local Dof in the element, -1 if not linked.
353 for instance if GeoEle=TriaP1 and RefEle=TriaP2: feIdLocalDof=[0,1,2,-1,-1,-1].
354 */
355 3224 void SupportDofMesh::m_linkNodesGeoAndRefEle(const CurBaseFiniteElement& fe, std::vector<int>& feIdLocalDof) const
356 {
357 // The vertex dof in the definition of the refElement have to be the first ones!
358 3224 int numDOFNodeVertex = fe.refEle().numDOFNodeVertex();
359 3224 int numDOF = fe.refEle().numDOF();
360 3224 const int* idSupportOfDOF = fe.refEle().idSupportOfDOF();
361
362
1/2
✓ Branch 1 taken 3224 times.
✗ Branch 2 not taken.
 3224 feIdLocalDof.resize(numDOF, -1);
363
364
2/2
✓ Branch 0 taken 9646 times.
✓ Branch 1 taken 3224 times.
 12870 for ( int iDOFNode = 0; iDOFNode < numDOFNodeVertex; iDOFNode++)
365 9646 feIdLocalDof[iDOFNode] = idSupportOfDOF[iDOFNode];
366 3224 }
367
368 /*
369 void SupportDofMesh::m_linkNodesGeoAndRefEle(const CurvilinearFiniteElement& fe, std::vector<int>& feIdLocalDof) const
370 {
371 int numDOFNode = fe.refEle().numDOFNode();
372 int numDOF = fe.refEle().numDOF();
373
374 feIdLocalDof.clear();
375 feIdLocalDof.resize( numDOF, -1 );
376
377 const int* idSupportOfDOF = fe.refEle().idSupportOfDOF();
378 double dist = 0.;
379
380 // TODO
381 // - Change the norm (use infinite or L1 norm).
382 // NB fe is initialized here such that the refEle and the geoEle
383 // have the same vertices (reference element).
384 // -> No need to use relative norm in the comparison.
385 for ( int iDOFNode = 0; iDOFNode < numDOFNode; iDOFNode++) {
386 for ( int i = 0; i < fe.geoEle().numPoint(); i++) {
387 dist = std::sqrt( (fe.refEle().node()[iDOFNode].x() - fe.geoEle().pointCoor(i, 0))*(fe.refEle().node()[iDOFNode].x() - fe.geoEle().pointCoor(i, 0)) +
388 (fe.refEle().node()[iDOFNode].y() - fe.geoEle().pointCoor(i, 1))*(fe.refEle().node()[iDOFNode].y() - fe.geoEle().pointCoor(i, 1)) +
389 (fe.refEle().node()[iDOFNode].z() - fe.geoEle().pointCoor(i, 2))*(fe.refEle().node()[iDOFNode].z() - fe.geoEle().pointCoor(i, 2)) );
390 if ( dist < 0.0000000001 ) {
391 feIdLocalDof[iDOFNode] = idSupportOfDOF[iDOFNode];
392 }
393 }
394 }
395 }
396 */
397
398 //! \brief copy the mesh points (all points in mesh) into m_listNode
399 769 void SupportDofMesh::m_copyNodes()
400 {
401
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 769 times.
 769 if ( m_copiedMeshPoint ) {
402 ✗ return;
403 } else {
404 // copy nodes only once!
405 769 m_listNode = m_mesh->listPoints();
406 769 m_copiedMeshPoint = true;
407 }
408 }
409
410 //! \brief create edge nodes if necessary
411 168 void SupportDofMesh::m_resizeEdgeNodes()
412 {
413
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 168 times.
 168 if ( m_resizedEdgeNode ) {
414 ✗ return;
415 } else { // resize only once!
416 168 felInt size = m_listNode.size() + m_numDofSupportedByEdge;
417 168 m_listNode.resize(size); // expand the std::vector
418 168 m_resizedEdgeNode = true;
419 }
420 }
421
422 //! \brief create internal faces if necessary
423 ✗ void SupportDofMesh::m_resizeFaceNodes(const RefElement& refEle)
424 {
425 ✗ if ( m_resizedFaceNode ) {
426 ✗ return;
427 } else { // resize only once!
428 ✗ int numDOFPerFace = refEle.numDOFPerFace();
429 ✗ felInt size = m_listNode.size() + m_mesh->numFaces() * numDOFPerFace;
430 ✗ m_listNode.resize(size); // expand the std::vector
431 ✗ m_resizedFaceNode = true;
432 }
433 }
434
435 //! \brief create internal volumes if necessary
436 ✗ void SupportDofMesh::m_resizeVolumeNodes(const RefElement& refEle)
437 {
438 ✗ if ( m_resizedVolumeNode ) {
439 ✗ return;
440 } else { // resize only once!
441 ✗ int numDOFPerVolume = refEle.numDOFVolume();
442 ✗ felInt size = m_listNode.size() + m_mesh->getNumElement3D() * numDOFPerVolume;
443 ✗ m_listNode.resize(size); // expand the std::vector
444 ✗ m_resizedVolumeNode = true;
445 }
446 }
447
448 //! build and add the edge and face NODES to m_listNode
449 //! (typically build the edge midpoints for constructing P2/Q2 from P1/Q1 mesh).
450 1538 void SupportDofMesh::m_buildNodesOfEdgeFaceVolumePerBag(const std::vector<ElementType>& theBagElt, felInt* countEltPerType, bool boundary_flag)
451 {
452 1538 std::vector<felInt> elem;
453 1538 std::vector<Point*> elePoint;
454 1538 std::vector<felInt> elemConnectivity;
455 1538 std::vector<int> feIdLocalDof;
456
2/2
✓ Branch 1 taken 1608 times.
✓ Branch 2 taken 1538 times.
 3146 for (std::size_t ieltype = 0; ieltype < theBagElt.size(); ++ieltype) {
457 1608 ElementType eltType = theBagElt[ieltype];
458 1608 felInt numPtsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
459
460 //! definition of the finite element type ( geoEle, RefEle )
461 1608 const GeoElement& geoEle = *GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
462
1/2
✓ Branch 1 taken 1608 times.
✗ Branch 2 not taken.
 1608 const RefElement& refEle = m_chooseRefEleFromFEType(geoEle, eltType);
463
464 1608 CurBaseFiniteElement::Pointer fe = nullptr;
465
466
2/2
✓ Branch 0 taken 839 times.
✓ Branch 1 taken 769 times.
 1608 if ( boundary_flag ) { // curvFE necessary for the boundary elements to have
467 // A correct coorMap for ALL coordinates (because numCoor(Bndry) = numCoor(Domain)-1)
468
1/2
✓ Branch 1 taken 839 times.
✗ Branch 2 not taken.
 839 fe = felisce::make_shared<CurvilinearFiniteElement>(refEle, geoEle, DegreeOfExactness_0);
469 } else {
470
1/2
✓ Branch 1 taken 769 times.
✗ Branch 2 not taken.
 769 fe = felisce::make_shared<CurrentFiniteElement>(refEle, geoEle, DegreeOfExactness_0);
471 }
472
473 1608 const int numDOFNodeVertex = refEle.numDOFNodeVertex();
474 1608 const int numDOFNodeEdge = refEle.numDOFNodeEdge(); // numDOFNodeEdge or numDOFEdge here ??? VM
475 1608 const int numDOFNodeFace = refEle.numDOFNodeFace();
476 1608 const int numDOFNodeVolume = refEle.numDOFNodeVolume();
477
478 // check that m_listNode is resized before inserting the nodes:
479
3/4
✓ Branch 0 taken 1172 times.
✓ Branch 1 taken 436 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1172 times.
 1608 FEL_ASSERT( m_resizedEdgeNode || numDOFNodeEdge == 0 ); //check: numDOFNodeEdge > 0 => m_resizedEdgeNode = true.
480
2/4
✓ Branch 0 taken 1608 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1608 times.
 1608 FEL_ASSERT( m_resizedFaceNode || numDOFNodeFace == 0 );
481
4/6
✓ Branch 0 taken 324 times.
✓ Branch 1 taken 1284 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 324 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1608 times.
 1608 if ( numDOFNodeEdge > 0 && m_mesh->statusEdges() == false ) {
482 ✗ FEL_ERROR("Error: you must build the Edges BEFORE the build of Edge nodes!");
483 }
484
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1608 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1608 times.
 1608 if ( numDOFNodeFace > 0 && m_mesh->statusFaces() == false ) {
485 ✗ FEL_ERROR("Error: you must build the Faces BEFORE the build of Face nodes!");
486 }
487
488 // Link geoEle and refElement nodes (P1->P2 for instance)
489
1/2
✓ Branch 2 taken 1608 times.
✗ Branch 3 not taken.
 1608 m_linkNodesGeoAndRefEle(*fe, feIdLocalDof);
490
491
1/2
✓ Branch 1 taken 1608 times.
✗ Branch 2 not taken.
 1608 elem.resize(GeometricMeshRegion::m_numPointsPerElt[eltType], 0);
492
1/2
✓ Branch 1 taken 1608 times.
✗ Branch 2 not taken.
 1608 elePoint.resize(GeometricMeshRegion::m_numPointsPerElt[eltType], nullptr);
493
1/2
✓ Branch 2 taken 1608 times.
✗ Branch 3 not taken.
 1608 elemConnectivity.resize(refEle.numDOF(),0);
494
495
2/2
✓ Branch 6 taken 3753 times.
✓ Branch 7 taken 1608 times.
 5361 for(auto itRef = m_mesh->intRefToBegEndMaps[eltType].begin();itRef != m_mesh->intRefToBegEndMaps[eltType].end(); itRef++) {
496 3753 const felInt numElemsPerRef = itRef->second.second;
497
498 // Build and add the edge middle points to the mesh.
499
500 //??? Question: pourquoi ne pas travailler par liste de faces/d'aretes au lieu de passer par les elements??? VM
501
2/2
✓ Branch 0 taken 5113818 times.
✓ Branch 1 taken 3753 times.
 5117571 for ( felInt iel = 0; iel < numElemsPerRef; iel++ ) {
502 //! get the elem and its connectivity
503
1/2
✓ Branch 1 taken 5113818 times.
✗ Branch 2 not taken.
 5113818 m_getElementSupportConnectivity(refEle, eltType, countEltPerType[eltType], feIdLocalDof, elem, elemConnectivity);
504
2/2
✓ Branch 0 taken 19331382 times.
✓ Branch 1 taken 5113818 times.
 24445200 for ( int ivert = 0; ivert < numPtsPerElt; ivert++) {
505 19331382 elePoint[ivert] = &m_mesh->listPoint( elem[ivert] );
506 }
507
1/2
✓ Branch 2 taken 5113818 times.
✗ Branch 3 not taken.
 5113818 fe->update(iel, elePoint); // geometric update only
508
509 // TODO : CHECK THAT IT WORKS FOR MULTIPLE DOF PER EDGE/FACE! VM
510
511 //! computing edge nodes
512
2/2
✓ Branch 0 taken 130996 times.
✓ Branch 1 taken 4982822 times.
 5113818 if ( numDOFNodeEdge > 0 ) {
513 Point middleEdge;
514
2/2
✓ Branch 0 taken 513544 times.
✓ Branch 1 taken 130996 times.
 644540 for ( int iDof = 0; iDof < numDOFNodeEdge; iDof++) { // numDOFNodeEdge or numDOFEdge here ??? VM
515
1/2
✓ Branch 1 taken 513544 times.
✗ Branch 2 not taken.
 513544 if ( feIdLocalDof[iDof + numDOFNodeVertex] == -1 ) { //add the midpoint if not yet created
516
517
1/2
✓ Branch 3 taken 513544 times.
✗ Branch 4 not taken.
 513544 fe->coorMap(middleEdge, refEle.node()[numDOFNodeVertex + iDof]);
518 513544 m_listNode[ elemConnectivity[numDOFNodeVertex + iDof] ] = middleEdge;
519 }
520 }
521 }
522
523 //! computing face nodes
524
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5113818 times.
 5113818 if ( numDOFNodeFace > 0 ) {
525 Point middleFace;
526 ✗ for ( int iDof = 0; iDof < numDOFNodeFace; iDof++) {
527 ✗ if ( feIdLocalDof[iDof + numDOFNodeVertex + numDOFNodeEdge] == -1 ) {
528 ✗ fe->coorMap(middleFace, refEle.node()[numDOFNodeVertex + numDOFNodeEdge + iDof]);
529 ✗ m_listNode[ elemConnectivity[numDOFNodeVertex + numDOFNodeEdge + iDof] ] = middleFace;
530 }
531 }
532 }
533
534 //! computing volume nodes
535
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5113818 times.
 5113818 if ( numDOFNodeVolume > 0 ) {
536 Point middleVolume;
537 ✗ int shift = numDOFNodeVertex + numDOFNodeEdge + numDOFNodeFace;
538 ✗ for ( int iDof = 0; iDof < numDOFNodeVolume; iDof++) {
539 ✗ if ( feIdLocalDof[shift + iDof] == -1 ) {
540 ✗ fe->coorMap(middleVolume, refEle.node()[shift + iDof]);
541 ✗ m_listNode[ elemConnectivity[shift + iDof] ] = middleVolume;
542 }
543 }
544 }
545 5113818 countEltPerType[eltType]++;
546 }
547 }
548 1608 }
549 1538 }
550
551 //! return:
552 //! elem : the element.
553 //! elemConnectivity: the ID of the nodes of the element. Ex. P2: [V0, V1, V2, Mid0, Mid1, Mid2]
554 10232244 void SupportDofMesh::m_getElementSupportConnectivity(const RefElement& refEle, ElementType eltType, felInt iel, const std::vector<int>& feIdLocalDof,
555 std::vector<felInt>& the_elem, std::vector<felInt>& elemConnectivity) const
556 {
557 //const int numDOFNode = refEle.numDOFNode();
558 10232244 const int numDOFNodeVertex = refEle.numDOFNodeVertex();
559 10232244 const int numDOFNodeEdge = refEle.numDOFNodeEdge();
560 10232244 const int numDOFNodeFace = refEle.numDOFNodeFace();
561 10232244 const int numDOFNodeVolume = refEle.numDOFNodeVolume();
562 10232244 const int numDOFEdge = refEle.numDOFEdge();
563 10232244 const int numDOFFace = refEle.numDOFFace();
564 10232244 const int numDOFVolume = refEle.numDOFVolume();
565 10232244 const int numDOF = refEle.numDOF();
566 10232244 const int numDOFPerFace = refEle.numDOFPerFace();
567
568 10232244 const GeoElement& geoEle = *GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
569 10232244 const int numPoint = geoEle.numPoint();
570
571
1/2
✗ Branch 3 not taken.
✓ Branch 4 taken 10232244 times.
 10232244 FEL_ASSERT( the_elem.size() == static_cast<std::size_t>(m_mesh->numPtsPerElement(eltType)) );
572
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 10232244 times.
 10232244 FEL_ASSERT( elemConnectivity.size() == static_cast<std::size_t>(numDOF) );
573
574 //! We assume that either the edge/face dofs are fully nodal, or fully non-nodal.
575
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10232244 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 10232244 FEL_ASSERT(numDOFNodeEdge == numDOFEdge || numDOFNodeEdge == 0);
576
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10232244 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 10232244 FEL_ASSERT(numDOFNodeFace == numDOFFace || numDOFNodeFace == 0);
577
578
1/2
✓ Branch 2 taken 10232244 times.
✗ Branch 3 not taken.
 10232244 m_mesh->getOneElement(eltType, iel, the_elem, 0);
579
580 //! IDs of the edges (or faces) of the element (as provided by listEdges, given edges first).
581
1/2
✓ Branch 2 taken 10232244 times.
✗ Branch 3 not taken.
 10232244 std::vector<felInt> elem2Edge(GeometricMeshRegion::m_numEdgesPerElt[eltType], -1);
582
1/2
✓ Branch 2 taken 10232244 times.
✗ Branch 3 not taken.
 10232244 std::vector<felInt> elem2Face(GeometricMeshRegion::m_numFacesPerElt[eltType], -1);
583
584 10232244 int shiftLocalDOF = 0;
585 10232244 felInt shiftGlobalDOF = 0;
586
587 //! 1/ VERTEX nodes
588
1/2
✓ Branch 0 taken 10232244 times.
✗ Branch 1 not taken.
 10232244 if ( numDOFNodeVertex != 0 ) {
589 // Is this correct for non nodal dof, such as RT0 or P0: ??? VM
590 // start with vertex nodes.
591
2/2
✓ Branch 1 taken 38676204 times.
✓ Branch 2 taken 10232244 times.
 48908448 for (std::size_t iv = 0; iv < the_elem.size(); iv++)
592 38676204 elemConnectivity[iv] = the_elem[iv];
593
594 10232244 shiftLocalDOF += numDOFNodeVertex;
595 10232244 shiftGlobalDOF += m_mesh->numPoints();
596 }
597
598 //! 2/ EDGE dofs
599 //! It assumes that for all elements, there is the same number of dof per edges.
600
4/4
✓ Branch 0 taken 261992 times.
✓ Branch 1 taken 9970252 times.
✓ Branch 2 taken 232296 times.
✓ Branch 3 taken 29696 times.
 10232244 if ( numDOFEdge != 0 && numDOF != numPoint) {
601 // look only at the 1rst DOFNodeEdge
602 232296 int idof = 0;
603
604 //! if there is a node dof supported by an edge (ex: refEle=P2, built from geoEle=P1)
605 //! or a non nodal dof on an edge
606
3/8
✓ Branch 0 taken 232296 times.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 232296 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 232296 times.
✗ Branch 8 not taken.
 232296 if ( (numDOFNodeEdge > 0 && feIdLocalDof[idof + shiftLocalDOF] == -1) || (numDOFNodeEdge == 0) ) {
607 // get all the edges of the element
608
1/2
✓ Branch 2 taken 232296 times.
✗ Branch 3 not taken.
 232296 m_mesh->getAllEdgeOfElement(eltType, iel, elem2Edge);
609 // elem2Edge contains the ids of the edges,
610 // i.e. m_listEdge[elem2Edge[i]] returns the pointer to the i-th edge of the element
611
612 // add the edge dofs after the vertex node dofs
613 232296 int localCount(0);
614
2/2
✓ Branch 1 taken 959808 times.
✓ Branch 2 taken 232296 times.
 1192104 for ( std::size_t iedg = 0; iedg < elem2Edge.size(); iedg++ ) {
615
2/2
✓ Branch 1 taken 941120 times.
✓ Branch 2 taken 959808 times.
 1900928 for ( int jdof = 0; jdof < refEle.numDOFPerEdge(iedg); jdof++ ) {
616
2/2
✓ Branch 1 taken 9856 times.
✓ Branch 2 taken 931264 times.
 941120 if (refEle.edgesCanHaveDifferentNumberOfDof()) {
617
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9856 times.
 9856 FEL_ASSERT( refEle.numDOFPerEdge(iedg) == 1 );
618 9856 elemConnectivity[shiftLocalDOF + localCount] = shiftGlobalDOF + m_mesh->listEdges()[elem2Edge[iedg]]->idOnlySupporting() + jdof;
619 9856 localCount++;
620 } else {
621
1/2
✓ Branch 2 taken 931264 times.
✗ Branch 3 not taken.
 931264 elemConnectivity[shiftLocalDOF + localCount] = jdof + elem2Edge[iedg]*refEle.numDOFPerEdge() + shiftGlobalDOF;
622 931264 localCount++;
623 }
624 }
625 }
626 }
627 232296 shiftLocalDOF += numDOFEdge;
628 232296 shiftGlobalDOF += m_mesh->listEdges().numEdgesSupportingADof();
629 }
630
631 //! 2/ FACE dofs
632 //! It assumes that for all elements, there is the same number of dof per faces
633
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10232244 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 10232244 if ( numDOFFace != 0 && numDOF != numPoint) {
634 // look only at the 1rst DOFNodeFace
635 ✗ int idof = 0;
636
637 //! If there is a node dof supported by a face or a non nodal dof on a face
638 ✗ if ( (numDOFNodeFace > 0 && feIdLocalDof[idof + shiftLocalDOF] == -1) || (numDOFNodeFace == 0) ) {
639 // get all the faces of the element
640 ✗ m_mesh->getAllFaceOfElement(eltType,iel,elem2Face);
641
642 // add the face dofs after the edge node dofs
643 ✗ for ( std::size_t ifac = 0; ifac < elem2Face.size(); ifac++ ) {
644 ✗ for ( int jdof = 0; jdof < numDOFPerFace; jdof++ ) {
645 ✗ elemConnectivity[shiftLocalDOF + ifac*numDOFPerFace + jdof] = jdof + elem2Face[ifac]*numDOFPerFace + shiftGlobalDOF;
646 }
647 }
648 }
649 ✗ shiftLocalDOF += numDOFFace; //or numDOFNodeFace; ???
650 ✗ shiftGlobalDOF += m_mesh->numFaces() * numDOFPerFace;
651 }
652
653 //! 3/ VOLUME dofs
654 //! It assumes that there is always the same number of volume dof for all elements
655
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10232244 times.
 10232244 if ( numDOFVolume != 0 ) {
656 // look only at the 1rst DOFNodeVolume
657 ✗ int idof = 0;
658
659 ✗ int elem2Volume = 0;
660 ✗ if ( numDOFNodeVolume > 0 && feIdLocalDof[idof + shiftLocalDOF] == -1 ) {
661 // compute elem2Volume
662 ✗ auto& bagElementType3D = GeometricMeshRegion::bagElementType3D;
663 ✗ for (auto it_elttype = bagElementType3D.begin(); it_elttype != bagElementType3D.end(); ++it_elttype) {
664 ✗ ElementType elementType = (ElementType)*it_elttype;
665 ✗ if(elementType == eltType) {
666 ✗ elem2Volume += iel;
667 } else {
668 ✗ elem2Volume += m_mesh->numElements(elementType);
669 }
670 }
671
672 // add the volume dofs after the face dofs.
673 ✗ for ( int jdof = 0; jdof < numDOFVolume; jdof++ )
674 ✗ elemConnectivity[shiftLocalDOF + jdof] = shiftGlobalDOF + elem2Volume*numDOFVolume + jdof;
675 }
676 }
677 10232244 }
678
679 //! helper to build the support arrays in CSR format
680 //!m_iEle, m_iSupportDof :
681 1538 void SupportDofMesh::m_buildSupportCSR(const std::vector<ElementType>& theBagElt, felInt& cptElt, felInt* countEltPerType)
682 {
683 1538 std::vector<felInt> elem;
684 1538 std::vector<felInt> elemConnectivity;
685 1538 std::vector<int> feIdLocalDof;
686
687
2/2
✓ Branch 1 taken 1608 times.
✓ Branch 2 taken 1538 times.
 3146 for (std::size_t ieltype = 0; ieltype < theBagElt.size(); ++ieltype) {
688 1608 ElementType eltType = theBagElt[ieltype];
689
690 //! definition of the finite element type ( geoEle, RefEle )
691 1608 const GeoElement& geoEle = *GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
692
1/2
✓ Branch 1 taken 1608 times.
✗ Branch 2 not taken.
 1608 const RefElement& refEle = m_chooseRefEleFromFEType(geoEle, eltType);
693
694
1/2
✓ Branch 1 taken 1608 times.
✗ Branch 2 not taken.
 1608 CurrentFiniteElement fe(refEle, geoEle, DegreeOfExactness_0);
695 // Link geo and ref nodes
696
1/2
✓ Branch 1 taken 1608 times.
✗ Branch 2 not taken.
 1608 m_linkNodesGeoAndRefEle(fe, feIdLocalDof);
697
698
1/2
✓ Branch 2 taken 1608 times.
✗ Branch 3 not taken.
 1608 elemConnectivity.resize(refEle.numDOF(),0);
699
1/2
✓ Branch 1 taken 1608 times.
✗ Branch 2 not taken.
 1608 elem.resize(GeometricMeshRegion::m_numPointsPerElt[eltType], 0);
700
701
2/2
✓ Branch 6 taken 3753 times.
✓ Branch 7 taken 1608 times.
 5361 for(auto itRef = m_mesh->intRefToBegEndMaps[eltType].begin(); itRef != m_mesh->intRefToBegEndMaps[eltType].end(); itRef++) {
702 // felInt currentLabel = itRef->first;
703 3753 felInt numElemsPerRef = itRef->second.second;
704
705
2/2
✓ Branch 0 taken 5113818 times.
✓ Branch 1 taken 3753 times.
 5117571 for ( felInt iel = 0; iel < numElemsPerRef; iel++ ) {
706 //! get the elem and its connectivity
707
1/2
✓ Branch 1 taken 5113818 times.
✗ Branch 2 not taken.
 5113818 m_getElementSupportConnectivity(refEle, eltType, countEltPerType[eltType], feIdLocalDof, elem, elemConnectivity);
708
709
3/4
✓ Branch 1 taken 5113818 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 70636 times.
✓ Branch 4 taken 5043182 times.
 5113818 if(FelisceParam::instance().FusionDof) {
710
2/2
✓ Branch 1 taken 262896 times.
✓ Branch 2 taken 70636 times.
 333532 for (std::size_t iDof = 0; iDof < elemConnectivity.size(); iDof++) {
711 262896 m_iSupportDof[m_iEle[cptElt] + iDof ] = m_listPerm[elemConnectivity[iDof]];
712 }
713
714 70636 m_iEle[cptElt + 1] = m_iEle[cptElt] + elemConnectivity.size();
715
1/2
✓ Branch 2 taken 70636 times.
✗ Branch 3 not taken.
 70636 m_vectorIdElementSupport[cptElt].push_back(cptElt);
716 70636 cptElt ++;
717 70636 countEltPerType[eltType]++;
718 } else {
719 // Link from geometric element to supportdofmesh: CSR format
720
2/2
✓ Branch 1 taken 19539046 times.
✓ Branch 2 taken 5043182 times.
 24582228 for (std::size_t iDof = 0; iDof < elemConnectivity.size(); iDof++) {
721 19539046 m_iSupportDof[m_iEle[cptElt] + iDof ] = elemConnectivity[iDof];
722 }
723
724 5043182 m_iEle[cptElt + 1] = m_iEle[cptElt] + elemConnectivity.size();
725
1/2
✓ Branch 2 taken 5043182 times.
✗ Branch 3 not taken.
 5043182 m_vectorIdElementSupport[cptElt].push_back(cptElt);
726 5043182 cptElt++;
727 5043182 countEltPerType[eltType]++;
728 }
729 }
730 }
731 1608 }
732 1538 }
733
734 //! to switch from geoEle P1 to refEle P2 for instance
735 5601 const RefElement& SupportDofMesh::m_chooseRefEleFromFEType(const GeoElement& geoEle, const ElementType eltType) const
736 {
737 // Choice of the reference finite element (0: default, 1: quadratic)
738
3/4
✓ Branch 1 taken 5089 times.
✓ Branch 2 taken 432 times.
✓ Branch 3 taken 80 times.
✗ Branch 4 not taken.
 5601 switch ( m_variable.finiteElementType() ) {
739 5089 case 0: { //! default
740 5089 return geoEle.defaultFiniteEle();
741 }
742 432 case 1: { //! quadratic
743 432 const GeoElement& geoEleBis = *GeometricMeshRegion::eltEnumToFelNameGeoEle[GeometricMeshRegion::eltLinearToEltQuad[eltType]].second;
744 432 return geoEleBis.defaultFiniteEle();
745 }
746 80 case 2: { //! bubble
747
2/2
✓ Branch 3 taken 32 times.
✓ Branch 4 taken 48 times.
 80 if(m_mesh->domainDim() == geoEle.numCoor()) {
748 32 const GeoElement& geoEleBis = *GeometricMeshRegion::eltEnumToFelNameGeoEle[GeometricMeshRegion::eltLinearToEltBubble[eltType]].second;
749 32 return geoEleBis.defaultFiniteEle();
750 } else {
751 // This looks ugly, I know, but the alternative would be to add another case after bubble finite element and add a lot of useless mappings.
752 // If someone ever tries to use P1bubble in prisms this has to be changed!!
753
6/8
✓ Branch 3 taken 48 times.
✗ Branch 4 not taken.
✓ Branch 8 taken 48 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 24 times.
✓ Branch 11 taken 24 times.
✓ Branch 12 taken 24 times.
✓ Branch 13 taken 24 times.
 48 if (m_mesh->bagElementTypeDomain().size()==1 && m_mesh->bagElementTypeDomain()[0] == GeometricMeshRegion::Prism6 && eltType == GeometricMeshRegion::Quad4)
754 24 return refElementQuadrangleP1xP2;
755 24 return geoEle.defaultFiniteEle();
756 }
757 }
758 ✗ default:
759 ✗ FEL_ERROR("Error: the finite Element type (linear, quadratic, ...) is not set properly.");
760 }
761 return refElementNode;
762 }
763
764 /*!
765 \brief Get the id of the element from its type and id with respect to its type
766 \param[in] mesh The mesh where to get the id from
767 \param[in] eltType The type of element
768 \param[in] ielGeom The id of the element
769 \param[out] ielSupportDof The id of the element
770 */
771 6761930 void SupportDofMesh::getIdElementSupport(const ElementType& eltType, felInt ielGeom, std::vector<felInt>& vectorSupport) const
772 {
773 6761930 felInt position = 0;
774
1/2
✓ Branch 2 taken 6761930 times.
✗ Branch 3 not taken.
 6761930 m_mesh->getIdElemFromTypeElemAndIdByType(eltType, ielGeom, position);
775
1/2
✓ Branch 2 taken 6761930 times.
✗ Branch 3 not taken.
 6761930 vectorSupport = m_vectorIdElementSupport[position];
776 6761930 }
777
778 /*!
779 \brief Get the id of the element from its type and id with respect to its type
780 \param[in] mesh The mesh where to get the id from
781 \param[in] eltType The type of element
782 \param[in] ielGeom The id of the element
783 \param[out] ielSupportDof The id of the element
784 */
785 27175653 void SupportDofMesh::getIdElementSupport(const ElementType& eltType, felInt ielGeom, felInt& ielSupportDof) const
786 {
787 27175653 felInt position = 0;
788
1/2
✓ Branch 2 taken 27175653 times.
✗ Branch 3 not taken.
 27175653 m_mesh->getIdElemFromTypeElemAndIdByType(eltType, ielGeom, position);
789 27175653 ielSupportDof = m_vectorIdElementSupport[position][0];
790 27175653 }
791
792 10378710 void SupportDofMesh::getIdElementSupport(felInt idEle, std::vector<felInt>& vecSupport) const
793 {
794 10378710 vecSupport = m_vectorIdElementSupport[idEle];
795 10378710 }
796
797 2500089 void SupportDofMesh::getIdElementSupport(felInt idEle, felInt& ielSupportDof) const
798 {
799 2500089 ielSupportDof = m_vectorIdElementSupport[idEle][0];
800 2500089 }
801
802 ✗ void SupportDofMesh::getIdPointElementSupport(felInt idEle, std::vector<felInt>& vecIdPoint) const
803 {
804 ✗ for ( int i = 0; i < getNumSupportDof(idEle); ++i)
805 ✗ vecIdPoint[i] = m_iSupportDof[m_iEle[idEle]+i];
806 }
807
808 20 void SupportDofMesh::fusionSupportDof(const Variable& variable)
809 {
810 20 std::size_t cpt = 0;
811 std::size_t label, labelIN, labelOUT;
812 std::size_t FusionNumLabPerVariable;
813 20 std::vector< std::pair<std::size_t, std::size_t> > listSuppDofMatched;
814 20 std::set<std::size_t> listSuppDofIN;
815 20 std::set<std::size_t> listSuppDofOUT;
816 20 std::set<std::size_t> listSuppDof;
817 20 std::set<std::size_t> listSuppDofToMatched;
818 double distNodes;
819 20 std::vector<int> elem, elemConnectivity, feIdLocalDof;
820 std::size_t countEltPerType[ GeometricMeshRegion::m_numTypesOfElement ];
821 20 std::size_t cptElt = 0;
822
823
2/2
✓ Branch 0 taken 500 times.
✓ Branch 1 taken 20 times.
 520 for (int itType=0; itType<GeometricMeshRegion::m_numTypesOfElement; itType++ )
824 500 countEltPerType[itType] = 0;
825
826
3/4
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✓ Branch 5 taken 20 times.
 40 for(std::size_t fvar = 0; fvar < FelisceParam::instance().FusionVariable.size(); fvar++) { // loop on fusion variables
827
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 FusionNumLabPerVariable = FelisceParam::instance().FusionNumLabel[fvar];
828
4/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 8 times.
✓ Branch 11 taken 12 times.
 20 if ( FelisceParam::instance().FusionVariable[fvar] == variable.name().c_str()) {
829
830
3/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 8 times.
 8 if ( FelisceParam::instance().AllToOne[fvar] ) { // all nodes fusioned into node 0
831 ✗ for ( std::size_t i= 1; i < m_listNode.size(); ++i )
832 ✗ listSuppDofMatched.emplace_back(0,i);
833 }
834
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 8 times.
 8 else if ( FelisceParam::instance().AllInLabelToOne ) { // all nodes of a label fusioned into a single node
835 // get the nodes
836 ✗ PetscPrintf(PETSC_COMM_WORLD,"##### For variable %s \n",variable.name().c_str() );
837 ✗ for(std::size_t flab = 0; flab < FusionNumLabPerVariable; flab++) { // loop on labels
838 ✗ label = FelisceParam::instance().FusionLabel[cpt + flab];
839 ✗ PetscPrintf(PETSC_COMM_WORLD,"Fusion dof for label %ld \n", static_cast<unsigned long>(label) );
840 ✗ getElemConnectivitySingleLabel(label, listSuppDof); // get list of all nodes with label = labelDarcy
841 ✗ for (auto it = listSuppDof.begin(); it!=listSuppDof.end(); ++it) {
842 ✗ if ( it != listSuppDof.begin() )
843 ✗ listSuppDofMatched.emplace_back(*listSuppDof.begin(),*it);
844 }
845 }
846 }
847
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 8 times.
 8 else if ( FelisceParam::instance().contactDarcy > 0 ) {
848 // =================================================
849 // MF : This part has to be reworked, it is to specifiy to the Darcy test case
850 //==================================================
851 // WARNING: done for only one type of label per variable
852 // get the list of darcy global nodes
853 ✗ PetscPrintf(PETSC_COMM_WORLD,"##### For variable %s \n",variable.name().c_str() );
854 ✗ for(std::size_t flab = 0; flab < FusionNumLabPerVariable; flab++) { // loop on labels
855 ✗ label = FelisceParam::instance().FusionLabel[cpt + flab];
856 ✗ PetscPrintf(PETSC_COMM_WORLD,"Fusion dof for label %ld \n", static_cast<unsigned long>(label) );
857 ✗ getElemConnectivitySingleLabel(label, listSuppDof); // get list of all nodes with label = labelDarcy
858 }
859
860 // std::set the nodeId for the fusion in m_idFusionDarcy
861 ✗ m_idFusionDarcy = *listSuppDof.begin();
862
863 // loop on the mesh nodes and match the nodes not on the list, with the first darcy node
864 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_mesh->bagElementTypeDomain();
865 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
866 ✗ ElementType eltType = bagElementTypeDomain[i];
867 ✗ const GeoElement& geoEle = *GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
868 ✗ const RefElement& refEle = m_chooseRefEleFromFEType(geoEle, eltType);
869 ✗ CurrentFiniteElement fe(refEle, geoEle, DegreeOfExactness_0);
870 ✗ m_linkNodesGeoAndRefEle(fe, feIdLocalDof);
871 ✗ elemConnectivity.resize(refEle.numDOF(),0);
872 ✗ elem.resize(GeometricMeshRegion::m_numPointsPerElt[eltType], 0);
873
874 // second loop on label of the mesh.
875 ✗ for(auto itRef = m_mesh->intRefToBegEndMaps[eltType].begin();
876 ✗ itRef != m_mesh->intRefToBegEndMaps[eltType].end(); itRef++) {
877 ✗ const std::size_t numEltPerLabel = itRef->second.second;
878 ✗ for (std::size_t iel = 0; iel < numEltPerLabel; iel++ ) {
879 ✗ m_getElementSupportConnectivity(refEle, eltType, countEltPerType[eltType], feIdLocalDof, elem, elemConnectivity);
880
881 ✗ for (std::size_t iDof = 0; iDof < elemConnectivity.size(); iDof++) {
882 ✗ auto kelin = std::find(listSuppDof.begin(), listSuppDof.end(),elemConnectivity[iDof]);
883 ✗ if ( kelin == listSuppDof.end() ) // if the node is not inside darcy_nodes_list we add it to the std::set of all nodes to be matched
884 ✗ listSuppDofToMatched.insert(elemConnectivity[iDof]); // we first save inside a std::set all the nodes to be matched after
885 //listSuppDofMatched.push_back(std::make_pair(*listSuppDof.begin(),elemConnectivity[iDof])); // the pair is not unique!
886 }
887 ✗ cptElt++;
888 ✗ countEltPerType[eltType]++;
889 }
890 }
891 }
892 // creation of std::vector<pair> listSuppDofMatched from the std::set listSuppDofToMatched
893 ✗ for (auto itMatch=listSuppDofToMatched.begin(); itMatch!=listSuppDofToMatched.end(); ++itMatch)
894 ✗ listSuppDofMatched.emplace_back(*listSuppDof.begin(),*itMatch);
895 }
896 else {
897
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
 8 PetscPrintf(PETSC_COMM_WORLD,"##### For variable %s \n",variable.name().c_str() );
898
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 16 for(std::size_t flab = 0; flab < FusionNumLabPerVariable; flab++) { // loop on labels
899
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 labelIN = FelisceParam::instance().FusionLabel[cpt + 2*flab];
900
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 labelOUT = FelisceParam::instance().FusionLabel[cpt + 2*flab +1];
901
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 PetscPrintf(PETSC_COMM_WORLD,"Fusion label %ld in label %ld \n", static_cast<unsigned long>(labelIN), static_cast<unsigned long>(labelOUT) );
902
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 matchingElemConnectivity(labelIN, labelOUT, listSuppDofIN, listSuppDofOUT);
903 }
904
905 //loop to delete the supporDofIN in common with the supporDofOUT
906
2/2
✓ Branch 4 taken 268 times.
✓ Branch 5 taken 8 times.
 276 for (auto kelout=listSuppDofOUT.begin(); kelout != listSuppDofOUT.end() ; ++kelout ) {
907
1/2
✓ Branch 4 taken 268 times.
✗ Branch 5 not taken.
 268 auto kelin = std::find(listSuppDofIN.begin(), listSuppDofIN.end(), *kelout);
908
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 268 times.
 268 if ( kelin != listSuppDofIN.end() )
909 ✗ listSuppDofIN.erase( kelin );
910 }
911
912 //loop to match the SuppDofIN with the SuppDofOUT
913
2/2
✓ Branch 4 taken 268 times.
✓ Branch 5 taken 8 times.
 276 for (auto kelout=listSuppDofOUT.begin(); kelout != listSuppDofOUT.end() ; ++kelout ) {
914 268 Point& pout = listNode()[*kelout];
915
2/2
✓ Branch 4 taken 13780 times.
✓ Branch 5 taken 268 times.
 14048 for (auto kelin=listSuppDofIN.begin(); kelin != listSuppDofIN.end() ; ++kelin ) {
916 13780 Point& pin = listNode()[*kelin];
917 13780 distNodes = std::sqrt( ( pin.x() - pout.x() ) * ( pin.x() - pout.x() ) + ( pin.y() - pout.y() ) * ( pin.y() - pout.y() ) + ( pin.z() - pout.z() ) * ( pin.z() - pout.z() ) );
918
3/4
✓ Branch 1 taken 13780 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 268 times.
✓ Branch 4 taken 13512 times.
 13780 if ( distNodes < FelisceParam::instance().FusionTolerance )
919
1/2
✓ Branch 3 taken 268 times.
✗ Branch 4 not taken.
 268 listSuppDofMatched.emplace_back(*kelin,*kelout);
920 }
921 }
922 }
923 8 m_numSupportDof = m_listNode.size() - listSuppDofMatched.size();
924
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 GetPermutationList(listSuppDofMatched);
925 }
926
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
 20 if ( FelisceParam::instance().AllInLabelToOne )
927 ✗ cpt++;
928 else
929 20 cpt += FusionNumLabPerVariable*2;
930 }
931 20 }
932
933 ✗ void SupportDofMesh::getElemConnectivitySingleLabel(std::size_t label, std::set<std::size_t>& listSuppDof)
934 {
935 ✗ std::vector<int> elem, elemConnectivity, feIdLocalDof;
936 ✗ std::size_t cptElt = 0;
937 std::size_t countEltPerType[ GeometricMeshRegion::m_numTypesOfElement ];
938
939 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ )
940 ✗ countEltPerType[ityp] = 0;
941
942 ✗ for (std::size_t ieltype = 0; ieltype < m_mesh->bagElementTypeDomainBoundary().size(); ++ieltype) {
943 ✗ ElementType eltType = m_mesh->bagElementTypeDomainBoundary()[ieltype];
944 ✗ const GeoElement& geoEle = *GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
945 ✗ const RefElement& refEle = m_chooseRefEleFromFEType(geoEle, eltType);
946 ✗ CurrentFiniteElement fe(refEle, geoEle, DegreeOfExactness_0);
947 ✗ m_linkNodesGeoAndRefEle(fe, feIdLocalDof);
948 ✗ elemConnectivity.resize(refEle.numDOF(),0); // contains the supportDofs points
949 ✗ elem.resize(GeometricMeshRegion::m_numPointsPerElt[eltType], 0); // contains the nodes' points (equal for P1)
950
951 ✗ for(auto itRef = m_mesh->intRefToBegEndMaps[eltType].begin();
952 ✗ itRef != m_mesh->intRefToBegEndMaps[eltType].end(); itRef++) {
953 ✗ std::size_t currentLabel = itRef->first;
954 ✗ std::size_t numElemsPerRef = itRef->second.second;
955 ✗ for (std::size_t iel = 0; iel < numElemsPerRef; iel++ ) {
956 ✗ m_getElementSupportConnectivity(refEle, eltType, countEltPerType[eltType], feIdLocalDof, elem, elemConnectivity);
957 ✗ for (std::size_t iDof = 0; iDof < elemConnectivity.size(); iDof++) {
958 ✗ if(currentLabel == label )
959 ✗ listSuppDof.insert(elemConnectivity[iDof]);
960 }
961 ✗ cptElt++;
962 ✗ countEltPerType[eltType]++;
963 }
964 }
965 }
966 }
967
968 8 void SupportDofMesh::matchingElemConnectivity(std::size_t labelIN, std::size_t labelOUT, std::set<std::size_t>& listSuppDofIN, std::set<std::size_t>& listSuppDofOUT)
969 {
970 8 std::vector<int> elem, elemConnectivity, feIdLocalDof;
971 8 std::size_t cptElt = 0;
972 std::size_t countEltPerType[ GeometricMeshRegion::m_numTypesOfElement ];
973
974
2/2
✓ Branch 0 taken 200 times.
✓ Branch 1 taken 8 times.
 208 for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ )
975 200 countEltPerType[ityp] = 0;
976
977
2/2
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 8 times.
 16 for (std::size_t ieltype = 0; ieltype < m_mesh->bagElementTypeDomainBoundary().size(); ++ieltype) {
978 8 ElementType eltType = m_mesh->bagElementTypeDomainBoundary()[ieltype];
979 8 const GeoElement& geoEle = *GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
980
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 const RefElement& refEle = m_chooseRefEleFromFEType(geoEle, eltType);
981
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 CurrentFiniteElement fe(refEle, geoEle, DegreeOfExactness_0);
982
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_linkNodesGeoAndRefEle(fe, feIdLocalDof);
983
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 elemConnectivity.resize(refEle.numDOF(),0);
984
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 elem.resize(GeometricMeshRegion::m_numPointsPerElt[eltType], 0);
985
986 8 for(auto itRef = m_mesh->intRefToBegEndMaps[eltType].begin();
987
2/2
✓ Branch 4 taken 76 times.
✓ Branch 5 taken 8 times.
 84 itRef != m_mesh->intRefToBegEndMaps[eltType].end(); itRef++) {
988 76 std::size_t currentLabel = itRef->first;
989 76 std::size_t numElemsPerRef = itRef->second.second;
990
2/2
✓ Branch 0 taken 4608 times.
✓ Branch 1 taken 76 times.
 4684 for (std::size_t iel = 0; iel < numElemsPerRef; iel++ ) {
991
1/2
✓ Branch 1 taken 4608 times.
✗ Branch 2 not taken.
 4608 m_getElementSupportConnectivity(refEle, eltType, countEltPerType[eltType], feIdLocalDof, elem, elemConnectivity);
992
2/2
✓ Branch 1 taken 13440 times.
✓ Branch 2 taken 4608 times.
 18048 for (std::size_t iDof = 0; iDof < elemConnectivity.size(); iDof++) {
993
2/2
✓ Branch 0 taken 1132 times.
✓ Branch 1 taken 12308 times.
 13440 if(currentLabel == labelIN )
994
1/2
✓ Branch 2 taken 1132 times.
✗ Branch 3 not taken.
 1132 listSuppDofIN.insert(elemConnectivity[iDof]);
995
2/2
✓ Branch 0 taken 1132 times.
✓ Branch 1 taken 12308 times.
 13440 if(currentLabel == labelOUT )
996
1/2
✓ Branch 2 taken 1132 times.
✗ Branch 3 not taken.
 1132 listSuppDofOUT.insert(elemConnectivity[iDof]);
997 }
998 4608 cptElt++;
999 4608 countEltPerType[eltType]++;
1000 }
1001 }
1002 8 }
1003 8 }
1004
1005 8 void SupportDofMesh::GetPermutationList(std::vector <std::pair< std::size_t, std::size_t> >& listSuppDofMatched)
1006 {
1007 std::size_t iin, iio, iaux;
1008
1009
2/2
✓ Branch 1 taken 268 times.
✓ Branch 2 taken 8 times.
 276 for (std::size_t ip = 0; ip < listSuppDofMatched.size(); ++ip) {
1010
1011 // retrieve the SuppDof to be linked
1012 268 iin = listSuppDofMatched[ip].first; // SuppDofs to be kept
1013 268 iio = listSuppDofMatched[ip].second; // SuppDofs to be removed
1014
1015 // 1st step - remove iioSuppDof in m_listPerm
1016
2/2
✓ Branch 1 taken 125760 times.
✓ Branch 2 taken 268 times.
 126028 for(std::size_t jp = iio+1; jp < m_listPerm.size(); ++jp)
1017 125760 m_listPerm[jp]-=1;
1018
2/2
✓ Branch 0 taken 6756 times.
✓ Branch 1 taken 268 times.
 7024 for (std::size_t jp = 0; jp < ip; ++jp) {
1019 6756 iaux = listSuppDofMatched[jp].second;
1020
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 6756 times.
 6756 if ( m_listPerm[iaux] > m_listPerm[iio] )
1021 ✗ m_listPerm[iaux]-=1;
1022 }
1023
1024 // 2nd step - assign new link to SuppDof
1025 268 m_listPerm[iio]=m_listPerm[iin];
1026 }
1027 8 }
1028
1029 // Embedded interface
1030 ✗ void SupportDofMesh::matchEmbeddedLabelPairs(const Variable& variable)
1031 {
1032 ✗ std::size_t cpt = 0;
1033 std::size_t MatchNumLabPerVariable;
1034 double distNodes;
1035 ✗ std::vector< std::size_t > checkboard;
1036
1037 std::size_t labelSide0, labelSide1;
1038 std::size_t idElemSide0, idElemSide1;
1039 ✗ std::set<std::size_t> listSuppDofSide0, listSuppDofSide1;
1040 ✗ std::vector< std::pair<std::size_t, std::set<std::size_t> > > listSuppElemSide0, listSuppElemSide1;
1041
1042 ✗ for(std::size_t mvar = 0; mvar < FelisceParam::instance().EmbeddedVariable.size(); mvar++) {
1043 ✗ MatchNumLabPerVariable = FelisceParam::instance().EmbeddedNumInterface[mvar];
1044 ✗ if ( FelisceParam::instance().EmbeddedVariable[mvar] == variable.name().c_str()) {
1045 ✗ PetscPrintf(PETSC_COMM_WORLD,"##### For variable %s \n",variable.name().c_str() );
1046
1047 ✗ for(std::size_t mlab = 0; mlab < MatchNumLabPerVariable; mlab++) {
1048 ✗ labelSide0 = FelisceParam::instance().EmbeddedLabelPairs[cpt + 2*mlab];
1049 ✗ labelSide1 = FelisceParam::instance().EmbeddedLabelPairs[cpt + 2*mlab +1];
1050 ✗ PetscPrintf(PETSC_COMM_WORLD,"Match label %ld (the computations are perfomed on this label) and label %ld \n", static_cast<unsigned long>(labelSide0), static_cast<unsigned long>(labelSide1) );
1051 ✗ matchingElemConnectivity(labelSide0, labelSide1, listSuppDofSide0, listSuppDofSide1);
1052 ✗ extractIdElemWithConnectivityFromSides(labelSide0, labelSide1, listSuppElemSide0, listSuppElemSide1);
1053 }
1054
1055 // loop to match the listSuppDofSide0 with the listSuppDofSide1
1056 ✗ for (auto kDofSide1=listSuppDofSide1.begin(); kDofSide1 != listSuppDofSide1.end() ; ++kDofSide1 ) {
1057 ✗ Point& pSide1 = listNode()[*kDofSide1];
1058 ✗ for (auto kDofSide0=listSuppDofSide0.begin(); kDofSide0 != listSuppDofSide0.end() ; ++kDofSide0 ) {
1059 ✗ Point& pSide0 = listNode()[*kDofSide0];
1060 ✗ distNodes = std::sqrt( ( pSide0.x() - pSide1.x() ) * ( pSide0.x() - pSide1.x() ) +
1061 ✗ ( pSide0.y() - pSide1.y() ) * ( pSide0.y() - pSide1.y() ) +
1062 ✗ ( pSide0.z() - pSide1.z() ) * ( pSide0.z() - pSide1.z() ) );
1063 ✗ if ( distNodes < FelisceParam::instance().FusionTolerance ){
1064 ✗ m_mapSuppDofMatched0.insert(std::make_pair(*kDofSide0,*kDofSide1));
1065 ✗ m_mapSuppDofMatched1.insert(std::make_pair(*kDofSide1,*kDofSide0));
1066 }
1067 }
1068 }
1069
1070 // loop to match listSuppElemSide0 with the listSuppElemSide1
1071 ✗ for (std::size_t kElemSide0= 0; kElemSide0 < listSuppElemSide0.size(); kElemSide0++ ){
1072 ✗ idElemSide0 = listSuppElemSide0[kElemSide0].first;
1073 ✗ for (std::size_t kElemSide1= 0; kElemSide1 < listSuppElemSide1.size(); kElemSide1++ ){
1074 ✗ idElemSide1 = listSuppElemSide1[kElemSide1].first;
1075 // check wether dofs of the elements shared the same support point
1076 ✗ checkboard.clear();
1077 ✗ checkboard.resize((listSuppElemSide0[kElemSide0].second).size());
1078 ✗ int i = 0;
1079 ✗ for (auto kDofSide0=(listSuppElemSide0[kElemSide0].second).begin(); kDofSide0 != (listSuppElemSide0[kElemSide0].second).end() ; ++kDofSide0 ) {
1080 ✗ if( (listSuppElemSide1[kElemSide1].second).count(m_mapSuppDofMatched0[*kDofSide0]) )
1081 ✗ checkboard[i] = 1;
1082 else
1083 ✗ checkboard[i] = 0;
1084 ✗ ++i;
1085 }
1086 ✗ std::size_t count = 0;
1087 ✗ for(std::size_t check= 0; check < checkboard.size(); check++ )
1088 ✗ count += checkboard[check];
1089 ✗ if ( count == (listSuppElemSide0[kElemSide0].second).size() )
1090 ✗ m_mapSuppElemMatched.insert(std::make_pair(idElemSide0,idElemSide1));
1091 }
1092 }
1093 ✗ cpt += MatchNumLabPerVariable*2;
1094 }
1095 }
1096 ✗ m_createdEmbeddedInterfaceMaps = true;
1097 }
1098
1099 /// CRACK INTERFACE
1100 ✗ void SupportDofMesh::matchCrackLabelPairs(const Variable& variable)
1101 {
1102 ✗ std::size_t cpt = 0;
1103 std::size_t MatchNumLabPerVariable;
1104 // double distNodes;
1105 ✗ std::vector< std::size_t > checkboard;
1106
1107 std::size_t labelSide0, labelSide1;
1108 std::size_t idElemSide0, idElemSide1;
1109 ✗ std::set<std::size_t> listSuppDofSide0, listSuppDofSide1;
1110 ✗ std::vector< std::pair<std::size_t, std::set<std::size_t> > > listSuppElemSide0, listSuppElemSide1;
1111
1112 ✗ for(std::size_t mvar = 0; mvar < FelisceParam::instance().CrackVariable.size(); mvar++) {
1113 ✗ MatchNumLabPerVariable = FelisceParam::instance().CrackNum[mvar];
1114 ✗ if ( FelisceParam::instance().CrackVariable[mvar] == variable.name().c_str()) {
1115 ✗ PetscPrintf(PETSC_COMM_WORLD,"##### For variable %s \n",variable.name().c_str() );
1116
1117 ✗ for(std::size_t mlab = 0; mlab < MatchNumLabPerVariable; mlab++) {
1118 ✗ labelSide0 = FelisceParam::instance().CrackLabelPairs[cpt + 2*mlab];
1119 ✗ labelSide1 = FelisceParam::instance().CrackLabelPairs[cpt + 2*mlab +1];
1120 ✗ PetscPrintf(PETSC_COMM_WORLD,"Match label %ld (the computations are perfomed on this label) and label %ld \n", static_cast<unsigned long>(labelSide0), static_cast<unsigned long>(labelSide1) );
1121 ✗ matchingElemConnectivity(labelSide0, labelSide1, listSuppDofSide0, listSuppDofSide1);
1122 ✗ extractIdElemWithConnectivityFromSides(labelSide0, labelSide1, listSuppElemSide0, listSuppElemSide1);
1123 }
1124
1125 // loop to match the listSuppDofSide0 with the listSuppDofSide1
1126 ✗ for (auto kDofSide1=listSuppDofSide1.begin(); kDofSide1 != listSuppDofSide1.end() ; ++kDofSide1 ) {
1127 // Point& pSide1 = listNode()[*kDofSide1];
1128 ✗ for (auto kDofSide0=listSuppDofSide0.begin(); kDofSide0 != listSuppDofSide0.end() ; ++kDofSide0 ) {
1129 // Point& pSide0 = listNode()[*kDofSide0];
1130 ✗ m_mapSuppDofMatchedCrack0.insert(std::make_pair(*kDofSide0,*kDofSide1));
1131 ✗ m_mapSuppDofMatchedCrack1.insert(std::make_pair(*kDofSide1,*kDofSide0));
1132 }
1133 }
1134
1135 // loop to match listSuppElemSide0 with the listSuppElemSide1
1136 ✗ for (std::size_t kElemSide0= 0; kElemSide0 < listSuppElemSide0.size(); kElemSide0++ ){
1137 ✗ idElemSide0 = listSuppElemSide0[kElemSide0].first;
1138 ✗ for (std::size_t kElemSide1= 0; kElemSide1 < listSuppElemSide1.size(); kElemSide1++ ){
1139 ✗ idElemSide1 = listSuppElemSide1[kElemSide1].first;
1140 // check wether dofs of the elements shared the same support point
1141 ✗ checkboard.clear();
1142 ✗ checkboard.resize((listSuppElemSide0[kElemSide0].second).size());
1143 ✗ int i = 0;
1144 ✗ for (auto kDofSide0=(listSuppElemSide0[kElemSide0].second).begin(); kDofSide0 != (listSuppElemSide0[kElemSide0].second).end() ; ++kDofSide0 ) {
1145 ✗ if( (listSuppElemSide1[kElemSide1].second).count(m_mapSuppDofMatchedCrack0[*kDofSide0]) )
1146 ✗ checkboard[i] = 1;
1147 else
1148 ✗ checkboard[i] = 0;
1149 ✗ ++i;
1150 }
1151 ✗ std::size_t count = 0;
1152 ✗ for(std::size_t check= 0; check < checkboard.size(); check++ )
1153 ✗ count += checkboard[check];
1154 ✗ if ( count == (listSuppElemSide0[kElemSide0].second).size() )
1155 ✗ m_mapSuppElemMatchedCrack.insert(std::make_pair(idElemSide0,idElemSide1));
1156 }
1157 }
1158 ✗ cpt += MatchNumLabPerVariable*2;
1159 }
1160 }
1161 ✗ m_createdCrackInterfaceMaps = true;
1162 }
1163
1164 // both CRACK AND EMBEDDED
1165 ✗ void SupportDofMesh::extractIdElemWithConnectivityFromSides(std::size_t labelSide0, std::size_t labelSide1,
1166 std::vector< std::pair<std::size_t, std::set<std::size_t> > >& listSuppElemSide0,
1167 std::vector< std::pair<std::size_t, std::set<std::size_t> > >& listSuppElemSide1)
1168 {
1169 ✗ std::vector<int> elem, elemConnectivity, feIdLocalDof;
1170 ✗ std::size_t cptElt = 0;
1171 std::size_t countEltPerType[ GeometricMeshRegion::m_numTypesOfElement ];
1172 ✗ std::vector<felInt> vecSupport;
1173 ✗ std::set<std::size_t> listSuppDof;
1174
1175 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ )
1176 ✗ countEltPerType[ityp] = 0;
1177
1178 ✗ for (std::size_t ieltype = 0; ieltype < m_mesh->bagElementTypeDomainBoundary().size(); ++ieltype) {
1179 ✗ ElementType eltType = m_mesh->bagElementTypeDomainBoundary()[ieltype];
1180 ✗ const GeoElement& geoEle = *GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
1181 ✗ const RefElement& refEle = m_chooseRefEleFromFEType(geoEle, eltType);
1182 ✗ CurrentFiniteElement fe(refEle, geoEle, DegreeOfExactness_0);
1183 ✗ m_linkNodesGeoAndRefEle(fe, feIdLocalDof);
1184 ✗ elemConnectivity.resize(refEle.numDOF(),0);
1185 ✗ elem.resize(GeometricMeshRegion::m_numPointsPerElt[eltType], 0);
1186
1187 ✗ for(auto itRef = m_mesh->intRefToBegEndMaps[eltType].begin();
1188 ✗ itRef != m_mesh->intRefToBegEndMaps[eltType].end(); itRef++) {
1189 ✗ std::size_t currentLabel = itRef->first;
1190 ✗ std::size_t numElemsPerRef = itRef->second.second;
1191
1192 ✗ for (std::size_t iel = 0; iel < numElemsPerRef; iel++ ) {
1193
1194 ✗ m_getElementSupportConnectivity(refEle, eltType, countEltPerType[eltType], feIdLocalDof, elem, elemConnectivity);
1195 ✗ getIdElementSupport(eltType, countEltPerType[eltType], vecSupport);
1196
1197 ✗ for (std::size_t iSupportEl = 0; iSupportEl < vecSupport.size(); iSupportEl++) {
1198 ✗ if(currentLabel == labelSide0 ){
1199 ✗ listSuppDof.clear();
1200 ✗ for (std::size_t iDof = 0; iDof < elemConnectivity.size(); iDof++)
1201 ✗ listSuppDof.insert(elemConnectivity[iDof]);
1202 ✗ listSuppElemSide0.emplace_back(vecSupport[iSupportEl],listSuppDof);
1203 }
1204 ✗ if(currentLabel == labelSide1 ){
1205 ✗ listSuppDof.clear();
1206 ✗ for (std::size_t iDof = 0; iDof < elemConnectivity.size(); iDof++)
1207 ✗ listSuppDof.insert(elemConnectivity[iDof]);
1208 ✗ listSuppElemSide1.emplace_back(vecSupport[iSupportEl],listSuppDof);
1209 }
1210 }
1211 ✗ cptElt++;
1212 ✗ countEltPerType[eltType]++;
1213 }
1214 }
1215 }
1216 }
1217
1218 26 void SupportDofMesh::duplicateSupportElements( std::map<felInt, felInt>& listIntersectedEltIdx,
1219 std::map<felInt, std::vector<felInt> >& listIntersectedVerSgn,
1220 std::map<felInt, std::map<felInt, felInt> >& mapVerticesIdxToDupl
1221 )
1222 {
1223
1224 // add new point to listNode and listPerm and set the id of each node marked for duplication
1225 // TODO in case of fusionDof we do it here, but it should be done calling fusionSupportDof after the duplication...
1226 26 felInt curSupport = 0;
1227
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 auto& r_instance = FelisceParam::instance();
1228 26 auto& fusionVariable = r_instance.FusionVariable;
1229
2/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 7 taken 26 times.
✗ Branch 8 not taken.
 26 auto doFusionDof = std::find(fusionVariable.begin(), fusionVariable.end(), m_variable.name().c_str());
1230
1231
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
 26 if( doFusionDof != fusionVariable.end() ) {
1232
1233 ✗ if ( r_instance.contactDarcy > 0 ) {
1234 // TODO D.C.
1235 // -------------
1236 // If we are modifying the supportDofUnknow for the darcyPressure we fusion also the new duplicated node
1237 // on the first node with label = darcyLabel
1238 ✗ for(auto it_map = mapVerticesIdxToDupl.begin(); it_map != mapVerticesIdxToDupl.end(); ++it_map) {
1239
1240 ✗ auto& nest_map = it_map->second;
1241 ✗ for(auto it_nest_map = nest_map.begin(); it_nest_map != nest_map.end(); ++it_nest_map) {
1242
1243 ✗ m_listNode.push_back(m_listNode[it_map->first]);
1244 ✗ m_listPerm.push_back(0); // because we do fusionDof with the first Darcy-label node
1245
1246 ✗ it_nest_map->second = getIdFusionDarcy();
1247 }
1248 }
1249
1250 // we don't need to change the number of supportDof for the Darcy Pressure in case of fusionDof for Darcy,
1251 // TODO D.C.
1252 } else {
1253
1254 ✗ FEL_ERROR("Not implemented yet");
1255 }
1256 } else {
1257
1258 // add the duplicated support dof at the end of listNode and update listPerm
1259
2/2
✓ Branch 4 taken 970 times.
✓ Branch 5 taken 26 times.
 996 for(auto it_map = mapVerticesIdxToDupl.begin(); it_map != mapVerticesIdxToDupl.end(); ++it_map) {
1260
1261 970 auto& nest_map = it_map->second;
1262
2/2
✓ Branch 3 taken 970 times.
✓ Branch 4 taken 970 times.
 1940 for(auto it_nest_map = nest_map.begin(); it_nest_map != nest_map.end(); ++it_nest_map){
1263
1264
1/2
✓ Branch 3 taken 970 times.
✗ Branch 4 not taken.
 970 m_listNode.push_back(m_listNode[it_map->first]);
1265
1/2
✓ Branch 2 taken 970 times.
✗ Branch 3 not taken.
 970 m_listPerm.push_back(m_listNode.size());
1266
1267 970 it_nest_map->second = m_numSupportDof + curSupport;
1268 970 ++curSupport;
1269 }
1270 }
1271
1272 // change m_numSupportDof
1273 26 m_numSupportDof = m_listNode.size();
1274 }
1275
1276 // for each intersected element
1277
2/2
✓ Branch 4 taken 970 times.
✓ Branch 5 taken 26 times.
 996 for(auto itListInt = listIntersectedEltIdx.begin(); itListInt != listIntersectedEltIdx.end(); ++itListInt) {
1278
1279 // get the id of the element
1280 970 const felInt idElement = itListInt->first;
1281
1282 // get the id of the interface
1283 970 const felInt idInterface = itListInt->second;
1284
1285 // get the id of the support element
1286 970 const felInt idSupportElt = m_vectorIdElementSupport[idElement][0];
1287
1288 // number of support dof of the element
1289 970 const felInt numSupportDofOfElem = getNumSupportDof(idSupportElt);
1290
1291 // insert new support dof (duplicate them) in m_iEle
1292
1/2
✓ Branch 6 taken 970 times.
✗ Branch 7 not taken.
 970 m_iEle.insert(m_iEle.begin()+idSupportElt+2, m_iEle[idSupportElt] + numSupportDofOfElem);
1293
2/2
✓ Branch 1 taken 1932516 times.
✓ Branch 2 taken 970 times.
 1933486 for(std::size_t jel = idSupportElt+2; jel < m_iEle.size(); ++jel)
1294 1932516 m_iEle[jel] += numSupportDofOfElem;
1295
1296 // add the new support element in vectorIdElementSupport
1297
2/2
✓ Branch 1 taken 1931546 times.
✓ Branch 2 taken 970 times.
 1932516 for(std::size_t jel = idElement+1; jel < m_vectorIdElementSupport.size(); ++jel)
1298
2/2
✓ Branch 2 taken 1931546 times.
✓ Branch 3 taken 1931546 times.
 3863092 for(std::size_t idSupport = 0; idSupport < m_vectorIdElementSupport[jel].size(); ++idSupport)
1299 1931546 m_vectorIdElementSupport[jel][idSupport]++;
1300
1/2
✓ Branch 2 taken 970 times.
✗ Branch 3 not taken.
 970 m_vectorIdElementSupport[idElement].push_back(idSupportElt+1);
1301
1302 // update m_iSupportDof, duplicate the support element and fill it with -1
1303
1/2
✓ Branch 5 taken 970 times.
✗ Branch 6 not taken.
 970 m_iSupportDof.insert(m_iSupportDof.begin()+m_iEle[idSupportElt+1], numSupportDofOfElem, -1);
1304
1305 // for each support dof
1306
2/2
✓ Branch 0 taken 2858 times.
✓ Branch 1 taken 970 times.
 3828 for(felInt iSupport = 0; iSupport < numSupportDofOfElem; ++iSupport) {
1307
1308 // look for the current dof in the list of duplicated vertices
1309
1/2
✓ Branch 3 taken 2858 times.
✗ Branch 4 not taken.
 2858 const auto ver_2_map = mapVerticesIdxToDupl.find( m_iSupportDof[m_iEle[idSupportElt]+iSupport] );
1310
1311 // duplicate the support if needed
1312
2/2
✓ Branch 2 taken 2806 times.
✓ Branch 3 taken 52 times.
 2858 if( ver_2_map != mapVerticesIdxToDupl.end() ) {
1313
1314 // get the duplicated vertex
1315
1/2
✓ Branch 2 taken 2806 times.
✗ Branch 3 not taken.
 2806 const auto itf_2_dup = ver_2_map->second.find( idInterface );
1316
1317
1/2
✓ Branch 3 taken 2806 times.
✗ Branch 4 not taken.
 2806 if ( itf_2_dup != ver_2_map->second.end() ) {
1318
1319 // get the side
1320
1/2
✓ Branch 1 taken 2806 times.
✗ Branch 2 not taken.
 2806 const auto side = listIntersectedVerSgn[idElement][iSupport];
1321
1322
2/2
✓ Branch 0 taken 1420 times.
✓ Branch 1 taken 1386 times.
 2806 if( side > 0 ) {
1323 // Point on the left part of the structure
1324 // set only the duplicated support dof
1325 1420 m_iSupportDof[m_iEle[idSupportElt+1]+iSupport] = itf_2_dup->second;
1326
1/2
✓ Branch 0 taken 1386 times.
✗ Branch 1 not taken.
 1386 } else if ( side < 0 ) {
1327 // Point on the right part of the structure
1328 // first set the duplicated support dof with the original one and second change the original.
1329 1386 m_iSupportDof[m_iEle[idSupportElt+1]+iSupport] = m_iSupportDof[m_iEle[idSupportElt]+iSupport];
1330 1386 m_iSupportDof[m_iEle[idSupportElt]+iSupport] = itf_2_dup->second;
1331 } else {
1332 // degenerated case, the mesh vertex is on the interface
1333 ✗ m_iSupportDof[m_iEle[idSupportElt+1]+iSupport] = itf_2_dup->second;
1334 }
1335 } else {
1336
1337 ✗ FEL_ERROR("Not really clear how the duplication should be in this case");
1338 }
1339 } else {
1340 // copy the one in the original support element
1341 52 m_iSupportDof[m_iEle[idSupportElt+1]+iSupport] = m_iSupportDof[m_iEle[idSupportElt]+iSupport];
1342 }
1343 }
1344 }
1345 26 }
1346
1347 //--------------------------------------
1348 ✗ void SupportDofMesh::getIndexesElementSupport(int iel, std::vector<int>& indexes) const
1349 {
1350 ✗ FEL_ASSERT_LT(iel, m_iEle.size() - 1);
1351 ✗ FEL_ASSERT(indexes.empty());
1352
1353 ✗ felInt first_index = m_iEle[iel];
1354 ✗ felInt last_index = m_iEle[iel + 1] - 1;
1355
1356 ✗ indexes.reserve(last_index - first_index + 1);
1357
1358 ✗ for (felInt i = first_index; i <= last_index; ++i) {
1359 ✗ FEL_ASSERT(i < static_cast<felInt>(m_iSupportDof.size()));
1360 ✗ indexes.push_back(m_iSupportDof[i]);
1361 }
1362 }
1363 //--------------------------------------
1364
1365 ✗ void SupportDofMesh::printIdElementSupport(std::ostream& outstr) const
1366 {
1367 ElementType eltType;
1368 felInt numElemsPerRef;
1369 felInt countEltPerType[ GeometricMeshRegion::m_numTypesOfElement ];
1370 ✗ std::vector<felInt> vecSupport;
1371
1372 ✗ for (int ityp = 0; ityp < GeometricMeshRegion::m_numTypesOfElement; ityp++) {
1373 ✗ eltType = (ElementType)ityp;
1374 ✗ countEltPerType[eltType] = 0;
1375 }
1376 ✗ const std::vector<ElementType>& bagElementTypeDomainBoundary = m_mesh->bagElementTypeDomainBoundary();
1377 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_mesh->bagElementTypeDomain();
1378
1379 ✗ outstr << "Mapping eltGeom to eltSupport in supportDof" << std::endl;
1380
1381 ✗ for (std::size_t ityp = 0; ityp < bagElementTypeDomain.size() ; ++ityp) {
1382 ✗ eltType = bagElementTypeDomain[ityp];
1383 ✗ if (m_mesh->numElements(eltType) != 0 ) {
1384 ✗ for(auto itRef = m_mesh->intRefToBegEndMaps[eltType].begin();
1385 ✗ itRef != m_mesh->intRefToBegEndMaps[eltType].end(); itRef++) {
1386 ✗ numElemsPerRef = itRef->second.second;
1387 ✗ for ( felInt iel = 0; iel < numElemsPerRef; iel++ ) {
1388 ✗ getIdElementSupport(eltType, countEltPerType[eltType], vecSupport);
1389 ✗ outstr << "( " << GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].first << ", " << iel << ")--->";
1390 ✗ for(std::size_t idSupport=0; idSupport<vecSupport.size(); ++idSupport)
1391 ✗ outstr << " " << vecSupport[idSupport];
1392 ✗ outstr << std::endl;
1393 ✗ countEltPerType[eltType]++;
1394 }
1395 }
1396 }
1397 }
1398
1399 ✗ for (std::size_t ityp = 0; ityp < bagElementTypeDomainBoundary.size() ; ++ityp) {
1400 ✗ eltType = bagElementTypeDomainBoundary[ityp];
1401 ✗ if (m_mesh->numElements(eltType) != 0 ) {
1402 ✗ for(auto itRef = m_mesh->intRefToBegEndMaps[eltType].begin();
1403 ✗ itRef != m_mesh->intRefToBegEndMaps[eltType].end(); itRef++) {
1404 ✗ numElemsPerRef = itRef->second.second;
1405 ✗ for ( felInt iel = 0; iel < numElemsPerRef; iel++ ) {
1406 ✗ getIdElementSupport(eltType, countEltPerType[eltType], vecSupport);
1407 ✗ outstr << "( " << GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].first << ", " << iel << ")--->";
1408 ✗ for(std::size_t idSupport=0; idSupport<vecSupport.size(); ++idSupport)
1409 ✗ outstr << " " << vecSupport[idSupport];
1410 ✗ outstr << std::endl;
1411 ✗ countEltPerType[eltType]++;
1412 }
1413 }
1414 }
1415 }
1416 }
1417
1418 /*!
1419 \brief Print the support of dof in CSR format
1420 \param[in] verbose Level of verbose
1421 \param[in] outstr The out stream
1422 */
1423 858 void SupportDofMesh::print(int verbose, std::ostream& outstr) const
1424 {
1425 int rank;
1426
1/2
✓ Branch 1 taken 858 times.
✗ Branch 2 not taken.
 858 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
1427
1428
1/2
✓ Branch 0 taken 858 times.
✗ Branch 1 not taken.
 858 if (verbose) {
1429
1/2
✓ Branch 2 taken 858 times.
✗ Branch 3 not taken.
 858 PetscPrintf(PETSC_COMM_WORLD,"[%d] Number of DOF supports: %d\n", rank, numSupportDof());
1430 }
1431
1432
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 858 times.
 858 if (verbose > 24) {
1433 ✗ outstr << "\n [" << rank << "] Nodes coordinates (size = " << m_listNode.size() <<"): \n";
1434 ✗ for (felInt iv = 0; iv < static_cast<felInt>(m_listNode.size()); iv++) {
1435 ✗ outstr << "[" << rank << "] ( " << m_listNode[iv].x() << ", " << m_listNode[iv].y() << ", " << m_listNode[iv].z() << ") "<< std::endl;
1436 }
1437 ✗ outstr << std::endl;
1438 ✗ outstr << "[" << rank << "] Supports of dof associated to elements (CSR format):\n";
1439 ✗ outstr << "[" << rank << "] The i th element have iEle[i+1] - iEle[i] supports of dof.\n";
1440 ✗ outstr << "[" << rank << "] iEle: ";
1441 ✗ for (felInt iel = 0; iel < static_cast<felInt>(m_iEle.size()); iel++) {
1442 ✗ outstr << m_iEle[iel] << " ";
1443 }
1444 ✗ outstr << std::endl << std::endl;
1445
1446 ✗ outstr << "[" << rank << "] Supports of dof of the i th element are stored from iSupportDof[ iEle[i] to (iSupportDof[ iEle[i+1] - 1):"<< std::endl;
1447 ✗ outstr << "[" << rank << "] iSupportDof: ";
1448 ✗ for (unsigned int iSupDof = 0; iSupDof < m_iSupportDof.size(); iSupDof++) {
1449 ✗ outstr << m_iSupportDof[iSupDof] << " ";
1450 }
1451 ✗ outstr << std::endl;
1452
1453 ✗ outstr << " m_vectorIdElementSupport is : " << std::endl;
1454 ✗ for (unsigned int iSupDof = 0; iSupDof < m_vectorIdElementSupport.size(); iSupDof++) {
1455 ✗ for (unsigned int jSupDof = 0; jSupDof < m_vectorIdElementSupport[iSupDof].size(); jSupDof++)
1456 ✗ outstr << m_vectorIdElementSupport[iSupDof][jSupDof] << " ";
1457 ✗ outstr << std::endl;
1458 }
1459 ✗ outstr << std::endl;
1460
1461 /* if (FelisceParam::instance().FusionDof) // to better arrange
1462 {
1463 std::cout<< "\n/======= Periodic boundary conditions. ======/\n"<< std::endl;
1464 std::cout<<"\nLabel IN: "<<std::endl;
1465 std::cout<<"Support Dof:"<<std::endl;
1466 for(int iin=0; iin<m_listSuppDofIN.size(); iin++){
1467 std::cout<<m_listSuppDofIN[iin]<< " ";
1468 }
1469 std::cout<<std::endl;
1470
1471
1472 std::cout<<"\nLabel OUT: " <<std::endl;
1473 std::cout<<"Support Dof:"<<std::endl;
1474 for(int iout=0; iout<m_listSuppDofOUT.size(); iout++){
1475 std::cout<<m_listSuppDofOUT[iout]<< " ";
1476 }
1477 std::cout<<std::endl;
1478 }*/
1479
1480
1481 ✗ if (FelisceParam::instance().EmbeddedInterface && m_createdEmbeddedInterfaceMaps){
1482 ✗ std::cout << "m_mapSuppElemMatched: {";
1483 ✗ for(auto idElem = m_mapSuppElemMatched.begin(); idElem != m_mapSuppElemMatched.end(); ++idElem)
1484 ✗ std::cout << " {" << idElem->first << " " << idElem->second << "}";
1485 ✗ std::cout << " }" << std::endl;
1486
1487 ✗ std::cout << "m_mapSuppDofMatched0: {";
1488 ✗ for(auto idDof = m_mapSuppDofMatched0.begin(); idDof != m_mapSuppDofMatched0.end(); ++idDof)
1489 ✗ std::cout << " {" << idDof->first << " " << idDof->second << "}";
1490 ✗ std::cout << " }" << std::endl;
1491
1492 ✗ std::cout << "m_mapSuppDofMatched1: {";
1493 ✗ for(auto idDof = m_mapSuppDofMatched1.begin(); idDof != m_mapSuppDofMatched1.end(); ++idDof)
1494 ✗ std::cout << " {" << idDof->first << " " << idDof->second << "}";
1495 ✗ std::cout << " }" << std::endl;
1496 }
1497
1498 }
1499 858 }
1500 }
1501