GCC Code Coverage Report

Directory: ./
File: Model/immersed.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 135 260 51.9%
Branches: 110 476 23.1%
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: M. Fernandez & L. Boilevin-Kayl
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Model/immersed.hpp"
21 #include "InputOutput/io.hpp"
22
23 namespace felisce
24 {
25 42 Immersed::Immersed() {
26
27 42 m_intfDisp = nullptr;
28 42 m_intfVelo = nullptr;
29 42 m_intfForc = nullptr;
30 42 m_velFluid = nullptr;
31 42 m_indexTime = 0;
32
33
1/2
✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
 42 m_allowProjection = FelisceParam::instance().allowProjectionImmersedStruct;
34 42 }
35
36 42 Immersed::~Immersed() = default;
37
38 8 void Immersed::initialize(GeometricMeshRegion::Pointer backMesh){
39
40 8 m_backMesh = backMesh;
41
42
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 6 times.
 8 if ( MpiInfo::rankProc() == 0 ) { // only the first processor possesses the full immersed mesh and will be responsible for the interpolation
43
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_readInterfaceMesh(); // read the immersed structure mesh and put it into m_interfMesh
44
45
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 const int verbose = FelisceParam::verbose();
46
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 const double tolrescaling = FelisceParam::instance().Geo.tolrescaling;
47
48
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 m_interpolator = felisce::make_shared<Interpolator>(m_backMesh.get(), m_allowProjection, tolrescaling, verbose);
49
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 m_intersector = felisce::make_shared<Intersector>(m_backMesh.get(), m_interfMesh.get(), tolrescaling, verbose);
50 }
51
52 8 MPI_Bcast(&m_nbPoints, 1, MPI_INT, 0, MpiInfo::petscComm());
53 8 MPI_Bcast(&m_nbCoor, 1, MPI_INT, 0, MpiInfo::petscComm());
54 8 }
55
56 8 void Immersed::setInterfaceExchangedData(
57 std::vector<double>* intfDisp,
58 std::vector<double>* intfVelo,
59 std::vector<double>* intfForc,
60 std::vector<double>* velFluid
61 ) {
62
63 8 m_intfDisp = intfDisp;
64 8 m_intfVelo = intfVelo;
65 8 m_intfForc = intfForc;
66 8 m_velFluid = velFluid;
67
68
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_intfDisp->resize(m_nbCoor*m_nbPoints,0.);
69
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_intfVelo->resize(m_nbCoor*m_nbPoints,0.);
70 8 }
71
72 320 void Immersed::updateInterface(GeometricMeshRegion& mesh) {
73
74
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 320 times.
 320 FEL_CHECK(m_nbPoints != 0,"Empty list of immersed points");
75
76
2/2
✓ Branch 1 taken 80 times.
✓ Branch 2 taken 240 times.
 320 if ( MpiInfo::rankProc() == 0 ) {
77 // WARNING: the interpolation only works so far for tetrahedral (in 3D) and triangular (in 2D) meshes
78
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 78 times.
 80 if ( m_nbVertices == 0 ) {
79 2 m_nbVertices = mesh.domainDim()+1; // number of coordinates = dimension + 1
80 }
81
82 // update interface points
83
2/2
✓ Branch 0 taken 1750 times.
✓ Branch 1 taken 80 times.
 1830 for (int l = 0; l < m_nbPoints; ++l) {
84
2/2
✓ Branch 0 taken 3500 times.
✓ Branch 1 taken 1750 times.
 5250 for (int ic = 0; ic < m_nbCoor; ++ic)
85 3500 m_listPoints[l].coor(ic) = m_listPointsRef[l].coor(ic) + (*m_intfDisp)[ic+l*m_nbCoor];
86 }
87
88 // update the mesh with m_intfDisp
89
1/2
✓ Branch 2 taken 80 times.
✗ Branch 3 not taken.
 80 m_interfMesh->moveMesh(*m_intfDisp, 1); // move the immersed structure mesh
90 80 m_buildLumpedMassMatrixOfStructure(); // update the lumped mass matrix of the immersed structure mesh
91
92 // interpolation of the immersed structure mesh
93 80 m_interpolator->interpolate(m_indElt, m_intOp_val, m_listPoints);
94 }
95 320 }
96
97 200 void Immersed::updateIntersection(GeometricMeshRegion& mesh, PetscVector& delta, const MPI_Comm& comm) {
98
99 // TODO remove mesh and use m_backMesh
100
101 // only first process will modify the value of the PetscVector delta and will write the files "intersection"
102
2/2
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 150 times.
 200 if ( MpiInfo::rankProc() == 0 ) {
103
104 // locate the elements near the structure interface and the ones in star of the previous elements where the SUPG-PSPG and div-div stabilization will be modified
105
1/2
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
 50 std::vector<felInt> elementID(m_interfMesh->numPoints());
106 50 std::vector<GeometricMeshRegion::seedElement> elementIDstar;
107
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 double coef = 1./FelisceParam::instance().localDivDivPenalty;
108
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 delta.set(1.0);
109
1/2
✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
 50 m_intersector->intersectMeshes(m_indElt, elementID);
110
111
2/2
✓ Branch 1 taken 1100 times.
✓ Branch 2 taken 50 times.
 1150 for (std::size_t i=0; i < elementID.size(); ++i) {
112
1/2
✓ Branch 2 taken 1100 times.
✗ Branch 3 not taken.
 1100 delta.setValue(elementID[i], coef, INSERT_VALUES);
113
1/2
✓ Branch 4 taken 1100 times.
✗ Branch 5 not taken.
 1100 mesh.getElementPatch(mesh.bagElementTypeDomain()[0], elementID[i], elementIDstar);
114
2/2
✓ Branch 1 taken 13500 times.
✓ Branch 2 taken 1100 times.
 14600 for (std::size_t j=0; j < elementIDstar.size(); ++j) {
115
1/2
✓ Branch 2 taken 13500 times.
✗ Branch 3 not taken.
 13500 delta.setValue(elementIDstar[j].second, coef, INSERT_VALUES);
116 }
117 }
118
119
2/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 50 times.
 50 if ( FelisceParam::instance().divDivStabPostProcess ) {
120 ✗ m_divDivStabPostProcess(elementID);
121 }
122
123 50 ++m_indexTime;
124 50 }
125
126 // once the PetscVector delta has been modified by the main process, broadcast it to all other processes of the communicator
127 200 delta.broadCastSequentialVector(comm, 0);
128 200 }
129
130 320 void Immersed::applyPenaltyToMatrixAndRHS(AO& ao, Dof& dof, int idVarVel, PetscMatrix& matrix, PetscVector& rhs, bool only_rhs) {
131
132 double lambda_i, lambda_j;
133
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 double eps = FelisceParam::instance().penaltyCoeffImmersedStruct;
134 int k, ig, jg;
135 320 double val = 0.0, coef = 1.0;
136
137
3/4
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 80 times.
✓ Branch 4 taken 240 times.
 320 if ( MpiInfo::rankProc() == 0 ) { // only the first processor add the "penalized" terms to the NS matrix
138
139 // matrix and rhs contributions for NS matrix
140
2/2
✓ Branch 0 taken 1750 times.
✓ Branch 1 taken 80 times.
 1830 for (int l=0; l< m_nbPoints; ++l) {
141 1750 k = m_indElt[l]; // the element containing the point
142
1/2
✓ Branch 0 taken 1750 times.
✗ Branch 1 not taken.
 1750 if ( k != -1 ) { // if the interpolation of the immersed structure mesh into the fluid mesh has been successful
143
144 // implicit method
145
2/4
✓ Branch 1 taken 1750 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1750 times.
✗ Branch 4 not taken.
 1750 if ( FelisceParam::instance().fsiRNscheme == -1000 ) {
146
2/4
✓ Branch 1 taken 1750 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1750 times.
 1750 if ( FelisceParam::instance().useMassMatrix == 1 ) { // with use of the lumped mass matrix in the fourth block of the fluid problem matrix
147 ✗ coef = 1.0/(eps*m_lumpedMassMatrix[l]);
148 }
149 else { // with only penalisation coefficient in the fourth block of the fluid problem matrix
150 1750 coef = 1.0/eps;
151 }
152 }
153
154 // explicit method without penalisation coefficient
155 ✗ else if ( FelisceParam::instance().fsiRNscheme == 1 ) {
156 ✗ coef = (m_lumpedMassMatrix[l]*(FelisceParam::instance().densitySolid)*(FelisceParam::instance().thickness))/(FelisceParam::instance().timeStep);
157 }
158
159
160 // if the specified fsiRNscheme has not been recognized, returns an error
161 else {
162 ✗ std::ostringstream tmp;
163 ✗ tmp << FelisceParam::instance().fsiRNscheme;
164 ✗ std::string fsiRNscheme_as_string(tmp.str());
165 ✗ FEL_ERROR(("Immersed::computePenaltyForce: value of fsiRNscheme = " + fsiRNscheme_as_string + " not recognized! STOP HERE.\n"));
166 }
167
168 //
169
2/2
✓ Branch 0 taken 5250 times.
✓ Branch 1 taken 1750 times.
 7000 for (int i=0; i < m_nbVertices; ++i) {
170 5250 lambda_i = m_intOp_val[m_nbVertices*l+i ];
171
172
1/2
✓ Branch 0 taken 5250 times.
✗ Branch 1 not taken.
 5250 if ( not only_rhs ) {
173 // matrix
174
2/2
✓ Branch 0 taken 15750 times.
✓ Branch 1 taken 5250 times.
 21000 for (int j=0; j < m_nbVertices; ++j) {
175 15750 lambda_j = m_intOp_val[m_nbVertices*l+j ];
176 15750 val = coef*lambda_i*lambda_j; ///(1 + eps*coef);
177
2/2
✓ Branch 0 taken 31500 times.
✓ Branch 1 taken 15750 times.
 47250 for (int ic=0; ic < m_nbCoor; ++ic){
178
1/2
✓ Branch 1 taken 31500 times.
✗ Branch 2 not taken.
 31500 dof.loc2glob(k, i, idVarVel, ic,ig);
179
1/2
✓ Branch 1 taken 31500 times.
✗ Branch 2 not taken.
 31500 dof.loc2glob(k, j, idVarVel, ic,jg);
180
1/2
✓ Branch 1 taken 31500 times.
✗ Branch 2 not taken.
 31500 AOApplicationToPetsc(ao,1,&ig);
181
1/2
✓ Branch 1 taken 31500 times.
✗ Branch 2 not taken.
 31500 AOApplicationToPetsc(ao,1,&jg);
182
1/2
✓ Branch 1 taken 31500 times.
✗ Branch 2 not taken.
 31500 matrix.setValues(1,&ig,1,&jg,&val,ADD_VALUES);
183 }
184 }
185 }
186
187 // rhs
188
2/2
✓ Branch 0 taken 10500 times.
✓ Branch 1 taken 5250 times.
 15750 for (int ic=0; ic < m_nbCoor; ++ic) {
189
190 // implicit method
191
2/4
✓ Branch 1 taken 10500 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10500 times.
✗ Branch 4 not taken.
 10500 if ( FelisceParam::instance().fsiRNscheme == -1000 ) {
192 10500 val = coef*lambda_i*(*m_intfVelo)[ic + l*m_nbCoor ];
193 }
194
195 // explicit method without penalisation coefficient
196 ✗ else if ( (FelisceParam::instance().fsiRNscheme == 1) || (FelisceParam::instance().fsiRNscheme == -1) ) {
197 ✗ val = coef*lambda_i*(*m_intfVelo)[ic + l*m_nbCoor ] - lambda_i*m_intfForcOld[ic + l*m_nbCoor ];
198 }
199
200 // if the specified fsiRNscheme has not been recognized, returns an error
201 else {
202 ✗ std::ostringstream tmp;
203 ✗ tmp << FelisceParam::instance().fsiRNscheme;
204 ✗ std::string fsiRNscheme_as_string(tmp.str());
205 ✗ FEL_ERROR(("Immersed::computePenaltyForce: value of fsiRNscheme = " + fsiRNscheme_as_string + " not recognized! STOP HERE.\n"));
206 }
207
208
1/2
✓ Branch 1 taken 10500 times.
✗ Branch 2 not taken.
 10500 dof.loc2glob(k,i, idVarVel, ic,ig);
209
1/2
✓ Branch 1 taken 10500 times.
✗ Branch 2 not taken.
 10500 AOApplicationToPetsc(ao,1,&ig);
210
1/2
✓ Branch 1 taken 10500 times.
✗ Branch 2 not taken.
 10500 rhs.setValues(1,&ig,&val,ADD_VALUES);
211 }
212 }
213 }
214 }
215 } // end of MpiInfo::rankProc() == 0
216
217
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 matrix.assembly(MAT_FINAL_ASSEMBLY);
218
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 rhs.assembly();
219 320 }
220
221 ✗ void Immersed::computePenaltyForce(AO& ao, Dof& dof, int idVarVel, PetscVector& seqVel) {
222
223 double lambda_i;
224 ✗ double eps = FelisceParam::instance().penaltyCoeffImmersedStruct;
225
226 int k, ig;
227 ✗ double val = 0.0, coef = 1.0;
228 ✗ std::vector<double> aux(m_nbPoints*m_nbCoor);
229
230 ✗ if ( MpiInfo::rankProc() == 0 ) { // only the first processor add the "penalized" terms to the NS matrix
231
232 // compute
233 ✗ for (int l=0; l< m_nbPoints; ++l) {
234 ✗ k= m_indElt[l]; // the element containing the point
235 ✗ if ( k != -1 ) {
236
237 // implicit method
238 ✗ if ( FelisceParam::instance().fsiRNscheme == -1000 ) {
239 ✗ if ( FelisceParam::instance().useMassMatrix == 1 ) { // with use of the lumped mass matrix in the fourth block of the fluid problem matrix
240 ✗ coef = 1.0/(eps*m_lumpedMassMatrix[l]);
241 }
242 else { // with only penalisation coefficient in the fourth block of the fluid problem matrix
243 ✗ coef = 1.0/eps;
244 }
245 }
246
247 // explicit method without penalisation coefficient
248 ✗ else if ( FelisceParam::instance().fsiRNscheme == 1 ) {
249 ✗ coef = (m_lumpedMassMatrix[l]*(FelisceParam::instance().densitySolid)*(FelisceParam::instance().thickness))/(FelisceParam::instance().timeStep);
250 }
251
252 // if the specified fsiRNscheme has not been recognized, returns an error
253 else {
254 ✗ std::ostringstream tmp;
255 ✗ tmp << FelisceParam::instance().fsiRNscheme;
256 ✗ std::string fsiRNscheme_as_string(tmp.str());
257 ✗ FEL_ERROR(("Immersed::computePenaltyForce: value of fsiRNscheme = " + fsiRNscheme_as_string + " not recognized! STOP HERE.\n"));
258 }
259
260 ✗ for (int ic=0; ic < m_nbCoor; ++ic) {
261 ✗ (*m_velFluid)[l*m_nbCoor + ic] = 0.0;
262
263 ✗ for (int i=0; i < m_nbVertices; ++i) {
264 ✗ lambda_i = m_intOp_val[m_nbVertices*l+i ];
265 ✗ dof.loc2glob(k,i, idVarVel, ic,ig);
266 ✗ AOApplicationToPetsc(ao,1,&ig);
267 ✗ seqVel.getValues(1,&ig,&val);
268 ✗ (*m_velFluid)[l*m_nbCoor + ic] += val*lambda_i;
269 }
270
271 // implicit method
272 ✗ if ( FelisceParam::instance().fsiRNscheme == -1000 ) {
273 ✗ (*m_intfForc)[l*m_nbCoor + ic] = coef*((*m_velFluid)[l*m_nbCoor + ic] - (*m_intfVelo)[l*m_nbCoor + ic]);
274 }
275
276 // explicit method without penalisation coefficient
277 ✗ else if ( FelisceParam::instance().fsiRNscheme == 1 ) {
278 ✗ (*m_intfForc)[l*m_nbCoor + ic] = coef*((*m_velFluid)[l*m_nbCoor + ic] - (*m_intfVelo)[l*m_nbCoor + ic]) + m_intfForcOld[l*m_nbCoor + ic] ;
279 }
280
281 // if the specified fsiRNscheme has not been recognized, returns an error
282 else {
283 ✗ std::ostringstream tmp;
284 ✗ tmp << FelisceParam::instance().fsiRNscheme;
285 ✗ std::string fsiRNscheme_as_string(tmp.str());
286 ✗ FEL_ERROR(("Immersed::computePenaltyForce: value of fsiRNscheme = " + fsiRNscheme_as_string + " not recognized! STOP HERE.\n"));
287 }
288 }
289 }
290
291 else {
292
293 ✗ std::ostringstream tmp;
294 ✗ tmp << l;
295 ✗ std::string l_as_string(tmp.str());
296 ✗ FEL_ERROR(("Immersed::computePenaltyForce: it seems that the interpolation table met some problems for localising the structure node of label " + l_as_string + ", hence the m_intForc std::vector has been std::set to 0! STOP HERE.\n"));
297 }
298 }
299
300 // explicit method without penalisation coefficient
301 ✗ if ( FelisceParam::instance().fsiRNscheme == 1 ) {
302 ✗ m_intfForcOld = (*m_intfForc);
303 }
304 }
305 }
306
307 2 void Immersed::m_readInterfaceMesh() {
308
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 std::string inputDirectory1 = "./";
309
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 std::string inputFile1 = "";
310
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::string inputMesh1 = FelisceParam::instance().meshFileImmersedStruct;
311
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 std::string outputMesh1 = "toto.geo";
312
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::string meshDir1 = FelisceParam::instance().meshDirImmersedStruct;
313
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 std::string resultDir1 = "";
314
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 std::string prefixName1 = "";
315
316
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 IO io(inputDirectory1, inputFile1, inputMesh1, outputMesh1, meshDir1, resultDir1, prefixName1);
317
318
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_interfMesh = felisce::make_shared<GeometricMeshRegion>();
319
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 io.readMesh(*m_interfMesh, 1.0);
320
321 2 std::vector<Point>& listPoints = m_interfMesh->listPoints();
322 2 m_nbPoints = listPoints.size();
323
324
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_listPointsRef = listPoints;
325
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_listPoints = listPoints;
326
327 2 m_nbCoor = m_interfMesh->domainDim()+1; // number of coordinates = dimension + 1 // TODO D.C. doesn't work for .geo mesh
328
329 2 m_nbVertices = 0; // set in method updateInterface
330
331 // whent the method used is the implicit one with use of the lumped mass matrix in the fourth block of the fluid problem matrix OR the explicit method
332
7/14
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2 times.
 2 if ( (( FelisceParam::instance().fsiRNscheme == -1000 ) && ( FelisceParam::instance().useMassMatrix == 1 ) ) || ( FelisceParam::instance().fsiRNscheme == 1 ) ) { //
333 ✗ m_buildLumpedMassMatrixOfStructure(); // build the lumped mass matrix of the immersed structure mesh
334
335 ✗ if ( FelisceParam::instance().fsiRNscheme == 1) { // if explicit method
336 ✗ m_intfForcOld.resize(m_nbCoor*m_nbPoints, 0); // std::vector to store the previous force required for the computation of the new force in the explicit scheme
337 }
338 }
339 2 }
340
341 80 void Immersed::m_buildLumpedMassMatrixOfStructure() {
342
343 // The purpose of this method is to build the lumped mass matrix of the structure used into the semi-implicit/explicit scheme
344
345
2/4
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 80 times.
 80 if ( FelisceParam::verbose() > 1 ) {
346 ✗ std::cout << " -- Building of the lumped mass matrix of the structure..." << std::endl;
347 }
348
349 // initialization of the lumped mass matrix of the structure
350
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 m_lumpedMassMatrix.resize(m_nbPoints);
351
2/2
✓ Branch 0 taken 1750 times.
✓ Branch 1 taken 80 times.
 1830 for(int l = 0; l < m_nbPoints; ++l) {
352 1750 m_lumpedMassMatrix[l] = 0.;
353 }
354
355 typedef felisce::GeometricMeshRegion::ElementType ElementType;
356 ElementType eltType; // geometric element type in the mesh.
357 80 int numPointsPerElt = 0; // number of points per geometric element.
358 80 std::vector<Point*> elemPoint; // points of the current element.
359 80 std::vector<felInt> elemIdPoint; // id of points of the element.
360
361 // first loop on geometry type.
362 80 const std::vector<GeometricMeshRegion::ElementType>& bagElementTypeDomain = m_interfMesh->bagElementTypeDomain(); // list of domain element types in the mesh
363
2/2
✓ Branch 1 taken 80 times.
✓ Branch 2 taken 80 times.
 160 for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) { // for each domain element type
364 80 eltType = bagElementTypeDomain[i];
365
366 // resize array.
367 80 numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
368
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 elemPoint.resize(numPointsPerElt, nullptr);
369
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 elemIdPoint.resize(numPointsPerElt, 0);
370
371 80 int typeOfFiniteElement = 0;
372 80 const GeoElement *geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
373
1/2
✓ Branch 2 taken 80 times.
✗ Branch 3 not taken.
 80 const RefElement* refEle = geoEle->defineFiniteEle(eltType, typeOfFiniteElement, *m_interfMesh);
374
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 CurvilinearFiniteElement fe(*refEle, *geoEle, DegreeOfExactness_0);
375
376 // second loop on the number of elements for each geometry type.
377
2/2
✓ Branch 2 taken 1670 times.
✓ Branch 3 taken 80 times.
 1750 for(int j = 0; j<m_interfMesh->numElements(eltType); ++j) { // for each element
378
1/2
✓ Branch 2 taken 1670 times.
✗ Branch 3 not taken.
 1670 m_interfMesh->getOneElement(eltType, j, elemIdPoint, 0);
379
2/2
✓ Branch 0 taken 3340 times.
✓ Branch 1 taken 1670 times.
 5010 for(int iPoint = 0; iPoint < numPointsPerElt; iPoint++) {
380 3340 elemPoint[iPoint] = &m_interfMesh->listPoints()[elemIdPoint[iPoint]];
381 }
382
1/2
✓ Branch 1 taken 1670 times.
✗ Branch 2 not taken.
 1670 fe.updateMeas(0, elemPoint); // update the measure
383
2/2
✓ Branch 0 taken 3340 times.
✓ Branch 1 taken 1670 times.
 5010 for(int iPoint = 0; iPoint < numPointsPerElt; ++iPoint) {
384
1/2
✓ Branch 1 taken 3340 times.
✗ Branch 2 not taken.
 3340 m_lumpedMassMatrix[elemIdPoint[iPoint]] += fe.measure()/numPointsPerElt;
385 }
386 }
387 80 }
388
389
2/4
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 80 times.
 80 if ( FelisceParam::verbose() > 1 ) {
390 ✗ std::cout << " -- The lumped mass matrix of the structure for the immersed method has been successfully built!" << std::endl;
391 }
392 80 }
393
394 ✗ void Immersed::m_divDivStabPostProcess(std::vector<felInt>& elementID){
395
396 ✗ double dt = FelisceParam::instance().timeStep;
397
398 ✗ std::vector<GeometricMeshRegion::seedElement> elementIDstar;
399 ✗ std::vector<felInt> value(m_backMesh->getNumDomainElement(),0.);
400
401 ✗ for (std::size_t i=0; i < elementID.size(); ++i) {
402 ✗ value[elementID[i]] = 1;
403 ✗ m_backMesh->getElementPatch(m_backMesh->bagElementTypeDomain()[0], elementID[i],elementIDstar);
404 ✗ for (std::size_t j=0; j < elementIDstar.size(); ++j) {
405 ✗ value[elementIDstar[j].second] = 1;
406 }
407 }
408
409 {
410 // write case
411 FILE * pFile;
412 ✗ std::string fileName, fileNameDir, fileNameCase;
413
414 ✗ fileNameDir = FelisceParam::instance().resultDir;
415 ✗ fileNameCase = fileNameDir+"intersection.case";
416
417 ✗ pFile = fopen (fileNameCase.c_str(),"w");
418 ✗ if (pFile==NULL) {
419 ✗ std::string command = "mkdir -p " + fileNameDir;
420 ✗ int ierr = system(command.c_str());
421 ✗ if( ierr > 0)
422 ✗ FEL_ERROR("Impossible to write "+fileNameCase+" case file");
423 ✗ pFile = fopen (fileNameCase.c_str(),"w");
424 }
425 ✗ fprintf( pFile,"FORMAT\n");
426 // fprintf( pFile,"type: ensight\n");
427 ✗ fprintf( pFile,"type: ensight gold\n");
428 ✗ fprintf( pFile,"GEOMETRY\n");
429 // fprintf( pFile,"model: 1 fluid.geo.*****.geo\n");
430 ✗ fprintf( pFile,"model: 1 fluid.geo\n");
431 ✗ fprintf( pFile,"VARIABLE\n");
432 ✗ fprintf( pFile,"scalar per element: 1 intersection intersection.*****.scl\n");
433 ✗ fprintf( pFile,"TIME\n");
434 ✗ fprintf( pFile, "time set: %d\n", 1);
435 ✗ fprintf( pFile,"number of steps: %d \n", m_indexTime+1);
436 ✗ fprintf( pFile, "filename start number: %d\n", 0);
437 ✗ fprintf( pFile, "filename increment: %d\n", 1);
438 ✗ fprintf( pFile, "time values:\n");
439
440 ✗ int count=0;
441 ✗ for (int i=0; i < m_indexTime+1; ++i) {
442 ✗ fprintf( pFile, "%12.5e",i*dt);
443 ✗ ++count;
444 ✗ if ( count == 6) {
445 ✗ fprintf( pFile, "\n" );
446 ✗ count=0;
447 }
448 }
449 ✗ fclose(pFile);
450 }
451
452 {
453 ✗ std::string fileName = FelisceParam::instance().resultDir + "intersection";
454 ✗ std::stringstream oss;
455 ✗ oss << m_indexTime;
456 ✗ std::string aux = oss.str();
457 ✗ if ( m_indexTime < 10 ) {
458 ✗ aux = "0000" + aux;
459 ✗ } else if ( m_indexTime < 100 ) {
460 ✗ aux = "000" + aux;
461 ✗ } else if ( m_indexTime < 1000 ) {
462 ✗ aux = "00" + aux;
463 ✗ } else if ( m_indexTime < 10000 ) {
464 ✗ aux = "0" + aux;
465 }
466
467 ✗ fileName = fileName + "." + aux + ".scl";
468 FILE * pFile;
469 ✗ pFile = fopen (fileName.c_str(),"w");
470
471 ✗ if (pFile==nullptr) {
472 ✗ FEL_ERROR("Impossible to write " + fileName + " scalar ensight solution");
473 }
474
475 ✗ fprintf(pFile, "Scalar per element\n");
476 ✗ fprintf(pFile, "part\n");
477 ✗ fprintf(pFile, "%10d\n",1); // change the value at the end of the line depending on the part in your .geo file which will receive the "intersection" variable for the post-processing part
478 ✗ if (m_backMesh->domainDim() == 2) { // if the fluid mesh is a 2D mesh
479 ✗ fprintf(pFile, "tria3\n"); // consider the triangle elements
480 }
481 ✗ else if (m_backMesh->domainDim() == 3) { // if the fluid mesh is a 3D mesh
482 ✗ fprintf(pFile, "tetra4\n"); // consider the tetraedra elements
483 }
484
485 ✗ for (std::size_t i=0; i < value.size(); ++i) {
486 ✗ fprintf(pFile,"%d\n", value[i]);
487 }
488
489 ✗ fprintf(pFile, "\n");
490 ✗ fclose(pFile);
491 }
492 }
493 }
494
495