GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemNS.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 359 1903 18.9%
Branches: 302 2867 10.5%
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 & J-F. Gerbeau & V. Martin
13 //
14
15 // System includes
16 #include <stack>
17 #include <functional>
18
19 // External includes
20
21 // Project includes
22 #include "DegreeOfFreedom/boundaryCondition.hpp"
23 #include "Core/felisceTools.hpp"
24 #include "FiniteElement/elementMatrix.hpp"
25 #include "FiniteElement/elementVector.hpp"
26 #include "InputOutput/io.hpp"
27 #include "Solver/linearProblemNS.hpp"
28
29 /*!
30 \file LinearProblemNS.cpp
31
32 \date 05/01/2011
33 \brief Monolithic algorithm for the Navier-Stokes equations
34 */
35
36 namespace felisce {
37
38 54 LinearProblemNS::LinearProblemNS():
39 108 LinearProblem("Navier-Stokes monolithic",FelisceParam::instance().nonLinearFluid == true ? 3 : 2), // TODO FelisceParam::instance(???) which one??
40 54 m_initializeMapPnt(true),
41 54 m_viscosity(0.),
42 54 m_density(0.),
43 54 m_theta(0.),
44 54 m_lumpedModelBC(nullptr),
45 54 m_forceConvAndStabComputation(false),
46 54 m_bdf(nullptr),
47 54 m_cardiacCycle(nullptr),
48 54 m_ris(nullptr),
49 54 allocateSeqVec(false),
50 54 allocateSeqVecExt(false),
51
24/46
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 46 times.
✓ Branch 5 taken 54 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 54 times.
✗ Branch 9 not taken.
✓ Branch 20 taken 54 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 54 times.
✗ Branch 24 not taken.
✓ Branch 29 taken 54 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 54 times.
✗ Branch 33 not taken.
✓ Branch 35 taken 54 times.
✗ Branch 36 not taken.
✓ Branch 38 taken 54 times.
✗ Branch 39 not taken.
✓ Branch 41 taken 54 times.
✗ Branch 42 not taken.
✓ Branch 44 taken 54 times.
✗ Branch 45 not taken.
✓ Branch 47 taken 54 times.
✗ Branch 48 not taken.
✓ Branch 50 taken 54 times.
✗ Branch 51 not taken.
✓ Branch 53 taken 54 times.
✗ Branch 54 not taken.
✓ Branch 56 taken 54 times.
✗ Branch 57 not taken.
✓ Branch 60 taken 54 times.
✗ Branch 61 not taken.
✓ Branch 63 taken 54 times.
✗ Branch 64 not taken.
✓ Branch 69 taken 54 times.
✗ Branch 70 not taken.
✓ Branch 73 taken 54 times.
✗ Branch 74 not taken.
✓ Branch 76 taken 54 times.
✗ Branch 77 not taken.
✓ Branch 79 taken 54 times.
✗ Branch 80 not taken.
✓ Branch 84 taken 54 times.
✗ Branch 85 not taken.
✓ Branch 87 taken 54 times.
✗ Branch 88 not taken.
 54 m_isCsrNoInterfaceStored(false)
52 {
53 // When using cvgraph, the following vectors are needed.
54 // If we use a D2N or a N2D some of them may be unused
55
2/4
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 54 times.
 54 if ( FelisceParam::instance().withCVG ) { // TODO FelisceParam::instance(???) which one??
56 #ifdef FELISCE_WITH_CVGRAPH
57 ✗ m_cvgDirichletVariable = "cvgraphVelocity";
58 #endif
59 ✗ m_seqVecs.Init("cvgraphSTRESS");
60 ✗ m_seqVecs.Init("externalDisplacementOld");
61 ✗ m_seqVecs.Init("externalDisplacement");
62 ✗ m_seqVecs.Init("cvgraphVelocity");
63 }
64 54 }
65
66 /***********************************************************************************/
67 /***********************************************************************************/
68
69 108 LinearProblemNS::~LinearProblemNS()
70 {
71 108 m_elemFieldLumpedModelBC.clear();
72 108 m_bdf = nullptr;
73
1/2
✓ Branch 0 taken 54 times.
✗ Branch 1 not taken.
 108 if (allocateSeqVec) {
74 108 m_seqBdfRHS.destroy();
75 }
76
1/2
✓ Branch 0 taken 54 times.
✗ Branch 1 not taken.
 108 if (allocateSeqVecExt) {
77 108 m_solExtrapol.destroy();
78 }
79 }
80
81 /***********************************************************************************/
82 /***********************************************************************************/
83
84 1436 PetscErrorCode NSJacobian(SNES snes, Vec evaluationState, Mat jacobian , Mat preconditioner, void* context_as_void)
85 {
86 (void) snes;
87 (void) evaluationState;
88 (void) jacobian;
89 (void) preconditioner;
90 (void) context_as_void;
91 PetscFunctionBegin;
92
93 1436 SnesContext* context_ptr = static_cast<SnesContext*>(context_as_void);
94 1436 SnesContext& context = *context_ptr;
95 1436 LinearProblemNS& lpb = *static_cast<LinearProblemNS*>(context.linearProblem);
96
97 // Retrieve instance
98 1436 auto& r_instance = FelisceParam::instance(lpb.instanceIndex());
99
100 // add matrices linear and non linear terms
101 1436 lpb.addMatrixRHSNl(FlagMatrixRHS::only_matrix);
102
103 // remove rows of nl matrices that have prescribed velocities (Dirichlet)
104 2872 lpb.applyBC(r_instance.essentialBoundaryConditionsMethod,
105 1436 MpiInfo::rankProc(), FlagMatrixRHS::only_matrix, FlagMatrixRHS::only_matrix,
106 0, true, ApplyNaturalBoundaryConditions::no);
107
108 1436 PetscFunctionReturn(0);
109 }
110
111 /***********************************************************************************/
112 /***********************************************************************************/
113
114 2396 PetscErrorCode NSFunction(SNES snes, Vec evaluationState, Vec residual, void* context_as_void)
115 {
116 PetscFunctionBegin;
117 2396 SnesContext* context_ptr = static_cast<SnesContext*>(context_as_void);
118 2396 SnesContext& context = *context_ptr;
119 2396 LinearProblemNS& lpb = *static_cast<LinearProblemNS*>(context.linearProblem);
120
121 // Retrieve instance
122
2/4
✓ Branch 1 taken 2396 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2396 times.
✗ Branch 5 not taken.
 2396 auto& r_instance = FelisceParam::instance(lpb.instanceIndex());
123
124 PetscInt nbEval;
125
1/2
✓ Branch 1 taken 2396 times.
✗ Branch 2 not taken.
 2396 SNESGetNumberFunctionEvals(snes,&nbEval);
126
127 // transfer evaluation state U^(n, k-1) to computeElementArray
128
2/4
✓ Branch 2 taken 2396 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2396 times.
✗ Branch 6 not taken.
 2396 VecCopy(evaluationState,lpb.evaluationState().toPetsc());
129
1/2
✓ Branch 3 taken 2396 times.
✗ Branch 4 not taken.
 2396 lpb.gatherVector(lpb.evaluationState(), lpb.seqEvaluationState());
130
131 // add solution at previous time step
132
4/4
✓ Branch 0 taken 960 times.
✓ Branch 1 taken 1436 times.
✓ Branch 2 taken 160 times.
✓ Branch 3 taken 800 times.
 2396 if (nbEval == 0 && r_instance.useALEformulation)
133
1/2
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
 160 lpb.addMatrixRHSNl(FlagMatrixRHS::only_rhs);
134
135 // clear non linear matrix
136
1/2
✓ Branch 2 taken 2396 times.
✗ Branch 3 not taken.
 2396 lpb.matrix(2).zeroEntries();
137
138
2/2
✓ Branch 0 taken 1436 times.
✓ Branch 1 taken 960 times.
 2396 if (nbEval > 0 ) {
139 // assembly after first inner iteration
140
1/2
✓ Branch 1 taken 1436 times.
✗ Branch 2 not taken.
 1436 lpb.clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs);
141
1/2
✓ Branch 1 taken 1436 times.
✗ Branch 2 not taken.
 1436 lpb.assembleMatrixRHS(context.rankProc);
142
143 // add static and dynamic matrices
144
1/2
✓ Branch 1 taken 1436 times.
✗ Branch 2 not taken.
 1436 lpb.addMatrixRHS();
145
146 // add solution at previous time step
147
2/2
✓ Branch 0 taken 316 times.
✓ Branch 1 taken 1120 times.
 1436 if (r_instance.useALEformulation)
148
1/2
✓ Branch 1 taken 316 times.
✗ Branch 2 not taken.
 316 lpb.addMatrixRHSNl(FlagMatrixRHS::only_rhs);
149
150 // apply boundary conditions.
151
1/2
✓ Branch 1 taken 1436 times.
✗ Branch 2 not taken.
 1436 lpb.finalizeEssBCTransient();
152
1/2
✓ Branch 1 taken 1436 times.
✗ Branch 2 not taken.
 1436 lpb.applyBC(r_instance.essentialBoundaryConditionsMethod,context.rankProc);
153 }
154
155 // compute residual
156
1/2
✓ Branch 1 taken 2396 times.
✗ Branch 2 not taken.
 2396 lpb.evalDynamicResidual();
157
158 // send residual to petsc
159
2/4
✓ Branch 2 taken 2396 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2396 times.
✗ Branch 6 not taken.
 2396 VecCopy(lpb.vector().toPetsc(),residual);
160
161 2396 PetscFunctionReturn(0);
162
163 }
164
165 2396 void LinearProblemNS::evalDynamicResidual()
166 {
167 // compute residual
168 2396 multAdd(matrix(0),evaluationState(),vector(),vector());
169 2396 }
170
171 54 void LinearProblemNS::InitializeDerivedAttributes()
172 {
173 // Retrieve instance
174 54 auto& r_instance = FelisceParam::instance(instanceIndex());
175
176
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 46 times.
 54 if (r_instance.nonLinearFluid) {
177 8 m_snesContext.initialize(*this, MpiInfo::rankProc(), MpiInfo::numProc(), ApplyNaturalBoundaryConditions::no, FlagMatrixRHS::matrix_and_rhs);
178 8 snesInterface().setFunction(vector(), NSFunction, &m_snesContext);
179 8 snesInterface().setJacobian(matrix(0), matrix(0), NSJacobian, &m_snesContext);
180 }
181 54 }
182
183 /***********************************************************************************/
184 /***********************************************************************************/
185
186 54 void LinearProblemNS::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES)
187 {
188 // Retrieve instance
189
2/4
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
 54 auto& r_instance = FelisceParam::instance(instanceIndex());
190
191
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 LinearProblem::initialize(mesh, comm, doUseSNES);
192 54 m_fstransient = fstransient;
193 54 std::vector<PhysicalVariable> listVariable;
194 54 std::vector<std::size_t> listNumComp;
195
196
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 listVariable.push_back(velocity);
197
1/2
✓ Branch 2 taken 54 times.
✗ Branch 3 not taken.
 54 listNumComp.push_back(dimension());
198
199
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 listVariable.push_back(pressure);
200
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 listNumComp.push_back(1);
201
202
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 36 times.
 54 if ( r_instance.useALEformulation ) {
203
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 listVariable.push_back(displacement);
204
1/2
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
 18 listNumComp.push_back(dimension());
205 }
206
207 // define unknown of the linear system.
208
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 m_listUnknown.push_back(velocity);
209
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 m_listUnknown.push_back(pressure);
210
211
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
 54 if ( r_instance.useFDlagrangeMult ) {
212 ✗ listVariable.push_back(lagMultiplier);
213 ✗ listNumComp.push_back(dimension());
214 ✗ m_listUnknown.push_back(lagMultiplier);
215
216 ✗ if ( r_instance.massConstraint ) {
217 ✗ listVariable.push_back(temperature);
218 ✗ listNumComp.push_back(1);
219 ✗ m_listUnknown.push_back(temperature);
220 }
221 }
222
223
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 userAddOtherUnknowns(listVariable,listNumComp);
224
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 userAddOtherVariables(listVariable,listNumComp);
225
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 definePhysicalVariable(listVariable,listNumComp);
226
227
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 m_iVelocity = m_listVariable.getVariableIdList(velocity);
228
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 m_iPressure = m_listVariable.getVariableIdList(pressure);
229
230 54 m_velocity = &m_listVariable[m_iVelocity];
231 54 m_pressure = &m_listVariable[m_iPressure];
232
233
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 m_iUnknownVel = m_listUnknown.getUnknownIdList(velocity);
234
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 m_iUnknownPre = m_listUnknown.getUnknownIdList(pressure);
235
236
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 m_velBlock = m_listUnknown.getBlockPosition(m_iUnknownVel,0);
237
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 m_preBlock = m_listUnknown.getBlockPosition(m_iUnknownPre,0);
238
239
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
 54 if (r_instance.useFDlagrangeMult) {
240
241 ✗ if ( r_instance.massConstraint ) {
242 ✗ if (r_instance.BHstab) {
243
244 ✗ if ( dimension() == 2 ) {
245 std::vector<int> mask{ 1, 1, /**/ 1, /**/ 1, 1, /**/ 1,
246 1, 1, /**/ 1, /**/ 1, 1, /**/ 1,
247 1, 1, /**/ 1, /**/ 0, 0, /**/ 0,
248 1, 1, /**/ 0, /**/ 1, 1, /**/ 1,
249 1, 1, /**/ 0, /**/ 1, 1, /**/ 1,
250 ✗ 1, 1, /**/ 0, /**/ 1, 1, /**/ 1};
251 ✗ m_listUnknown.setMask(mask);
252 }
253 ✗ else if ( dimension() == 3 ) {
254 std::vector<int> mask{ 1, 1, 1, /**/ 1, /**/ 1, 1, 1, /**/ 1,
255 1, 1, 1, /**/ 1, /**/ 1, 1, 1, /**/ 1,
256 1, 1, 1, /**/ 1, /**/ 1, 1, 1, /**/ 1,
257 1, 1, 1, /**/ 1, /**/ 0, 0, 0, /**/ 0,
258 1, 1, 1, /**/ 0, /**/ 1, 1, 1, /**/ 1,
259 1, 1, 1, /**/ 0, /**/ 1, 1, 1, /**/ 1,
260 1, 1, 1, /**/ 0, /**/ 1, 1, 1, /**/ 1,
261 ✗ 1, 1, 1, /**/ 0, /**/ 1, 1, 1, /**/ 1};
262 ✗ m_listUnknown.setMask(mask);
263 }
264 }
265 ✗ else if (r_instance.BPstab) {
266
267 ✗ if ( dimension() == 2 ) {
268 std::vector<int> mask{ 1, 1, /**/ 1, /**/ 1, 1, /**/ 1,
269 1, 1, /**/ 1, /**/ 1, 1, /**/ 1,
270 1, 1, /**/ 1, /**/ 0, 0, /**/ 0,
271 1, 1, /**/ 0, /**/ 1, 1, /**/ 0,
272 1, 1, /**/ 0, /**/ 1, 1, /**/ 0,
273 ✗ 1, 1, /**/ 0, /**/ 0, 0, /**/ 1};
274 ✗ m_listUnknown.setMask(mask);
275 }
276 ✗ else if ( dimension() == 3 ) {
277 std::vector<int> mask{ 1, 1, 1, /**/ 1, /**/ 1, 1, 1, /**/ 1,
278 1, 1, 1, /**/ 1, /**/ 1, 1, 1, /**/ 1,
279 1, 1, 1, /**/ 1, /**/ 1, 1, 1, /**/ 1,
280 1, 1, 1, /**/ 1, /**/ 0, 0, 0, /**/ 0,
281 1, 1, 1, /**/ 0, /**/ 1, 1, 1, /**/ 0,
282 1, 1, 1, /**/ 0, /**/ 1, 1, 1, /**/ 0,
283 1, 1, 1, /**/ 0, /**/ 1, 1, 1, /**/ 0,
284 ✗ 1, 1, 1, /**/ 0, /**/ 0, 0, 0, /**/ 1};
285 ✗ m_listUnknown.setMask(mask);
286 }
287 }
288 else {
289
290 ✗ if ( dimension() == 2 ) {
291 std::vector<int> mask{ 1, 1, /**/ 1, /**/ 1, 1, /**/ 1,
292 1, 1, /**/ 1, /**/ 1, 1, /**/ 1,
293 1, 1, /**/ 1, /**/ 0, 0, /**/ 0,
294 1, 1, /**/ 0, /**/ 0, 0, /**/ 0,
295 1, 1, /**/ 0, /**/ 0, 0, /**/ 0,
296 ✗ 1, 1, /**/ 0, /**/ 0, 0, /**/ 0};
297 ✗ m_listUnknown.setMask(mask);
298 }
299 ✗ else if ( dimension() == 3 ) {
300 std::vector<int> mask{ 1, 1, 1, /**/ 1, /**/ 1, 1, 1, /**/ 1,
301 1, 1, 1, /**/ 1, /**/ 1, 1, 1, /**/ 1,
302 1, 1, 1, /**/ 1, /**/ 1, 1, 1, /**/ 1,
303 1, 1, 1, /**/ 1, /**/ 0, 0, 0, /**/ 0,
304 1, 1, 1, /**/ 0, /**/ 0, 0, 0, /**/ 0,
305 1, 1, 1, /**/ 0, /**/ 0, 0, 0, /**/ 0,
306 1, 1, 1, /**/ 0, /**/ 0, 0, 0, /**/ 0,
307 ✗ 1, 1, 1, /**/ 0, /**/ 0, 0, 0, /**/ 0};
308 ✗ m_listUnknown.setMask(mask);
309 }
310 }
311 }
312 else {
313 ✗ if (r_instance.BHstab) {
314
315 ✗ if ( dimension() == 2 ) {
316 std::vector<int> mask{ 1, 1, /**/ 1, /**/ 1, 1,
317 1, 1, /**/ 1, /**/ 1, 1,
318 1, 1, /**/ 1, /**/ 0, 0,
319 1, 1, /**/ 0, /**/ 1, 1,
320 ✗ 1, 1, /**/ 0, /**/ 1, 1};
321 ✗ m_listUnknown.setMask(mask);
322 }
323 ✗ else if ( dimension() == 3 ) {
324 std::vector<int> mask{ 1, 1, 1, /**/ 1, /**/ 1, 1, 1,
325 1, 1, 1, /**/ 1, /**/ 1, 1, 1,
326 1, 1, 1, /**/ 1, /**/ 1, 1, 1,
327 1, 1, 1, /**/ 1, /**/ 0, 0, 0,
328 1, 1, 1, /**/ 0, /**/ 1, 1, 1,
329 1, 1, 1, /**/ 0, /**/ 1, 1, 1,
330 ✗ 1, 1, 1, /**/ 0, /**/ 1, 1, 1};
331 ✗ m_listUnknown.setMask(mask);
332 }
333 }
334 ✗ else if (r_instance.BPstab) {
335
336 ✗ if ( dimension() == 2 ) {
337 std::vector<int> mask{ 1, 1, /**/ 1, /**/ 1, 1,
338 1, 1, /**/ 1, /**/ 1, 1,
339 1, 1, /**/ 1, /**/ 0, 0,
340 1, 1, /**/ 0, /**/ 1, 1,
341 ✗ 1, 1, /**/ 0, /**/ 1, 1};
342 ✗ m_listUnknown.setMask(mask);
343 }
344 ✗ else if ( dimension() == 3 ) {
345 std::vector<int> mask{ 1, 1, 1, /**/ 1, /**/ 1, 1, 1,
346 1, 1, 1, /**/ 1, /**/ 1, 1, 1,
347 1, 1, 1, /**/ 1, /**/ 1, 1, 1,
348 1, 1, 1, /**/ 1, /**/ 0, 0, 0,
349 1, 1, 1, /**/ 0, /**/ 1, 1, 1,
350 1, 1, 1, /**/ 0, /**/ 1, 1, 1,
351 ✗ 1, 1, 1, /**/ 0, /**/ 1, 1, 1};
352 ✗ m_listUnknown.setMask(mask);
353 }
354 }
355 else {
356
357 ✗ if ( dimension() == 2 ) {
358 std::vector<int> mask{ 1, 1, /**/ 1, /**/ 1, 1,
359 1, 1, /**/ 1, /**/ 1, 1,
360 1, 1, /**/ 1, /**/ 0, 0,
361 1, 1, /**/ 0, /**/ 0, 0,
362 ✗ 1, 1, /**/ 0, /**/ 0, 0};
363 ✗ m_listUnknown.setMask(mask);
364 }
365 ✗ else if ( dimension() == 3 ) {
366 std::vector<int> mask{ 1, 1, 1, /**/ 1, /**/ 1, 1, 1,
367 1, 1, 1, /**/ 1, /**/ 1, 1, 1,
368 1, 1, 1, /**/ 1, /**/ 1, 1, 1,
369 1, 1, 1, /**/ 1, /**/ 0, 0, 0,
370 1, 1, 1, /**/ 0, /**/ 0, 0, 0,
371 1, 1, 1, /**/ 0, /**/ 0, 0, 0,
372 ✗ 1, 1, 1, /**/ 0, /**/ 0, 0, 0};
373 ✗ m_listUnknown.setMask(mask);
374 }
375 }
376 }
377
378 ✗ m_iLagrange = m_listVariable.getVariableIdList(lagMultiplier);
379 ✗ m_iUnknownLag = m_listUnknown.getUnknownIdList(lagMultiplier);
380 ✗ m_lagBlock = m_listUnknown.getBlockPosition(m_iUnknownLag,0);
381
382 ✗ if ( r_instance.massConstraint ) {
383 ✗ m_iJapanCon = m_listVariable.getVariableIdList(temperature);
384 ✗ m_iUnknownJap = m_listUnknown.getUnknownIdList(temperature);
385 ✗ m_japBlock = m_listUnknown.getBlockPosition(m_iUnknownJap,0);
386 }
387
388 ✗ m_listUnknown.setUnknownsRows(std::vector<int>{m_iUnknownVel,m_iUnknownPre});
389 ✗ m_listUnknown.setUnknownsCols(std::vector<int>{m_iUnknownVel,m_iUnknownPre});
390 }
391
392 54 m_viscosity = r_instance.viscosity;
393 54 m_density = r_instance.density;
394 54 m_theta = r_instance.theta;
395 54 m_useSymmetricStress = r_instance.useSymmetricStress;
396
397
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 const int verbose = FelisceParam::verbose();
398 54 const double tolrescaling = r_instance.Geo.tolrescaling;
399
400
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
 54 if ( r_instance.characteristicMethod )
401 ✗ m_pLocator = felisce::make_shared<Locator>(mesh[m_currentMesh].get(), false, tolrescaling, verbose);
402
403
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
 54 if ( r_instance.useFDlagrangeMult )
404 ✗ m_pLocator = felisce::make_shared<Locator>(mesh[m_currentMesh].get(), true, tolrescaling, verbose);
405 54 }
406
407 /***********************************************************************************/
408 /***********************************************************************************/
409
410 3940 void LinearProblemNS::initPerElementType(ElementType /* eltType */, FlagMatrixRHS /* flagMatrixRHS */)
411 {
412 // Retrieve instance
413 3940 auto& r_instance = FelisceParam::instance(instanceIndex());
414
415
1/2
✓ Branch 2 taken 3940 times.
✗ Branch 3 not taken.
 3940 if ( m_currentMesh == m_listVariable[m_iVelocity].idMesh() ) {
416 3940 m_feVel = m_listCurrentFiniteElement[m_iVelocity];
417 3940 m_fePres = m_listCurrentFiniteElement[m_iPressure];
418
419 3940 m_elemFieldAdv.initialize(DOF_FIELD,*m_feVel,dimension());
420 3940 m_elemFieldRHS.initialize(DOF_FIELD,*m_feVel,dimension());
421 3940 m_elemFieldRHSbdf.initialize(DOF_FIELD,*m_feVel,dimension());
422 3940 m_elemFieldTotalPressureFormulation.initialize(QUAD_POINT_FIELD,*m_feVel,dimension());
423 3940 m_elemFieldCharact.initialize(QUAD_POINT_FIELD,*m_feVel,dimension());
424
425
2/2
✓ Branch 0 taken 2720 times.
✓ Branch 1 taken 1220 times.
 3940 if (r_instance.nonLinearFluid) {
426 2720 m_elemFieldVelSeq.initialize(DOF_FIELD,*m_feVel,dimension());
427 2720 m_elemFieldPresSeq.initialize(DOF_FIELD, *m_fePres, 1);
428 }
429
430 // theta-method RHS
431
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3940 times.
 3940 if ( m_theta < 1 ) {
432 ✗ m_elemFieldPres.initialize(DOF_FIELD,*m_fePres,1);
433 ✗ std::size_t SizeElMatRHSThetaMethod = 2;
434 ✗ if(r_instance.useALEformulation)
435 ✗ SizeElMatRHSThetaMethod = 1;
436
437 ✗ const std::vector<const CurBaseFiniteElement*> finiteElt({m_feVel, m_fePres});
438 ✗ std::vector<std::size_t> numberCmp({m_listVariable[m_iVelocity].numComponent(),m_listVariable[m_iPressure].numComponent()});
439 ✗ for (std::size_t i = 0; i < SizeElMatRHSThetaMethod; i++)
440 ✗ m_elMatRHSThetaMethod.push_back(new ElementMatrix(finiteElt,numberCmp, numberCmp));
441 }
442
443 //=========
444 // ALE
445 //=========
446
2/2
✓ Branch 0 taken 1396 times.
✓ Branch 1 taken 2544 times.
 3940 if(r_instance.useALEformulation) {
447
448 1396 m_iDisplacement = m_listVariable.getVariableIdList(displacement);
449 1396 m_displacement = &m_listVariable[m_iDisplacement];
450
451 1396 m_elemFieldVelMesh.initialize(DOF_FIELD,*m_feVel,dimension());
452
453 1396 numDofExtensionProblem = m_externalDof[0]->numDof();
454
455 1396 m_petscToGlobal1.resize(numDofExtensionProblem);
456 1396 m_petscToGlobal2.resize(numDofExtensionProblem);
457 1396 m_auxvec.resize(numDofExtensionProblem);
458
459
2/2
✓ Branch 0 taken 24291396 times.
✓ Branch 1 taken 1396 times.
 24292792 for (int i=0; i<numDofExtensionProblem; i++) {
460 24291396 m_petscToGlobal1[i]=i;
461 24291396 m_petscToGlobal2[i]=i;
462 }
463 1396 AOApplicationToPetsc(m_externalAO[0],numDofExtensionProblem,m_petscToGlobal1.data());
464 1396 AOApplicationToPetsc(m_ao,numDofExtensionProblem,m_petscToGlobal2.data());
465
466
1/2
✓ Branch 1 taken 1396 times.
✗ Branch 2 not taken.
 1396 if ( m_fstransient->iteration > 0) {
467
2/2
✓ Branch 1 taken 58 times.
✓ Branch 2 taken 1338 times.
 1396 if ( m_fstransient->iteration == 1) {
468 58 m_beta.duplicateFrom(m_solExtrapol);
469
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 42 times.
 58 if (r_instance.nonLinearFluid)
470 16 m_betaSeq.duplicateFrom(seqEvaluationState());
471 }
472 //=========
473 // ALE
474 //=========
475 // Construction of m_beta = u^{n-1} - w^{n}
476 // Recall that externalVec(0) = - w^{n}
477 1396 m_beta.copyFrom(m_solExtrapol);
478 1396 externalVec(0).getValues(numDofExtensionProblem,m_petscToGlobal1.data(), m_auxvec.data());
479 1396 m_beta.setValues(numDofExtensionProblem, m_petscToGlobal2.data() ,m_auxvec.data(), ADD_VALUES);
480
2/2
✓ Branch 0 taken 796 times.
✓ Branch 1 taken 600 times.
 1396 if (r_instance.nonLinearFluid){
481 796 m_auxvec.resize(numDofExtensionProblem);
482 796 externalVec(0).getValues(numDofExtensionProblem,m_petscToGlobal1.data(), m_auxvec.data());
483 796 m_betaSeq.copyFrom(seqEvaluationState());
484 796 m_betaSeq.setValues(numDofExtensionProblem, m_petscToGlobal2.data() ,m_auxvec.data(), ADD_VALUES);
485 }
486 }
487 }
488 }
489 3940 }
490
491 /***********************************************************************************/
492 /***********************************************************************************/
493
494 ✗ void LinearProblemNS::initializeRISModel(RISModel* risModel)
495 {
496 ✗ m_ris = risModel;
497 }
498
499 /***********************************************************************************/
500 /***********************************************************************************/
501
502 ✗ void LinearProblemNS::initializeCardiacCycle(CardiacCycle* cardiacCycle)
503 {
504 ✗ m_cardiacCycle = cardiacCycle;
505 }
506
507 /***********************************************************************************/
508 /***********************************************************************************/
509
510 8 void LinearProblemNS::initializeLumpedModelBC(LumpedModelBC* lumpedModelBC)
511 {
512 8 m_lumpedModelBC = lumpedModelBC;
513 8 m_elemFieldLumpedModelBC.clear();
514 8 m_elemFieldLumpedModelBC.resize(m_lumpedModelBC->size());
515 8 }
516
517 /***********************************************************************************/
518 /***********************************************************************************/
519
520 2508 void LinearProblemNS::initPerElementTypeBoundaryCondition(ElementType& /* eltType */, FlagMatrixRHS /* flagMatrixRHS */)
521 {
522 // Retrieve instance
523 2508 auto& r_instance = FelisceParam::instance(instanceIndex());
524
525
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2508 times.
 2508 if (r_instance.CardiacCycle) {
526 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
527 ✗ m_elemFieldSistole.initialize(CONSTANT_FIELD, *m_curvFeVel, dimension());
528 ✗ m_elemFieldDiastole.initialize(CONSTANT_FIELD, *m_curvFeVel);
529 }
530
531
2/2
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 2424 times.
 2508 if (m_lumpedModelBC) {
532 // allocate m_elemFieldLumpedModelBC
533
2/2
✓ Branch 1 taken 168 times.
✓ Branch 2 taken 84 times.
 252 for (std::size_t iLumpedModelBC=0; iLumpedModelBC <m_lumpedModelBC->size() ; iLumpedModelBC++) {
534
2/8
✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 168 times.
✗ Branch 10 not taken.
 168 if(r_instance.lumpedModelBCAlgo[iLumpedModelBC] == 1 || (r_instance.lumpedModelBCAlgo[iLumpedModelBC] == 2 && r_instance.model == "NS"))
535 168 m_elemFieldLumpedModelBC[iLumpedModelBC].initialize(CONSTANT_FIELD,1);
536 }
537 }
538
539 #ifdef FELISCE_WITH_CVGRAPH
540
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2508 times.
 2508 if (r_instance.withCVG ) {
541 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
542 ✗ if ( slave()->thereIsAtLeastOneRobinCondition() ){
543 ✗ m_robinAux.initialize( DOF_FIELD, *m_listCurvilinearFiniteElement[m_iVelocity], dimension() );
544 }
545 }
546 #endif
547 2508 }
548
549 /***********************************************************************************/
550 /***********************************************************************************/
551
552 7330 void LinearProblemNS::initPerDomainBoundaryCondition(std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint,
553 int label, felInt numEltPerLabel,felInt* ielSupportDof,
554 ElementType& eltType, felInt numElement_eltType, FlagMatrixRHS flagMatrixRHS)
555 {
556 // Retrieve instance
557 7330 auto& r_instance = FelisceParam::instance(instanceIndex());
558
559 7330 m_currentLabel = label;
560
561 //================================
562 // Explicit/implicit lumpedModel boundary conditions with model NS
563 //================================
564
5/6
✓ Branch 0 taken 312 times.
✓ Branch 1 taken 7018 times.
✓ Branch 3 taken 312 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 312 times.
✓ Branch 6 taken 7018 times.
 7330 if (m_lumpedModelBC && (r_instance.model == "NS") ) {
565
2/2
✓ Branch 1 taken 624 times.
✓ Branch 2 taken 312 times.
 936 for (std::size_t iLumpedModelBC=0; iLumpedModelBC <m_lumpedModelBC->size() ; iLumpedModelBC++) {
566
2/6
✗ Branch 1 not taken.
✓ Branch 2 taken 624 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 624 times.
✗ Branch 7 not taken.
 624 if(r_instance.lumpedModelBCAlgo[iLumpedModelBC] == 1 || r_instance.lumpedModelBCAlgo[iLumpedModelBC] == 2) {
567
2/2
✓ Branch 1 taken 42 times.
✓ Branch 2 taken 582 times.
 624 if (r_instance.lumpedModelBCLabel[iLumpedModelBC] == label) {
568 42 applyNaturalBoundaryConditionLumpedModelBC(elemPoint,elemIdPoint,label,numEltPerLabel, ielSupportDof,iLumpedModelBC,eltType, numElement_eltType, flagMatrixRHS);
569 }
570 }
571 }
572 }
573 7330 }
574
575 /***********************************************************************************/
576 /***********************************************************************************/
577
578 338616 void LinearProblemNS::computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /* elemIdPoint */, felInt& iel, FlagMatrixRHS /* flagMatrixRHS */)
579 {
580 // Retrieve instance
581 338616 auto& r_instance = FelisceParam::instance(instanceIndex());
582
583 // General updates
584
2/2
✓ Branch 1 taken 744224 times.
✓ Branch 2 taken 338616 times.
 1082840 for (std::size_t iFe = 0; iFe < m_listCurvilinearFiniteElement.size(); iFe++) {
585 744224 m_listCurvilinearFiniteElement[iFe]->updateMeasNormal(0, elemPoint);
586 }
587
588 // Saving ptr to velocity curvilinear finite element
589 338616 CurvilinearFiniteElement* curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
590 // Position of the velocity block
591 338616 int iblockVel = m_listUnknown.getBlockPosition(m_iUnknownVel,0);
592
593 // Windkessels model!
594
5/6
✓ Branch 0 taken 195812 times.
✓ Branch 1 taken 142804 times.
✓ Branch 3 taken 195812 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 195812 times.
✓ Branch 6 taken 142804 times.
 338616 if (m_lumpedModelBC && r_instance.model == "NS") {
595
2/2
✓ Branch 1 taken 391624 times.
✓ Branch 2 taken 195812 times.
 587436 for (std::size_t iLumpedModelBC=0; iLumpedModelBC <m_lumpedModelBC->size() ; iLumpedModelBC++) {
596
2/2
✓ Branch 1 taken 1120 times.
✓ Branch 2 taken 390504 times.
 391624 if ( r_instance.lumpedModelBCLabel[iLumpedModelBC] == m_currentLabel) {
597
2/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1120 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1120 times.
✗ Branch 7 not taken.
 1120 if(r_instance.lumpedModelBCAlgo[iLumpedModelBC] == 1 || r_instance.lumpedModelBCAlgo[iLumpedModelBC] == 2) {
598
599 1120 m_elementVectorBD[0]->f_phi_i_scalar_n(1.,*curvFeVel,m_elemFieldLumpedModelBC[iLumpedModelBC],iblockVel);
600 }
601 // to be continued // Justine 7/11/13
602 //if ( r_instance.lumpedModelBCAlgoImprovedExplicit == 1 ) {
603 // we add \int_Gamma (u^{n+1} \cdot n) (v \cdot n)
604 // and \int_Gamma (u^n \cdot n) (v \cdot n)
605 //}
606 }
607 }
608 }
609
610
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 338616 times.
 338616 if(r_instance.CardiacCycleLabel == m_currentLabel) {
611 ✗ double penalisationParam = 0.;
612 ✗ penalisationParam = m_cardiacCycle->penalisationParam(); //= 0 if valve closed, = TGV if valve open
613 ✗ m_elemFieldSistole.setValue(m_cardiacCycle->inflow());
614 ✗ m_elemFieldDiastole.setValue(m_cardiacCycle->Pressure_in());
615
616 ✗ m_elementMatBD[0]->phi_i_phi_j(penalisationParam, *curvFeVel, iblockVel, iblockVel, dimension()); // epsilon * int_{Gamma_in} [ u^{n+1} * v ]
617 ✗ m_elementVectorBD[0]->source(penalisationParam, *curvFeVel, m_elemFieldSistole, iblockVel, dimension()); // epsilon * int_{Gamma_in} [ u_{in} * v ]
618 ✗ m_elementVectorBD[0]->f_phi_i_scalar_n(-1., *curvFeVel, m_elemFieldDiastole, iblockVel); // int_{Gamma_in} [-p * N* v ] curvFeVel referred to test function v
619 }
620
621
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 338616 times.
 338616 if (r_instance.RISModels) {
622 ✗ if (r_instance.PenalizationMethod)
623 ✗ FEL_ERROR("RIS and penalization models are activated at the same time!");
624 ✗ double Res = 0.;
625 ✗ Res = m_ris->resistance(m_currentLabel);
626 ✗ if( Tools::notEqual(Res,0.) ) {
627 //LHS term of the penalisation
628 ✗ if ( std::find(r_instance.forbiddenLabelsForCorrection.begin(), r_instance.forbiddenLabelsForCorrection.end(), m_currentLabel) == r_instance.forbiddenLabelsForCorrection.end() ) {
629 ✗ m_elementMatBD[0]->phi_i_phi_j(Res,*curvFeVel, iblockVel, iblockVel, dimension());
630 }
631 //RHS term of the penalisation
632 //Correct the constraint with the velocity of the mesh (R*[u-w] = 0)
633 //This term is not added at the first iteration as the first iteration is not physical in our case (A.This, M. Fernandez, J-F. Gerbeau : Zygot case with full deformation).
634 ✗ if(m_fstransient->iteration>1){
635 ✗ if ( std::find(r_instance.forbiddenLabelsForCorrection.begin(), r_instance.forbiddenLabelsForCorrection.end(), m_currentLabel) == r_instance.forbiddenLabelsForCorrection.end() && r_instance.useALEformulation ) {
636 ✗ m_elementVectorBD[0]->source(-Res, *curvFeVel, m_elemFieldVelMesh, iblockVel, dimension());
637 }
638 }
639 }
640 }
641
642
2/2
✓ Branch 0 taken 46080 times.
✓ Branch 1 taken 292536 times.
 338616 if (r_instance.PenalizationMethod) {
643
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 46080 times.
 46080 if (r_instance.RISModels)
644 ✗ FEL_ERROR("RIS and penalization models are activated at the same time!");
645 46080 double Pen = r_instance.Gamma_penalization;
646
5/6
✓ Branch 1 taken 42240 times.
✓ Branch 2 taken 3840 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 42240 times.
✓ Branch 6 taken 3840 times.
✓ Branch 7 taken 42240 times.
 46080 if ( r_instance.closed_penalization[0] == m_currentLabel || r_instance.closed_penalization[1] == m_currentLabel ){
647
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3840 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 3840 times.
 3840 if (r_instance.useALEformulation && m_fstransient->iteration>1 )
648 ✗ m_elemFieldVelSeq.setValue(m_betaSeq, *m_feVel, iel, m_iVelocity, m_ao, dof());
649 else
650 3840 m_elemFieldVelSeq.setValue(seqEvaluationState(), *m_feVel, iel, m_iVelocity, m_ao, dof());
651 //LHS term of penalization
652 3840 m_elementMatBD[2]->penalty_phi_i_scalar_n_phi_j_scalar_n(Pen, *curvFeVel, m_elemFieldVelSeq, iblockVel, iblockVel);
653 //RHS term of penalization
654 3840 m_elementVectorBD[0]->penalty_phi_i_scalar_n(-Pen, *curvFeVel, m_elemFieldVelSeq, iblockVel);
655 }
656 }
657
658 #ifdef FELISCE_WITH_CVGRAPH
659
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 338616 times.
 338616 if (r_instance.withCVG )
660 ✗ cvgraphNaturalBC(iel);
661 #endif
662 338616 }
663
664 /***********************************************************************************/
665 /***********************************************************************************/
666
667 13890684 void LinearProblemNS::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /* elemIdPoint */, felInt& iel, FlagMatrixRHS flagMatrixRHS)
668 {
669 // Retrieve instance
670 13890684 auto& r_instance = FelisceParam::instance(instanceIndex());
671
672
1/2
✓ Branch 2 taken 13890684 times.
✗ Branch 3 not taken.
 13890684 if ( m_currentMesh == m_listVariable[m_iVelocity].idMesh() ) {
673 13890684 m_feVel->updateFirstDeriv(0, elemPoint);
674 13890684 m_fePres->updateFirstDeriv(0, elemPoint);
675
676 13890684 PetscScalar coef = m_density/m_fstransient->timeStep;
677
678
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 13890684 times.
 13890684 if ( m_feVel->measure() < 0. ) {
679 int iglo;
680 ✗ ISLocalToGlobalMappingApply(m_mappingElem[m_currentMesh],1,&iel,&iglo);
681 ✗ std::cout << " id global = " << iglo << std::endl;
682 ✗ FEL_ERROR("Error: Element of negative measure!");
683 }
684
685 std::size_t indexMatrix;
686
2/2
✓ Branch 0 taken 12276990 times.
✓ Branch 1 taken 1613694 times.
 13890684 if( r_instance.useALEformulation )
687 12276990 indexMatrix = 0; // The matrix actually used in the solver, which is not constant (because of advection, ALE,...)
688 else
689 1613694 indexMatrix = 1; // The "static" part of the matrix, when there is no ALE it is basically everything, but the advection.
690
691
5/6
✓ Branch 1 taken 134158 times.
✓ Branch 2 taken 13756526 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 134158 times.
✓ Branch 5 taken 12411148 times.
✓ Branch 6 taken 1479536 times.
 27647210 if ( ( (m_fstransient->iteration == 0) && (r_instance.useALEformulation != 1) ) ||
692
3/4
✓ Branch 1 taken 13756526 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12276990 times.
✓ Branch 4 taken 1479536 times.
 13756526 ( (m_fstransient->iteration > 0) && r_instance.useALEformulation ) ) {
693 /*! In this function we compute the term of the matrix that, without ALE does not change in time, so they should go to the matrix 1.
694 *
695 * We enter in the following cases:
696 *
697 * - Before the first iteration, when ALE is not 1 ( which means either no ALE or ALE = 2 )
698 * indexMatrix = 1 if ALE=0,
699 * indexMatrix = 0 if ALE>=2
700 *
701 * - For each iteration, but not before when ALE>0
702 * indexMatrix = 0
703 *
704 * When ALE=0
705 * you enter only before the first iteration.
706 * When ALE = 1
707 * you do not enter before the first iteration.
708 * When ALE = 2
709 * you enter the first iteration and the following
710 */
711
4/4
✓ Branch 0 taken 9243390 times.
✓ Branch 1 taken 3167758 times.
✓ Branch 2 taken 3081130 times.
✓ Branch 3 taken 6162260 times.
 12411148 if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix) ) {
712 // First part of system matrix.
713 //
714 // flagMatrixRHS :
715 // matrix_and_rhs = 0,
716 // only_matrix = 1,
717 // only_rhs = 2,
718 //
719 //================================
720 // matrix
721 //================================
722
2/2
✓ Branch 0 taken 6179816 times.
✓ Branch 1 taken 69072 times.
 6248888 if (m_useSymmetricStress) {
723 // \eta \int \epsilon u^{n+1} : \epsilon v
724 6179816 m_elementMat[indexMatrix]->eps_phi_i_eps_phi_j(m_theta*2*m_viscosity, *m_feVel, 0, 0, dimension());
725 } else {
726 // \eta \int \nabla u^{n+1} \nabla v
727 69072 m_elementMat[indexMatrix]->grad_phi_i_grad_phi_j(m_theta*m_viscosity, *m_feVel, 0, 0, dimension());
728 }
729
730 // \int q \nabla \cdot u^{n+1} and \int p \nabla \cdot v
731 6248888 m_elementMat[indexMatrix]->psi_j_div_phi_i_and_psi_i_div_phi_j (*m_feVel, *m_fePres, 0, m_velocity->numComponent(), -1, -1);
732
733
1/2
✓ Branch 0 taken 6248888 times.
✗ Branch 1 not taken.
 6248888 if ( ! r_instance.quasistatic ) { //normally you pass this if, you do not enter for instance when using stokesLinearElasticityCoupledModel
734 // \dfrac{\rho}{\dt} u^{n+1} v
735
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6248888 times.
 6248888 if ( r_instance.useALEformulation == 2 ) {
736 ✗ m_elementMat[indexMatrix+1]->phi_i_phi_j(0.5*m_bdf->coeffDeriv0()*coef, *m_feVel, 0, 0, dimension()); // indexMatrix + 1 = 1
737 ✗ m_elementMat[indexMatrix ]->phi_i_phi_j(0.5*m_bdf->coeffDeriv0()*coef, *m_feVel, 0, 0, dimension()); // indexMatrix = 0
738 } else {
739 6248888 m_elementMat[indexMatrix]->phi_i_phi_j(m_bdf->coeffDeriv0()*coef, *m_feVel, 0, 0, dimension());
740 }
741 }
742 }
743 }
744
745
5/6
✓ Branch 1 taken 134158 times.
✓ Branch 2 taken 13756526 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 134158 times.
✓ Branch 5 taken 13756526 times.
✓ Branch 6 taken 134158 times.
 13890684 if ( m_fstransient->iteration > 0 || m_forceConvAndStabComputation) {
746
4/4
✓ Branch 0 taken 9303520 times.
✓ Branch 1 taken 4453006 times.
✓ Branch 2 taken 3081130 times.
✓ Branch 3 taken 6222390 times.
 13756526 if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix)) {
747 //================================
748 // Convective Term
749 //================================
750
2/5
✓ Branch 0 taken 733696 times.
✓ Branch 1 taken 6800440 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 7534136 switch(r_instance.NSequationFlag) {
751 733696 case 0:
752 //nothing to do: Stokes' model
753 733696 break;
754
755 6800440 case 1:
756 // Method 1: \rho u \grad(u)
757 // Natural condition: \sigma(u,p) \cdot n = 0
758 // DEFAULT METHOD
759
760 // Temam's stabilization trick
761 6800440 m_elemFieldAdv.setValue(m_solExtrapol, *m_feVel, iel, m_iVelocity, m_ao, dof());
762 6800440 m_elementMat[0]->div_u_phi_j_phi_i(0.5*m_density,m_elemFieldAdv,*m_feVel,0,0,dimension());
763
764
2/2
✓ Branch 0 taken 6114730 times.
✓ Branch 1 taken 685710 times.
 6800440 if(r_instance.useALEformulation) {
765 //=========
766 // ALE
767 //=========
768 // m_elemFieldAdv -> m_beta where m_beta = u^{n-1} - w^{n} is defined in userElementInit()
769 // m_elemFielVelMesh -> -w^{n}
770 6114730 m_elemFieldAdv.setValue(m_beta, *m_feVel, iel, m_iVelocity, m_ao, dof());
771 6114730 m_elemFieldVelMesh.setValue(externalVec(0), *m_feVel, iel, m_iDisplacement, m_externalAO[0], *m_externalDof[0]);
772
773
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6114730 times.
 6114730 if ( r_instance.useALEformulation == 2 )
774 ✗ m_elementMat[0]->div_u_phi_j_phi_i(0.5*m_density,m_elemFieldVelMesh,*m_feVel,0,0,dimension());
775 else
776 6114730 m_elementMat[0]->div_u_phi_j_phi_i(m_density,m_elemFieldVelMesh,*m_feVel,0,0,dimension());
777 }
778
779 6800440 m_elementMat[0]->u_grad_phi_j_phi_i(m_theta*m_density,m_elemFieldAdv,*m_feVel,0,0,dimension());
780 6800440 break;
781
782 ✗ case 2:
783 ✗ m_elemFieldAdv.setValue(m_solExtrapol, *m_feVel, iel, m_iVelocity, m_ao, dof());
784 // Method 2: rotational formulation
785 ✗ m_elemFieldTotalPressureFormulation.setValue(sequentialSolution(), *m_feVel, iel, m_iVelocity, m_ao, dof());
786
787 // Natural condition: \sigma(u,p) \cdot n + \rho * 0.5 * |u|^2 \cdot n = 0 (total pressure)
788 ✗ m_elementMat[0]->convective_rot(m_density,m_elemFieldTotalPressureFormulation,*m_feVel,0,0,dimension());
789 ✗ break;
790
791 ✗ case 3:
792 // Method 3: Total pressure skew-symmetric formulation
793 ✗ m_elemFieldAdv.setValue(m_solExtrapol, *m_feVel, iel, m_iVelocity, m_ao, dof());
794 // \int_{Sigma} u^{n+1} Grad u^{n} \dot v - v Grad u^{n} \dot u^{n+1}
795 ✗ m_elementMat[0]->phi_grad_u_minus_grad_u_phi(m_density, m_elemFieldAdv, *m_feVel, 0, 0, dimension());
796 ✗ break;
797
798 ✗ default:
799 ✗ FEL_ERROR("Error: No Convection Method found for Navier-Stokes equations");
800 }
801
802 //================================
803 // FE stabilization
804 //================================
805
1/2
✓ Branch 0 taken 7534136 times.
✗ Branch 1 not taken.
 7534136 if(r_instance.NSStabType == 0) { // SUPG (1: face-oriented)
806
2/2
✓ Branch 2 taken 7480596 times.
✓ Branch 3 taken 53540 times.
 7534136 if ( m_velocity->finiteElementType() == m_pressure->finiteElementType()) {
807 7480596 double stab1 = r_instance.stabSUPG;
808 7480596 double stab2 = r_instance.stabdiv;
809 7480596 int type = r_instance.typeSUPG;
810 double Rem_elem,tau;
811 double delta;
812
813
4/4
✓ Branch 0 taken 256000 times.
✓ Branch 1 taken 7224596 times.
✓ Branch 2 taken 51700 times.
✓ Branch 3 taken 204300 times.
 7480596 if ( r_instance.hasImmersedData && r_instance.hasDivDivStab ) {
814
815 // Boilevin-Kayl 02/2018: for the moment, delta is the first external vector which has been pushed back in NSModel for the std::set of methods (hasImmersedData + hasDivDivStab + !useALEformulation)
816 // This should be updated when additional methods will be considered
817 // Boilevin-Kayl 11/2018: the following has been updated in order to handle the case where hasImmersedData and useALEformulation are used at the same time
818
819 // here iel is already relative to the global numbering, so no need to use an additional mapping or a global index like ielG
820
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 51700 times.
 51700 if ( r_instance.useALEformulation ) {
821 ✗ externalVec(1).getValues(1, &iel, &delta);
822 }
823 else {
824
1/2
✓ Branch 2 taken 51700 times.
✗ Branch 3 not taken.
 51700 externalVec(0).getValues(1, &iel, &delta);
825 }
826
827 51700 stab1 *= delta;
828 51700 stab2 *= 1.0/delta;
829 }
830
831 // Beware: if elemFieldRHS is modified in userNS, it wont be considered here (usually elemFieldRHS = 0 (no body force))...
832 // if the matrix and the rhs are not assembled together
833
834 //todo Jeremie 01 02 2012: assemble multiple matrix
835
1/2
✓ Branch 3 taken 7480596 times.
✗ Branch 4 not taken.
 14961192 m_elementMat[0]->stab_supg(m_fstransient->timeStep, stab1, stab2, m_density, m_viscosity, *m_feVel,
836 7480596 m_elemFieldAdv, m_elemFieldRHS, *m_elementVector[0], Rem_elem, tau, type);
837 }
838 }
839 }
840
841
6/8
✓ Branch 0 taken 13756526 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 6222390 times.
✓ Branch 3 taken 7534136 times.
✓ Branch 4 taken 3141260 times.
✓ Branch 5 taken 3081130 times.
✓ Branch 6 taken 3141260 times.
✗ Branch 7 not taken.
 13756526 if( (r_instance.NSStabType == 0) && (flagMatrixRHS==FlagMatrixRHS::only_rhs) && (m_meshUpdateState==FlagMeshUpdated::PreviousMesh) && !m_forceConvAndStabComputation) {
842
1/2
✓ Branch 2 taken 3141260 times.
✗ Branch 3 not taken.
 3141260 if ( m_velocity->finiteElementType() == m_pressure->finiteElementType()) {
843
1/2
✓ Branch 2 taken 3141260 times.
✗ Branch 3 not taken.
 3141260 m_elemFieldAdv.setValue(m_solExtrapol, *m_feVel, iel, m_iVelocity, m_ao, dof());
844 3141260 double stab1 = r_instance.stabSUPG;
845 3141260 int type = r_instance.typeSUPG;
846 double Rem_elem,tau;
847
848
1/2
✓ Branch 3 taken 3141260 times.
✗ Branch 4 not taken.
 6282520 m_elementVector[0]->stab_supg(m_fstransient->timeStep, stab1, m_density, m_viscosity, *m_feVel,
849 3141260 m_elemFieldAdv, m_elemFieldRHS, Rem_elem, tau, type);
850 }
851 }
852
853
6/6
✓ Branch 0 taken 9303520 times.
✓ Branch 1 taken 4453006 times.
✓ Branch 2 taken 6222390 times.
✓ Branch 3 taken 3081130 times.
✓ Branch 4 taken 3141260 times.
✓ Branch 5 taken 3081130 times.
 13756526 if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || ((flagMatrixRHS == FlagMatrixRHS::only_rhs) &&(m_meshUpdateState==FlagMeshUpdated::PreviousMesh))){
854 //================================
855 // RHS
856 //================================
857 // 1/ bdf (with order k from data)
858 // \frac{\rho}{\dt} \sum_{i=1}^k \alpha_i u_{n+1-i}
859
1/2
✓ Branch 0 taken 7594266 times.
✗ Branch 1 not taken.
 7594266 if ( !r_instance.quasistatic ) {
860 7594266 m_elemFieldRHSbdf.setValue(m_seqBdfRHS, *m_feVel, iel, m_iVelocity, m_ao, dof());
861
2/2
✓ Branch 0 taken 4896000 times.
✓ Branch 1 taken 2698266 times.
 7594266 if ( r_instance.nonLinearFluid ){
862
6/6
✓ Branch 0 taken 4569600 times.
✓ Branch 1 taken 326400 times.
✓ Branch 2 taken 3033600 times.
✓ Branch 3 taken 1536000 times.
✓ Branch 4 taken 326400 times.
✓ Branch 5 taken 3033600 times.
 4896000 if ( (r_instance.useALEformulation && m_meshUpdateState==FlagMeshUpdated::PreviousMesh) || (!r_instance.useALEformulation) )
863 1862400 m_elementVector[0]->source(-1.*m_density,*m_feVel,m_elemFieldRHSbdf,0,dimension());
864 }
865 else {
866 2698266 m_elementVector[0]->source(m_density,*m_feVel,m_elemFieldRHSbdf,0,dimension());
867 }
868 }
869 // 'm_density' and not 'coef' because the time step is integrated into the bdf term
870
871 // 2/ no bdf method
872 // \dfrac{\rho}{\dt} u^{n} v
873 //m_elemFieldAdv.setValue(m_seqSol, *m_feVel, iel, m_iVelocity, m_ao, dof());
874 //m_elementVector[0]->source(coef,*m_feVel,m_elemFieldAdv,0,dimension());
875
876 // theta-method RHS
877
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7594266 times.
 7594266 if ( m_theta < 1 ) {
878 ✗ computeRHSThetaMethod(iel);
879 }
880 }
881 }
882 }
883 13890684 }
884
885 /***********************************************************************************/
886 /***********************************************************************************/
887
888 ✗ void LinearProblemNS::computeElementArrayCharact(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /* elemIdPoint */, felInt& iel, ElementType& eltType, felInt& ielGeo, FlagMatrixRHS /* flagMatrixRHS */)
889 {
890
891 // Retrieve instance
892 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
893
894 ✗ m_feVel->updateFirstDeriv(0, elemPoint);
895 ✗ m_fePres->updateFirstDeriv(0, elemPoint);
896 ✗ PetscScalar coef = m_density/m_fstransient->timeStep;
897
898 ✗ if ( m_fstransient->iteration == 0) {
899 //================================
900 // matrix
901 //================================
902
903 ✗ if (m_useSymmetricStress) {
904 // \eta \int \epsilon u^{n+1} : \epsilon v
905 ✗ m_elementMat[1]->eps_phi_i_eps_phi_j(2*m_viscosity,*m_feVel, 0, 0, dimension());
906 } else {
907 // \eta \int \nabla u^{n+1} \nabla v
908 ✗ m_elementMat[1]->grad_phi_i_grad_phi_j(m_viscosity, *m_feVel, 0, 0, dimension());
909 }
910
911 // \int q \nabla \cdot u^{n+1} and \int p \nabla \cdot v
912 ✗ m_elementMat[1]->psi_j_div_phi_i_and_psi_i_div_phi_j (*m_feVel,*m_fePres,0,m_velocity->numComponent(), -1, -1);
913
914 // \dfrac{\rho}{\dt} u^{n+1} v
915 ✗ FEL_CHECK(r_instance.orderBdfNS == 1, "Error: Higher bdf order not yet implemented for method of characteristics");
916 ✗ if (r_instance.orderBdfNS == 1)
917 ✗ m_elementMat[1]->phi_i_phi_j(coef,*m_feVel,0,0,dimension());
918 } else {
919 //================================
920 // RHS
921 //================================
922
923 //method of characteristics: \dfrac{\rho}{\dt} u^n o X^n v
924 ✗ if (r_instance.orderBdfNS == 1) {
925 ✗ setValueCharact(sequentialSolution(),m_elemFieldCharact, *m_feVel, iel, m_iVelocity, eltType, ielGeo, m_fstransient->timeStep);
926 ✗ assert(!m_elementVector.empty()); // TODO D.C. FELASSERT?
927 ✗ m_elementVector[0]->source(coef, *m_feVel, m_elemFieldCharact, 0, dimension());
928 }
929
930 //================================
931 // FE stabilization
932 //================================
933
934 ✗ if ( m_velocity->finiteElementType() == m_pressure->finiteElementType()) {
935 ✗ double stab1 = r_instance.stabSUPG;
936 ✗ double stab2 = r_instance.stabdiv;
937 ✗ int type = r_instance.typeSUPG;
938 double Rem_elem,tau;
939
940 // Beware: if elemFieldRHS is modified in userNS, it wont be considered here...
941 // (usually elemFieldRHS = 0 (no body force))
942
943 //todo Jeremie 01 02 2012: assemble multiple matrix
944 ✗ m_elementMat[0]->stab_supg(m_fstransient->timeStep,
945 ✗ stab1,stab2, m_density, m_viscosity, *m_feVel,
946 ✗ m_elemFieldAdv, m_elemFieldRHS,
947 ✗ *m_elementVector[0], Rem_elem, tau, type);
948 }
949 }
950 }
951
952 /***********************************************************************************/
953 /***********************************************************************************/
954
955 ✗ void LinearProblemNS::setValueCharact(PetscVector& v, ElementField& elemFieldCharact, CurrentFiniteElement& fe, felInt iel,
956 int idVariable, const ElementType& eltType, felInt ielGeo, const double dt)
957 {
958 // Retrieve instance
959 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
960
961 Point pt, localCoord, vel;
962 ✗ ElementField elemFieldVelOnDep;
963 ElementType eltTypeArr;
964 felInt ielGeoArr;
965 ✗ double dtLeft, dtTol = r_instance.timeStepCharact;
966 ✗ int numComp = elemFieldCharact.numComp(), numStep = 1 + static_cast<int>(dt/dtTol), effStep;
967
968 ✗ switch (elemFieldCharact.type()) {
969 ✗ case QUAD_POINT_FIELD:
970 ✗ if ( !fe.hasQuadPtCoor() ) {
971 ✗ fe.computeCurrentQuadraturePoint();
972 }
973 //velocity field in the departure element, in the ordering of felisce
974 ✗ elemFieldVelOnDep.initialize(DOF_FIELD, fe, numComp);
975 ✗ elemFieldVelOnDep.setValue(v, fe, iel, idVariable, m_ao, dof());
976
977 //loop on quadrature points of the departure element
978 ✗ for(int ig=0; ig<fe.numQuadraturePoint(); ig++) {
979 //coordinates and local coordinates of each quadrature point pt
980 ✗ pt.setCoor( fe.currentQuadPoint[ig].coor() );
981 ✗ localCoord = fe.quadratureRule().quadraturePoint(ig);
982
983 //velocity at pt: v(pt) = \sum v(x_iDof) phi_iDof(localCoord)
984 ✗ vel.x() = vel.y() = vel.z() = 0.;
985 ✗ for(int iComp=0; iComp < numComp; iComp++) {
986 ✗ for (int iDof=0; iDof < fe.numDof(); iDof++) {
987 ✗ vel[iComp] += elemFieldVelOnDep.val(iComp, iDof) * fe.phi[ig](iDof);
988 }
989 }
990
991 //backtrack the characteristic line a length dt
992 ✗ eltTypeArr = eltType;
993 ✗ ielGeoArr = ielGeo;
994
995 ✗ if (r_instance.characteristicMethod == 1) {
996 //RK4 scheme
997 ✗ if ( !(charactLine(v, idVariable, numComp, dt, numStep, effStep, pt, localCoord, eltTypeArr, ielGeoArr, vel)) ) {
998 ✗ dtLeft = dt - effStep * (dt/numStep);
999 ✗ while ( charactLineElem(v, idVariable, numComp, dtTol, pt, localCoord, eltTypeArr, ielGeoArr, vel, dtLeft) );
1000 }
1001 } else {
1002 //Euler scheme
1003 ✗ dtLeft = dt;
1004 ✗ while ( charactLineElem(v, idVariable, numComp, dtTol, pt, localCoord, eltTypeArr, ielGeoArr, vel, dtLeft) );
1005 }
1006
1007 //compute elemFieldCharact
1008 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1009 ✗ elemFieldCharact.val(iComp, ig) = vel[iComp];
1010 }
1011 }
1012 ✗ break;
1013 ✗ case CONSTANT_FIELD:
1014 case DOF_FIELD:
1015 ✗ FEL_ERROR("setValueCharact from Petsc std::vector is not yet implemented for this type of element field");
1016 ✗ break;
1017 }
1018 }
1019
1020 /***********************************************************************************/
1021 /***********************************************************************************/
1022
1023 ✗ int LinearProblemNS::charactLine(PetscVector& v, int idVariable, int numComp, const double dt, int numStep, int& effStep,
1024 Point& pt, Point& localCoord, ElementType& eltType, felInt& ielGeo, Point& vel)
1025 {
1026 //given a point pt, backtrack the characteristic line a length dt
1027 ✗ double step = - dt / numStep;
1028 ✗ for(int iStep=1; iStep <= numStep; iStep++) {
1029 ✗ if ( !nextPt(v, idVariable, numComp, step, pt, localCoord, eltType, ielGeo, vel) ) {
1030 /* This case occurs when for example the point pt got at iStep is outside the domain,
1031 then the number of effectively done steps before stopping is effStep = iStep-1 and
1032 the output pt, localCoord... are those at (iStep-1)-th step.
1033 */
1034 ✗ effStep = iStep-1;
1035 ✗ return(0);
1036 }
1037 }
1038 ✗ effStep = numStep;
1039
1040 ✗ return(1);
1041 }
1042
1043 /***********************************************************************************/
1044 /***********************************************************************************/
1045
1046 ✗ int LinearProblemNS::nextPt(PetscVector& v, int idVariable, int numComp, const double step,
1047 Point& pt, Point& localCoord, ElementType& eltType, felInt& ielGeo, Point& vel)
1048 {
1049 //given a point pt and initial velocity vel, backtrack the characteristic line a length step (<0) by a variant of RK4 scheme
1050 //input: point pt, its local coordinates, element containing pt, velocity at pt
1051 //output: next point pt with local coordinates localCoord, element containing pt, vel is velocity at pt
1052 Point ptTmp, localCoordTmp;
1053 ElementType eltTypeTmp;
1054 felInt ielGeoTmp;
1055 ✗ GeometricMeshRegion::seedElement seedElt, seedEltTmp;
1056 ✗ double s1 = step / 6.;
1057 int base;
1058
1059 std::array<Point, 4> veltmp;
1060 ✗ seedElt = std::make_pair(eltType, ielGeo);
1061
1062 // first, second and third integration point, element containing this point, velocity at this point
1063 ✗ veltmp[0] = vel;
1064 ✗ for (int intPnt = 0; intPnt < 3; ++intPnt) {
1065
1066 ✗ for(int iComp = 0; iComp < numComp; iComp++)
1067 ✗ ptTmp[iComp] = pt[iComp] + 0.5 * step * veltmp[intPnt][iComp];
1068
1069 ✗ base = ++m_pLocator->getMark();
1070
1071 ✗ seedEltTmp = m_pLocator->localizePoint(seedElt, ptTmp, localCoordTmp, base);
1072
1073 ✗ if ( seedEltTmp.first != GeometricMeshRegion::ELEMENT_TYPE_COUNTER && seedEltTmp.second >= 0 ) {
1074 ✗ eltTypeTmp = seedEltTmp.first;
1075 ✗ ielGeoTmp = seedEltTmp.second;
1076
1077 ✗ seedElt = seedEltTmp;
1078 } else {
1079
1080 ✗ return 0;
1081 }
1082
1083 ✗ veltmp[intPnt+1] = vecInterp(v, idVariable, numComp, localCoordTmp, eltTypeTmp, ielGeoTmp);
1084 }
1085
1086 // final point
1087 ✗ for(int iComp=0; iComp < numComp; iComp++)
1088 ✗ ptTmp[iComp] = pt[iComp] + s1 * ( veltmp[0][iComp] + 2*(veltmp[1][iComp] + veltmp[2][iComp]) + veltmp[3][iComp] );
1089
1090 ✗ base = ++m_pLocator->getMark();
1091
1092 ✗ seedEltTmp = m_pLocator->localizePoint(seedElt, ptTmp, localCoordTmp, base);
1093
1094 ✗ if ( seedEltTmp.first != GeometricMeshRegion::ELEMENT_TYPE_COUNTER && seedEltTmp.second >= 0 ) {
1095 ✗ eltTypeTmp = seedEltTmp.first;
1096 ✗ ielGeoTmp = seedEltTmp.second;
1097
1098 ✗ seedElt = seedEltTmp;
1099 } else {
1100
1101 ✗ return 0;
1102 }
1103
1104 // if succeed, update pt, localCoord, vel
1105 ✗ pt = ptTmp;
1106 ✗ localCoord = localCoordTmp;
1107
1108 ✗ vel = vecInterp(v, idVariable, numComp, localCoord, eltType, ielGeo);
1109
1110 ✗ return 1;
1111 }
1112
1113 /***********************************************************************************/
1114 /***********************************************************************************/
1115
1116 ✗ int LinearProblemNS::charactLineElem(PetscVector& v, int idVariable, int numComp, double dtTol, Point& pt, Point& localCoord, ElementType& eltType, felInt& ielGeo, Point& vel, double& dtLeft)
1117 {
1118 //given a point pt, backtrack the characteristic line in one element a length dtTol at most
1119 //update pt, localCoord, eltType, ielGeo, vel and the time remaining dtLeft
1120 //return 0 if dtLeft = 0 or reach to the boundary of the domain, else return 1
1121 const GeoElement *geoEle;
1122 ✗ geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
1123 ✗ TypeShape eltShape = (TypeShape)(geoEle->shape().typeShape());
1124
1125 Point pt1, localCoord1;
1126 double barCoord[4], vecBarCoord[4];
1127 ✗ double dtAdvec = 0.;
1128 ✗ ElementType eltTypeTmp = eltType;
1129 ✗ felInt ielGeoTmp = ielGeo;
1130 int i0;//, i1, i2;
1131
1132 ✗ if (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2] < 1e-20) //almost null velocity, no advection
1133 ✗ return(0);
1134
1135 ✗ switch (eltShape) {
1136
1137 ✗ case Triangle: {
1138 ✗ barCoord[2] = localCoord[0];
1139 ✗ barCoord[0] = localCoord[1];
1140 ✗ barCoord[1] = 1. - localCoord[0] - localCoord[1];
1141 //endpoint of velocity std::vector starting from pt
1142 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1143 ✗ pt1[iComp] = pt[iComp] + vel[iComp];
1144 }
1145 //local coordinates of pt1, barycentric coordinates of vel
1146 ✗ m_mesh[m_currentMesh]->findLocalCoord(eltTypeTmp, ielGeoTmp, pt1, localCoord1);
1147 ✗ vecBarCoord[2] = localCoord1[0] - localCoord[0];
1148 ✗ vecBarCoord[0] = localCoord1[1] - localCoord[1];
1149 ✗ vecBarCoord[1] = - (vecBarCoord[0] + vecBarCoord[2]);
1150
1151 //dtAdvec = time length from pt to the first intersection of [pt pt1] with an edge
1152 // = min{barCoord[i]/vecBarCoord[i]} s.t. vecBarCoord[i]>0 and not too small
1153 ✗ dtAdvec = dtLeft * 1e+6;
1154 ✗ i0 = 0;
1155 ✗ for (int ind=0; ind < 3; ind++) {
1156 ✗ if (vecBarCoord[ind] > 1e-10) {
1157 ✗ if (dtAdvec > barCoord[ind]/vecBarCoord[ind]) {
1158 ✗ dtAdvec = barCoord[ind]/vecBarCoord[ind];
1159 ✗ i0 = ind;
1160 }
1161 }
1162 }
1163
1164 //case when advection stays in the current element, update only pt, localCoord, vel, dtLeft
1165 ✗ if ( (dtLeft <= dtAdvec) && (dtLeft <= dtTol) ) {
1166 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1167 ✗ pt[iComp] = pt[iComp] - dtLeft * vel[iComp];
1168 ✗ localCoord[iComp] = localCoord[iComp] - dtLeft * vecBarCoord[(iComp+2)%3];
1169 }
1170 ✗ vel = vecInterp(v, idVariable, numComp, localCoord, eltType, ielGeo);
1171 ✗ dtLeft = 0.;
1172 ✗ return(0);
1173 }
1174
1175 //case when advection stays in the current element but dtLeft > dtTol
1176 ✗ if ( (dtLeft <= dtAdvec) && (dtLeft > dtTol) ) {
1177 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1178 ✗ pt[iComp] = pt[iComp] - dtTol * vel[iComp];
1179 ✗ localCoord[iComp] = localCoord[iComp] - dtTol * vecBarCoord[(iComp+2)%3];
1180 }
1181 ✗ vel = vecInterp(v, idVariable, numComp, localCoord, eltType, ielGeo);
1182 ✗ dtLeft -= dtTol;
1183 ✗ return(1);
1184 }
1185
1186 //case when advection stays in the current element but dtAdvec > dtTol
1187 ✗ if ( (dtLeft > dtAdvec) && (dtAdvec >= dtTol) ) {
1188 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1189 ✗ pt[iComp] = pt[iComp] - dtTol * vel[iComp];
1190 ✗ localCoord[iComp] = localCoord[iComp] - dtTol * vecBarCoord[(iComp+2)%3];
1191 }
1192 ✗ vel = vecInterp(v, idVariable, numComp, localCoord, eltType, ielGeo);
1193 ✗ dtLeft -= dtTol;
1194 ✗ return(1);
1195 }
1196
1197 //case when advection goes out of the current element through the edge i0
1198 ✗ if ( (dtLeft > dtAdvec) && (dtAdvec < dtTol) ) {
1199 //update pt, vel, dtLeft (not yet localCoord)
1200 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1201 ✗ pt[iComp] = pt[iComp] - dtAdvec * vel[iComp]; //this point is on the edge i0 of the current element
1202 ✗ localCoord[iComp] = localCoord[iComp] - dtAdvec * vecBarCoord[(iComp+2)%3];
1203 }
1204 ✗ vel = vecInterp(v, idVariable, numComp, localCoord, eltType, ielGeo);
1205 ✗ dtLeft -= dtAdvec;
1206 //update element = edge i0 -> adjacency and localCoord in updated element
1207 ✗ if ( m_mesh[m_currentMesh]->getAdjElement(eltType, ielGeo, i0) ) { //element is already updated to the adjacent one
1208 ✗ eltTypeTmp = eltType;
1209 ✗ ielGeoTmp = ielGeo;
1210 ✗ m_mesh[m_currentMesh]->findLocalCoord(eltTypeTmp, ielGeoTmp, pt, localCoord);
1211 ✗ return(1);
1212 } else //reach to the boundary of the domain
1213 ✗ return(0);
1214 }
1215 ✗ break;
1216 }
1217
1218 ✗ case Tetrahedron: {
1219 ✗ barCoord[3] = localCoord[0];
1220 ✗ barCoord[1] = localCoord[1];
1221 ✗ barCoord[0] = localCoord[2];
1222 ✗ barCoord[2] = 1. - localCoord[0] - localCoord[1] - localCoord[2];
1223 //endpoint of velocity std::vector starting from pt
1224 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1225 ✗ pt1[iComp] = pt[iComp] + vel[iComp];
1226 }
1227 //local coordinates of pt1, barycentric coordinates of vel
1228 ✗ m_mesh[m_currentMesh]->findLocalCoord(eltTypeTmp, ielGeoTmp, pt1, localCoord1);
1229 ✗ vecBarCoord[3] = localCoord1[0] - localCoord[0];
1230 ✗ vecBarCoord[1] = localCoord1[1] - localCoord[1];
1231 ✗ vecBarCoord[0] = localCoord1[2] - localCoord[2];
1232 ✗ vecBarCoord[2] = - (vecBarCoord[0] + vecBarCoord[1] + vecBarCoord[3]);
1233
1234 //dtAdvec = time length from pt to the first intersection of [pt pt1] with a face
1235 // = min{barCoord[i]/vecBarCoord[i]} s.t. vecBarCoord[i]>0 and not too small
1236 ✗ dtAdvec = dtLeft * 1e+6;
1237 ✗ i0 = 0;
1238 ✗ for (int ind=0; ind < 4; ind++) {
1239 ✗ if (vecBarCoord[ind] > 1e-10) {
1240 ✗ if (dtAdvec > barCoord[ind]/vecBarCoord[ind]) {
1241 ✗ dtAdvec = barCoord[ind]/vecBarCoord[ind];
1242 ✗ i0 = ind;
1243 }
1244 }
1245 }
1246
1247 //case when advection stays in the current element, update only pt, localCoord, vel, dtLeft
1248 ✗ if ( (dtLeft <= dtAdvec) && (dtLeft <= dtTol) ) {
1249 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1250 ✗ pt[iComp] = pt[iComp] - dtLeft * vel[iComp];
1251 }
1252 ✗ localCoord[0] = localCoord[0] - dtLeft * vecBarCoord[3];
1253 ✗ localCoord[1] = localCoord[1] - dtLeft * vecBarCoord[1];
1254 ✗ localCoord[2] = localCoord[2] - dtLeft * vecBarCoord[0];
1255 ✗ vel = vecInterp(v, idVariable, numComp, localCoord, eltType, ielGeo);
1256 ✗ dtLeft = 0.;
1257 ✗ return(0);
1258 }
1259
1260 //case when advection stays in the current element but dtLeft > dtTol
1261 ✗ if ( (dtLeft <= dtAdvec) && (dtLeft > dtTol) ) {
1262 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1263 ✗ pt[iComp] = pt[iComp] - dtTol * vel[iComp];
1264 }
1265 ✗ localCoord[0] = localCoord[0] - dtTol * vecBarCoord[3];
1266 ✗ localCoord[1] = localCoord[1] - dtTol * vecBarCoord[1];
1267 ✗ localCoord[2] = localCoord[2] - dtTol * vecBarCoord[0];
1268 ✗ vel = vecInterp(v, idVariable, numComp, localCoord, eltType, ielGeo);
1269 ✗ dtLeft -= dtTol;
1270 ✗ return(1);
1271 }
1272
1273 //case when advection stays in the current element but dtAdvec > dtTol
1274 ✗ if ( (dtLeft > dtAdvec) && (dtAdvec >= dtTol) ) {
1275 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1276 ✗ pt[iComp] = pt[iComp] - dtTol * vel[iComp];
1277 }
1278 ✗ localCoord[0] = localCoord[0] - dtTol * vecBarCoord[3];
1279 ✗ localCoord[1] = localCoord[1] - dtTol * vecBarCoord[1];
1280 ✗ localCoord[2] = localCoord[2] - dtTol * vecBarCoord[0];
1281 ✗ vel = vecInterp(v, idVariable, numComp, localCoord, eltType, ielGeo);
1282 ✗ dtLeft -= dtTol;
1283 ✗ return(1);
1284 }
1285
1286 //case when advection goes out of the current element through the face i0
1287 ✗ if ( (dtLeft > dtAdvec) && (dtAdvec < dtTol) ) {
1288 //update pt, vel, dtLeft (not yet localCoord)
1289 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1290 ✗ pt[iComp] = pt[iComp] - dtAdvec * vel[iComp]; //this point is on the face i0 of the current element
1291 }
1292 ✗ localCoord[0] = localCoord[0] - dtAdvec * vecBarCoord[3];
1293 ✗ localCoord[1] = localCoord[1] - dtAdvec * vecBarCoord[1];
1294 ✗ localCoord[2] = localCoord[2] - dtAdvec * vecBarCoord[0];
1295 ✗ vel = vecInterp(v, idVariable, numComp, localCoord, eltType, ielGeo);
1296 ✗ dtLeft -= dtAdvec;
1297 //update element = face i0 -> adjacency and localCoord in updated element
1298 ✗ if ( m_mesh[m_currentMesh]->getAdjElement(eltType, ielGeo, i0) ) { //element is already updated to the adjacent one
1299 ✗ eltTypeTmp = eltType;
1300 ✗ ielGeoTmp = ielGeo;
1301 ✗ m_mesh[m_currentMesh]->findLocalCoord(eltTypeTmp, ielGeoTmp, pt, localCoord);
1302 ✗ return(1);
1303 } else //reach to the boundary of the domain
1304 ✗ return(0);
1305 }
1306 ✗ break;
1307 }
1308
1309 ✗ case NullShape:
1310 case Node:
1311 case Segment:
1312 case Quadrilateral:
1313 case Hexahedron:
1314 case Prism:
1315 case Pyramid: {
1316 ✗ FEL_ERROR("Unknown (or not implemented) element Shape");
1317 ✗ break;
1318 }
1319 }
1320
1321 ✗ return(0);
1322 }
1323
1324
1325 /***********************************************************************************/
1326 /***********************************************************************************/
1327
1328 ✗ Point LinearProblemNS::vecInterp(PetscVector& v, int idVariable, int numComp,
1329 const Point& localCoord, ElementType& eltType, felInt ielGeo)
1330 {
1331 //given a point in an element with its local coordinates, compute the interpolated std::vector at this point
1332 ✗ Point vecInt(0.);
1333 ✗ ElementField elemField;
1334
1335 //define finite element
1336 const RefElement *refEle;
1337 const GeoElement *geoEle;
1338 ✗ int typeOfFiniteElement = 0;
1339 CurrentFiniteElement* fe;
1340 ✗ int numPointPerElt = 0;
1341 ✗ std::vector<Point*> elemPoint;
1342 ✗ std::vector<felInt> elemIdPoint;
1343 felInt ielSupportDof;
1344
1345 ✗ geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
1346 ✗ typeOfFiniteElement = m_listVariable[idVariable].finiteElementType();
1347 ✗ refEle = geoEle->defineFiniteEle(eltType,typeOfFiniteElement, *m_mesh[m_currentMesh]);
1348 ✗ fe = new CurrentFiniteElement(*refEle,*geoEle,m_listVariable[idVariable].degreeOfExactness());
1349
1350 //resize array
1351 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
1352 ✗ elemPoint.resize(numPointPerElt, nullptr);
1353 ✗ elemIdPoint.resize(numPointPerElt, 0);
1354 ✗ setElemPoint(eltType, ielGeo, elemPoint, elemIdPoint, &ielSupportDof);
1355
1356 //interpolation
1357 ✗ elemField.initialize(DOF_FIELD, *fe, numComp);
1358 ✗ elemField.setValue(v, *fe, ielSupportDof, idVariable, m_ao, dof());
1359 ✗ for(int iComp=0; iComp < numComp; iComp++) {
1360 ✗ for (int iDof=0; iDof < fe->numDof(); iDof++) {
1361 ✗ vecInt[iComp] += elemField.val(iComp, iDof) * refEle->basisFunction().phi(iDof, localCoord);
1362 }
1363 }
1364
1365 ✗ delete fe;
1366 ✗ elemPoint.clear();
1367 ✗ elemIdPoint.clear();
1368 ✗ return(vecInt);
1369 }
1370
1371 /***********************************************************************************/
1372 /***********************************************************************************/
1373
1374 2508 void LinearProblemNS::allocateElemFieldBoundaryConditionDerivedLinPb(const int idBCforLinCombMethod)
1375 {
1376 // Retrieve instance
1377 2508 auto& r_instance = FelisceParam::instance(instanceIndex());
1378
1379 // implicit lumpedModelBC but applied as Neumann Normal BC
1380
4/6
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 2424 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 84 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2508 times.
 2508 if(m_lumpedModelBC && r_instance.model == "NSlinComb") {
1381 ✗ bool isFullImplicit = true;
1382 ✗ bool isFullExplicit = true;
1383 ✗ for(std::size_t i=0; i<m_lumpedModelBC->size(); i++) {
1384 ✗ if(r_instance.lumpedModelBCAlgo[i] != 2)
1385 ✗ isFullImplicit = false;
1386 ✗ if(r_instance.lumpedModelBCAlgo[i] != 1)
1387 ✗ isFullImplicit = false;
1388 }
1389 ✗ if(!isFullImplicit && !isFullExplicit)
1390 ✗ FEL_WARNING("NSlinComb not implemented with mixed explicit and implicit algorithm for the windkessel");
1391
1392 ✗ if (isFullImplicit) {
1393 BoundaryCondition* BC;
1394 ✗ std::vector <double> valueOfBC(1);
1395 ✗ valueOfBC[0] = 0.;
1396 // 1st step: std::set std::vector m_elemFieldNeumannNormal = 0
1397 ✗ for (std::size_t i=0; i < m_boundaryConditionList.numNeumannNormalBoundaryCondition(); i++) {
1398 ✗ BC = m_boundaryConditionList.NeumannNormal(i);
1399 ✗ switch (BC->typeValueBC()) {
1400 ✗ case Constant:
1401 ✗ m_elemFieldNeumannNormal[i].setValue(valueOfBC);
1402 ✗ break;
1403 ✗ case FunctionT:
1404 ✗ m_elemFieldNeumannNormal[i].setValue(valueOfBC);
1405 ✗ break;
1406 ✗ case Vector:
1407 case FunctionS:
1408 case FunctionTS:
1409 case EnsightFile:
1410 ✗ FEL_ERROR("This kind of boundary conditions are not available.");
1411 ✗ break;
1412 }
1413 }
1414
1415 // 2nd step: std::set m_elemFieldNeumannNormal[idBCforLinCombMethod] = -1 alternatively while preProcessing
1416 ✗ if ( m_fstransient->iteration == 0 ) {
1417 ✗ valueOfBC[0] = -1.;
1418 ✗ BC = m_boundaryConditionList.NeumannNormal(idBCforLinCombMethod);
1419 ✗ switch (BC->typeValueBC()) {
1420 ✗ case Constant:
1421 ✗ m_elemFieldNeumannNormal[idBCforLinCombMethod].setValue(valueOfBC);
1422 ✗ break;
1423 ✗ case FunctionT:
1424 ✗ m_elemFieldNeumannNormal[idBCforLinCombMethod].setValue(valueOfBC);
1425 ✗ break;
1426 ✗ case FunctionS:
1427 case Vector:
1428 case FunctionTS:
1429 case EnsightFile:
1430 ✗ FEL_ERROR("This kind of boundary conditions are not available.");
1431 ✗ break;
1432 }
1433 }
1434 ✗ valueOfBC.clear();
1435 }
1436 }
1437 2508 }
1438
1439 96 void LinearProblemNS::initExtrapol(PetscVector& V_1)
1440 {
1441
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 42 times.
 96 if (allocateSeqVecExt == false) {
1442 54 m_solExtrapol.create(PETSC_COMM_SELF);
1443 54 m_solExtrapol.setType(VECSEQ);
1444 54 m_solExtrapol.setSizes(PETSC_DECIDE,m_numDof);
1445 54 allocateSeqVecExt = true;
1446 }
1447 96 gatherVector(V_1, m_solExtrapol);
1448 96 }
1449
1450 1554 void LinearProblemNS::updateExtrapol(PetscVector& V_1)
1451 {
1452 1554 gatherVector(V_1, m_solExtrapol);
1453 1554 }
1454
1455 1608 void LinearProblemNS::gatherVectorBeforeAssembleMatrixRHS()
1456 {
1457
2/2
✓ Branch 0 taken 54 times.
✓ Branch 1 taken 1554 times.
 1608 if (allocateSeqVec == false) {
1458 54 m_seqBdfRHS.create(PETSC_COMM_SELF);
1459 54 m_seqBdfRHS.setType(VECSEQ);
1460 54 m_seqBdfRHS.setSizes(PETSC_DECIDE,m_numDof);
1461 54 allocateSeqVec = true;
1462 }
1463 1608 gatherVector(m_bdf->vector(), m_seqBdfRHS);
1464 1608 }
1465
1466 /***********************************************************************************/
1467 /***********************************************************************************/
1468
1469 3044 void LinearProblemNS::addMatrixRHS()
1470 {
1471
1472 //! The not constant matrix (0) is summed with the constant matrix (1).
1473 //! The result is stored into the not constant matrix (0)
1474 //! When using ALE there is not a constant matrix (1) and therefore is useless to do it.
1475
1476 // Retrieve instance
1477 3044 auto& r_instance = FelisceParam::instance(instanceIndex());
1478
1479
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3044 times.
 3044 if ( r_instance.useFDlagrangeMult )
1480 ✗ matrix(0).axpy(1,matrix(1),SUBSET_NONZERO_PATTERN);
1481
2/2
✓ Branch 0 taken 2368 times.
✓ Branch 1 taken 676 times.
 3044 else if( r_instance.useALEformulation == 0 )
1482 2368 matrix(0).axpy(1,matrix(1),SAME_NONZERO_PATTERN);
1483 3044 }
1484
1485 /***********************************************************************************/
1486 /***********************************************************************************/
1487
1488 1912 void LinearProblemNS::addMatrixRHSNl(const FlagMatrixRHS flagMatrixRHS)
1489 {
1490 // Retrieve instance
1491 1912 auto& r_instance = FelisceParam::instance(instanceIndex());
1492
1493
2/3
✓ Branch 0 taken 1436 times.
✓ Branch 1 taken 476 times.
✗ Branch 2 not taken.
 1912 switch (flagMatrixRHS) {
1494 1436 case FlagMatrixRHS::only_matrix:
1495 1436 matrix(0).axpy(1, matrix(2), SAME_NONZERO_PATTERN);
1496 1436 break;
1497 476 case FlagMatrixRHS::only_rhs:
1498
1/2
✓ Branch 0 taken 476 times.
✗ Branch 1 not taken.
 476 if (r_instance.useALEformulation)
1499 476 vector().axpy(1, externalVec(1));
1500 476 break;
1501 ✗ default :
1502 ✗ break;
1503 }
1504 1912 }
1505
1506 /***********************************************************************************/
1507 /***********************************************************************************/
1508
1509 5330 void LinearProblemNS::applyEssentialBoundaryConditionDerivedProblem(int rank, FlagMatrixRHS flagMatrixRHS)
1510 {
1511 // Retrieve instance
1512 5330 auto& r_instance = FelisceParam::instance(instanceIndex());
1513
1514
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5330 times.
 5330 if ( r_instance.zeroDisplacementAtBoarders ) {
1515 ✗ if (FelisceParam::verbose() > 2)
1516 ✗ std::cout<<"["<<rank<<"]--Applying zero-Dirichlet BCs on the velocity at the interface between inlet outlet and surface."<<std::endl;
1517 ✗ getRings();
1518 ✗ double bigValue=r_instance.penValueForEmbedFSI*1e10;
1519
1520 ✗ for(auto itDofBC = m_idDofRings.begin(); itDofBC != m_idDofRings.end(); ++itDofBC) {
1521 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix))
1522 ✗ matrix(0).setValue(*itDofBC, *itDofBC, bigValue, INSERT_VALUES);
1523 ✗ if((flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) ||( (flagMatrixRHS == FlagMatrixRHS::only_rhs)
1524 ✗ &&(m_meshUpdateState==FlagMeshUpdated::PreviousMesh) ) )
1525 ✗ vector().setValue(*itDofBC,0.0, INSERT_VALUES);
1526 }
1527 }
1528 5330 }
1529
1530 /***********************************************************************************/
1531 /***********************************************************************************/
1532
1533 42 void LinearProblemNS::applyNaturalBoundaryConditionLumpedModelBC(std::vector<Point*>& /* elemPoint */, const std::vector<felInt>& /* elemIdPoint */, int /* label */, felInt /* numEltPerLabel */, felInt* /* ielSupportDof */, std::size_t iLumpedModelBC, ElementType& /* eltType */, felInt /* numElement_eltType */, FlagMatrixRHS /* flagMatrixRHS */)
1534 {
1535 // Retrieve instance
1536 42 auto& r_instance = FelisceParam::instance(instanceIndex());
1537
1538 // to apply "NeumannNormal" BC with -p0
1539
1/2
✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
 42 if (r_instance.lumpedModelBCAlgo[iLumpedModelBC] == 1)
1540
1/2
✓ Branch 2 taken 42 times.
✗ Branch 3 not taken.
 42 PlumpedModelBC.push_back(-(m_lumpedModelBC->Pp(iLumpedModelBC)));
1541 ✗ else if (r_instance.lumpedModelBCAlgo[iLumpedModelBC] == 2) {
1542 // Add the explicit part on the RHS
1543 ✗ if(r_instance.lumpedModelBCType[iLumpedModelBC] == 1) {
1544 // RCR or R type
1545 ✗ double alpha = r_instance.lumpedModelBC_Rdist[iLumpedModelBC] * r_instance.lumpedModelBC_C[iLumpedModelBC];
1546 ✗ double beta = r_instance.timeStep;
1547 ✗ double Pv = r_instance.lumpedModelBC_Pvenous[iLumpedModelBC];
1548 ✗ double Pd = m_lumpedModelBC->Pd(iLumpedModelBC);
1549 ✗ PlumpedModelBC.push_back(-(alpha*Pd + beta*Pv)/(alpha + beta));
1550 ✗ } else if(r_instance.lumpedModelBCType[iLumpedModelBC] == 2) {
1551 // RC type
1552 ✗ PlumpedModelBC.push_back(-(m_lumpedModelBC->Pd(iLumpedModelBC)));
1553 }
1554 }
1555 42 m_elemFieldLumpedModelBC[iLumpedModelBC].setValue(PlumpedModelBC);
1556
1557
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 42 times.
 42 if (verbosity()>2) {
1558 ✗ std::cout << "LumpedModelBC on label " << m_currentLabel << std::endl;
1559 ✗ std::cout << "Applied PlumpedModelBC = " << -PlumpedModelBC[0] << std::endl;
1560 }
1561
1562 42 PlumpedModelBC.clear();
1563 42 }
1564
1565 /***********************************************************************************/
1566 /***********************************************************************************/
1567
1568 ✗ void LinearProblemNS::assembleMatrixImplLumpedModelBC(int rank, const std::size_t iLabel)
1569 {
1570 // Retrieve instance
1571 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
1572
1573 ✗ if (m_fstransient->iteration == 0) {
1574 ✗ int label = r_instance.lumpedModelBCLabel[iLabel];
1575 ✗ PetscScalar coef = 0.;
1576 ✗ if ( r_instance.lumpedModelBCType[iLabel] == 1 ) {
1577 // R lumpedModel bc, or RCR lumpedModel bc
1578 ✗ coef = r_instance.lumpedModelBC_Rdist[iLabel]/(1.0 + r_instance.lumpedModelBC_Rdist[iLabel] * r_instance.lumpedModelBC_C[iLabel]/r_instance.timeStep) + r_instance.lumpedModelBC_Rprox[iLabel];
1579 ✗ } else if ( r_instance.lumpedModelBCType[iLabel] == 2 ) {
1580 // RC lumpedModel bc
1581 ✗ coef = r_instance.lumpedModelBC_Rprox[iLabel] + r_instance.timeStep / r_instance.lumpedModelBC_C[iLabel];
1582 }
1583
1584 ✗ double *tmp_phi_i_dot_n = new double [m_numDof];
1585 ✗ double *phi_i_dot_n = new double [m_numDof]; // to store globally \int_Gamma phi_i n
1586 ✗ int *tmp_ind = new int [m_numDof];
1587 ✗ int *ind = new int [m_numDof]; // to mark the position of matrix assembling
1588 ✗ felInt node = 0;
1589 ✗ int idVar = m_listVariable.getVariableIdList(velocity);
1590 // int iUnknown = m_listVariable.listIdUnknownOfVariable(idVar);
1591 ✗ int numComp = m_listVariable[idVar].numComponent();
1592 ElementType eltType;
1593 felInt ielGeom;
1594 ✗ int numPointPerElt = 0;
1595 ✗ int currentLabel = 0;
1596 ✗ felInt numEltPerLabel = 0;
1597 ✗ std::vector<Point*> elemPoint;
1598 ✗ std::vector<felInt> elemIdPoint;
1599 felInt ielSupportDof;
1600
1601 // initialize tmp_phi_i_dot_n and phi_i_dot_n to 0
1602 ✗ for (int iDof = 0; iDof < m_numDof; iDof++) {
1603 ✗ tmp_phi_i_dot_n[iDof] = 0.;
1604 ✗ phi_i_dot_n[iDof] = 0.;
1605 ✗ tmp_ind[iDof] = 0;
1606 ✗ ind[iDof] = 0;
1607 }
1608 ✗ const std::vector<ElementType>& bagElementTypeDomainBoundary = m_meshLocal[m_currentMesh]->bagElementTypeDomainBoundary();
1609 ✗ for (std::size_t iTyp = 0; iTyp < bagElementTypeDomainBoundary.size(); ++iTyp) {
1610 ✗ eltType = bagElementTypeDomainBoundary[iTyp];
1611 ✗ ielGeom = 0;
1612 ✗ defineCurvilinearFiniteElement(eltType);
1613 ✗ initElementArrayBD();
1614 ✗ CurvilinearFiniteElement& fe = *m_listCurvilinearFiniteElement[idVar];
1615 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
1616 ✗ elemPoint.resize(numPointPerElt, nullptr);
1617 ✗ elemIdPoint.resize(numPointPerElt, 0);
1618
1619 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin();
1620 ✗ itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
1621 ✗ currentLabel = itRef->first;
1622 ✗ numEltPerLabel = itRef->second.second;
1623 ✗ if (currentLabel == label) {
1624 // compute tmp_phi_i_dot_n = local \int_Gamma phi_i n = local \sum \int_{element of Gamma} phi_i n
1625 ✗ for (felInt iel1 = 0; iel1 < numEltPerLabel; iel1++) {
1626 ✗ setElemPoint(eltType, ielGeom, elemPoint, elemIdPoint, &ielSupportDof);
1627 ✗ fe.updateMeasNormal(0, elemPoint);
1628 ✗ for (int iDof = 0; iDof < fe.numDof(); iDof++) {
1629 ✗ for (int iComp = 0; iComp < numComp; iComp++) {
1630 ✗ dof().loc2glob(ielSupportDof,iDof,idVar,iComp,node);
1631 ✗ tmp_ind[node] = 1;
1632 ✗ for (int ig=0; ig < fe.numQuadraturePoint(); ig++) {
1633 ✗ tmp_phi_i_dot_n[node] += fe.phi[ig](iDof) * fe.normal[ig](iComp) * fe.weightMeas(ig);
1634 }
1635 }
1636 }
1637 ✗ ielGeom++;
1638 } // end of loop on iel1 of each label
1639 } // end of if currentLabel == label
1640 else
1641 ✗ ielGeom += numEltPerLabel;
1642 } //end of loop on itRef
1643 } // end of loop on iTyp
1644
1645 // synchronise value of tmp_phi_i_dot_n and tmp_ind
1646 ✗ MPI_Reduce(&tmp_phi_i_dot_n[0], &phi_i_dot_n[0], m_numDof, MPI_DOUBLE, MPI_SUM,0,MpiInfo::petscComm());
1647 ✗ MPI_Reduce(&tmp_ind[0], &ind[0], m_numDof, MPI_INT, MPI_SUM,0,MpiInfo::petscComm());
1648
1649 ✗ if (rank == 0) {
1650 ✗ felInt cptDof = 0;
1651 ✗ for (int iDof = 0; iDof < m_numDof; iDof++) {
1652 ✗ if (ind[iDof] > 0) {
1653 ✗ tmp_ind[cptDof] = iDof;
1654 ✗ tmp_phi_i_dot_n[cptDof] = phi_i_dot_n[iDof];
1655 ✗ cptDof++;
1656 }
1657 }
1658
1659 ✗ std::vector<double*> valueMat;
1660 ✗ for (unsigned int i = 0; i < numberOfMatrices(); i++) {
1661 ✗ valueMat.push_back(new double[cptDof*cptDof]);
1662 }
1663 ✗ for (int iRow = 0; iRow < cptDof; iRow++) {
1664 ✗ for (int iCol = 0; iCol < cptDof; iCol++) {
1665 ✗ valueMat[0][iCol + iRow*cptDof] = 0.;
1666 ✗ valueMat[1][iCol + iRow*cptDof] = coef * tmp_phi_i_dot_n[iRow] * tmp_phi_i_dot_n[iCol];
1667 }
1668 }
1669 ✗ AOApplicationToPetsc(m_ao,cptDof,tmp_ind);
1670 ✗ for (unsigned int i = 0; i < numberOfMatrices(); i++) {
1671 ✗ matrix(i).setValues(cptDof,tmp_ind,cptDof,tmp_ind,valueMat[i],ADD_VALUES);
1672 ✗ delete [] valueMat[i];
1673 }
1674 ✗ valueMat.clear();
1675 }
1676
1677 ✗ delete [] tmp_phi_i_dot_n;
1678 ✗ delete [] phi_i_dot_n;
1679 ✗ delete [] tmp_ind;
1680 ✗ delete [] ind;
1681 }
1682 }
1683
1684 /***********************************************************************************/
1685 /***********************************************************************************/
1686
1687 54 void LinearProblemNS::userChangePattern(int /* numProc */, int /* rankProc */)
1688 {
1689 // Retrieve instance
1690 54 auto& r_instance = FelisceParam::instance(instanceIndex());
1691
1692 54 bool isImplicitLumpedModelBCused = false;
1693
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 54 times.
 70 for(std::size_t i=0; i<r_instance.lumpedModelBCAlgo.size(); i++) {
1694
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
 16 if(r_instance.lumpedModelBCAlgo[i] == 2)
1695 ✗ isImplicitLumpedModelBCused = true;
1696 }
1697
1698 // ------------------------------------------- //
1699 // monolithic model and implicit LumpedModelBC //
1700 // ------------------------------------------- //
1701
4/8
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 46 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 54 times.
 54 if (r_instance.lumpedModelBCLabel.size() > 0 && isImplicitLumpedModelBCused && r_instance.model == "NS") {
1702 // re-build random repartition of dof on processes
1703 ✗ felInt numDofByProc = m_numDof/MpiInfo::numProc();
1704 ✗ felInt *numDofPerProc = new felInt [MpiInfo::numProc()];
1705 ✗ for (int i = 0; i < MpiInfo::numProc(); i++) {
1706 ✗ if (i == MpiInfo::numProc() - 1)
1707 ✗ numDofPerProc[i] = m_numDof - i*numDofByProc;
1708 else
1709 ✗ numDofPerProc[i] = numDofByProc;
1710 }
1711 ✗ felInt scaleDof = 0;
1712 ✗ for (int i = 0; i < MpiInfo::rankProc(); i++) {
1713 ✗ scaleDof += numDofPerProc[i];
1714 }
1715 ✗ felInt *rankDof = new felInt[m_numDof];
1716 ✗ for (int i = 0; i < m_numDof; i++) {
1717 ✗ rankDof[i] = -1;
1718 }
1719 ✗ for (felInt i = 0; i < numDofPerProc[MpiInfo::rankProc()]; i++) {
1720 ✗ rankDof[i+scaleDof] = MpiInfo::rankProc();
1721 }
1722
1723 // loop to connect all added dof
1724 ✗ std::vector< std::set<felInt> > nodesNeighborhood(m_numDof);
1725 ✗ felInt node1 = 0;
1726 ✗ felInt node2 = 0;
1727 ✗ const int idVar = m_listVariable.getVariableIdList(velocity);
1728 ✗ const int iUnknown = m_listVariable.listIdUnknownOfVariable(idVar);
1729 ✗ const int numComp = m_listVariable[idVar].numComponent();
1730 ✗ const std::size_t iMesh = m_listVariable[idVar].idMesh();
1731 ✗ int searchLabel = -1;
1732 ElementType eltType;
1733 ✗ int currentLabel = 0;
1734 ✗ int numEltPerLabel = 0;
1735 felInt ielSupportDof1, ielSupportDof2;
1736 felInt ielGeom1, ielGeom2, ielGeomTmp;
1737
1738 ✗ for (std::size_t iLabel = 0; iLabel < r_instance.lumpedModelBCLabel.size(); iLabel++) {
1739 ✗ if(r_instance.lumpedModelBCAlgo[iLabel] == 2) {
1740 ✗ searchLabel = r_instance.lumpedModelBCLabel[iLabel];
1741 ✗ const std::vector<ElementType>& bagElementTypeDomainBoundary = m_mesh[iMesh]->bagElementTypeDomainBoundary();
1742 ✗ for (std::size_t iTyp = 0; iTyp < bagElementTypeDomainBoundary.size(); ++iTyp) {
1743 ✗ eltType = bagElementTypeDomainBoundary[iTyp];
1744 ✗ ielGeom1 = 0;
1745
1746 ✗ for(auto itRef = m_mesh[iMesh]->intRefToBegEndMaps[eltType].begin(); itRef != m_mesh[iMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
1747 ✗ currentLabel = itRef->first;
1748 ✗ numEltPerLabel = itRef->second.second;
1749
1750 ✗ if (searchLabel == currentLabel) {
1751 ✗ ielGeomTmp = ielGeom1;
1752 ✗ for (felInt iel1 = 0; iel1 < numEltPerLabel; iel1++) {
1753 // find ielSupportDof1
1754 ✗ supportDofUnknown(iUnknown).getIdElementSupport(eltType, ielGeom1, ielSupportDof1);
1755 ✗ for (int iSupport = 0; iSupport < supportDofUnknown(iUnknown).getNumSupportDof(ielSupportDof1); iSupport++) {
1756 ✗ for (int iComp = 0; iComp < numComp; iComp++) {
1757 ✗ dof().loc2glob(ielSupportDof1,iSupport,idVar,iComp,node1);
1758 ✗ if (rankDof[node1] == MpiInfo::rankProc()) {
1759 ✗ ielGeom2 = ielGeomTmp;
1760 ✗ for (felInt iel2 = 0; iel2 < numEltPerLabel; iel2++) {
1761 ✗ supportDofUnknown(iUnknown).getIdElementSupport(eltType, ielGeom2, ielSupportDof2);
1762 ✗ for (int jSupport = 0; jSupport < supportDofUnknown(iUnknown).getNumSupportDof(ielSupportDof2); jSupport++) {
1763 ✗ for (int jComp = 0; jComp < numComp; jComp++) {
1764 ✗ dof().loc2glob(ielSupportDof2,jSupport,idVar,jComp,node2);
1765 ✗ if ( node1 != node2)
1766 ✗ nodesNeighborhood[node1].insert(node2);
1767 }
1768 }
1769 ✗ ielGeom2++;
1770 } // end of loop on iel2 of each label
1771 } // end of if rankDof[node1] == MpiInfo::rankProc()
1772 } // end of loop on iComp
1773 } // end of loop on iSupport
1774 ✗ ielGeom1++;
1775 } // end of loop on iel1 of each label
1776 } // end of if searchLabel == currentLabel
1777 else
1778 ✗ ielGeom1 += numEltPerLabel;
1779 } //end of loop on itRef
1780 } //end of loop on bagElementTypeDomainBoundary
1781 } // end of if lumpedModelBCAlgo[iLabel] == 2
1782 } //end of loop on lumpedModelBCLabel
1783
1784 // Storage the added dof in CSR format
1785 ✗ felInt dofSize = 0;
1786 ✗ for (unsigned int iNode = 0; iNode < nodesNeighborhood.size(); iNode++) {
1787 ✗ dofSize += nodesNeighborhood[iNode].size();
1788 }
1789
1790 ✗ std::vector<felInt> m_iCSR1(dof().pattern().numRows()+1, 0);
1791 ✗ std::vector<felInt> m_jCSR1(dofSize, 0);
1792 ✗ felInt pos = 0;
1793 ✗ felInt cptDof = 0;
1794 ✗ for (unsigned int iNode = 0; iNode < nodesNeighborhood.size(); iNode++) {
1795 ✗ if (rankDof[iNode] == MpiInfo::rankProc()) {
1796 ✗ m_iCSR1[cptDof+1 ] = m_iCSR1[cptDof] + nodesNeighborhood[iNode].size();
1797 ✗ pos = 0;
1798 ✗ for(auto iCon = nodesNeighborhood[iNode].begin(); iCon != nodesNeighborhood[iNode].end(); iCon++) {
1799 ✗ m_jCSR1[m_iCSR1[cptDof] + pos] = *iCon;
1800 ✗ pos++;
1801 }
1802 ✗ cptDof++;
1803 }
1804 }
1805
1806 // merge pattern
1807 ✗ dof().mergeGlobalPattern(m_iCSR1, m_jCSR1);
1808
1809 ✗ delete [] numDofPerProc;
1810 ✗ delete [] rankDof;
1811 ✗ nodesNeighborhood.clear();
1812 ✗ m_iCSR1.clear();
1813 ✗ m_jCSR1.clear();
1814 }
1815
1816 // ---------------------------------------------- //
1817 // Modify pattern for face-oriented stabilization //
1818 // ---------------------------------------------- //
1819
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
 54 if (r_instance.NSStabType == 1) {
1820 // If the face-oriented stabilization is used, the dofs of an element an all its neighbor
1821 // are linked together.
1822
1823 // compute initial (uniform) repartition of dof on processes
1824 ✗ felInt numDofByProc = m_numDof/MpiInfo::numProc();
1825 ✗ std::vector<felInt> numDofPerProc(MpiInfo::numProc());
1826 ✗ for(felInt i=0; i<MpiInfo::numProc(); ++i) {
1827 ✗ if(i == MpiInfo::numProc() - 1)
1828 ✗ numDofPerProc[i] = m_numDof - i*numDofByProc;
1829 else
1830 ✗ numDofPerProc[i] = numDofByProc;
1831 }
1832
1833 ✗ felInt shiftDof = 0;
1834 ✗ for(felInt i=0; i<MpiInfo::rankProc(); ++i)
1835 ✗ shiftDof += numDofPerProc[i];
1836
1837 ✗ std::vector<felInt> rankDof(m_numDof, -1);
1838 ✗ for(felInt i=0; i<numDofPerProc[MpiInfo::rankProc()]; ++i)
1839 ✗ rankDof[i + shiftDof] = MpiInfo::rankProc();
1840
1841 // variables
1842 GeometricMeshRegion::ElementType eltType, eltTypeOpp;
1843 ✗ const int idVar = m_listVariable.getVariableIdList(velocity);
1844 ✗ const std::size_t iMesh = m_listVariable[idVar].idMesh();
1845 felInt idVar1, idVar2;
1846 felInt node1, node2;
1847 felInt numFacesPerElement, numEltPerLabel;
1848 felInt ielSupport, jelSupport;
1849 felInt ielCurrentGeo, ielOppositeGeo;
1850 ✗ std::vector<felInt> vecSupport, vecSupportOpposite;
1851 ✗ std::vector<felInt> numElement(GeometricMeshRegion::m_numTypesOfElement, 0);
1852 ✗ std::vector< std::set<felInt> > nodesNeighborhood(m_numDof);
1853
1854 // zeroth loop on the unknown of the linear problem
1855 ✗ for (std::size_t iUnknown1 = 0; iUnknown1 < m_listUnknown.size(); ++iUnknown1) {
1856 ✗ idVar1 = m_listUnknown.idVariable(iUnknown1);
1857
1858 ✗ if ( iMesh != m_listVariable[idVar1].idMesh() )
1859 ✗ continue;
1860
1861 // build the edges or faces depending on the dimension (for the global mesh)
1862 // In this function, "faces" refers to either edges or faces.
1863 ✗ if(dimension() == 2)
1864 ✗ m_mesh[iMesh]->buildEdges();
1865 else
1866 ✗ m_mesh[iMesh]->buildFaces();
1867
1868 // first loop on element type
1869 ✗ ielCurrentGeo = 0;
1870 ✗ for (std::size_t i=0; i<m_mesh[iMesh]->bagElementTypeDomain().size(); ++i) {
1871 ✗ eltType = m_mesh[iMesh]->bagElementTypeDomain()[i];
1872 ✗ const GeoElement* geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
1873 ✗ numElement[eltType] = 0;
1874 ✗ numFacesPerElement = geoEle->numBdEle();
1875
1876 // second loop on region of the mesh.
1877 ✗ for(auto itRef = m_mesh[iMesh]->intRefToBegEndMaps[eltType].begin(); itRef != m_mesh[iMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
1878 ✗ numEltPerLabel = itRef->second.second;
1879
1880 // Third loop on element
1881 ✗ for (felInt iel = 0; iel < numEltPerLabel; iel++) {
1882 ✗ m_supportDofUnknown[iUnknown1].getIdElementSupport(eltType, numElement[eltType], vecSupport);
1883
1884 ✗ for(std::size_t ielSup=0; ielSup<vecSupport.size(); ++ielSup) {
1885 ✗ ielSupport = vecSupport[ielSup];
1886
1887 // for all support dof in this support element
1888 ✗ for (felInt iSup = 0; iSup < m_supportDofUnknown[iUnknown1].getNumSupportDof(ielSupport); ++iSup) {
1889 // for all component of the unknown
1890 ✗ for (std::size_t iComp = 0; iComp < m_listVariable[idVar1].numComponent(); iComp++) {
1891 // get the global id of the support dof
1892 ✗ dof().loc2glob(ielSupport, iSup, idVar1, iComp, node1);
1893
1894 // if this node is on this process
1895 ✗ if(rankDof[node1] == MpiInfo::rankProc()) {
1896 // loop over the boundaries of the element
1897 ✗ for(felInt iface=0; iface < numFacesPerElement; ++iface) {
1898 // check if this face is an inner face or a boundary
1899 ✗ ielOppositeGeo = ielCurrentGeo;
1900 ✗ eltTypeOpp = eltType;
1901 ✗ bool isAnInnerFace = m_mesh[iMesh]->getAdjElement(eltTypeOpp, ielOppositeGeo, iface);
1902
1903 ✗ if(isAnInnerFace) {
1904 // for all unknown
1905 ✗ for (std::size_t iUnknown2 = 0; iUnknown2 < m_listUnknown.size(); ++iUnknown2) {
1906 ✗ idVar2 = m_listUnknown.idVariable(iUnknown2);
1907
1908 ✗ if ( iMesh != m_listVariable[idVar2].idMesh() )
1909 ✗ continue;
1910
1911 // get the support element of the opposite element
1912 ✗ supportDofUnknown(iUnknown2).getIdElementSupport(ielOppositeGeo, vecSupportOpposite);
1913
1914 // for all support element
1915 ✗ for(std::size_t jelSup=0; jelSup<vecSupportOpposite.size(); ++jelSup) {
1916 ✗ jelSupport = vecSupportOpposite[jelSup];
1917
1918 // for all support dof in this support element
1919 ✗ for (felInt jSup = 0; jSup < supportDofUnknown(iUnknown2).getNumSupportDof(jelSupport); ++jSup) {
1920
1921 // for all component of the second unknown
1922 ✗ for (std::size_t jComp = 0; jComp < m_listVariable[idVar2].numComponent(); jComp++) {
1923
1924 // If the two current components are connected
1925 ✗ const felInt iConnect = dof().getNumGlobComp(iUnknown1, iComp);
1926 ✗ const felInt jConnect = dof().getNumGlobComp(iUnknown2, jComp);
1927 ✗ if (m_listUnknown.mask()(iConnect, jConnect) > 0) {
1928 // get the global id of the second support dof
1929 ✗ dof().loc2glob(jelSupport, jSup, idVar2, jComp, node2);
1930
1931 // remove diagonal term to define pattern and use it in Parmetis.
1932 // if the two global ids are different, they are neighboors
1933 ✗ if(node1 != node2)
1934 ✗ nodesNeighborhood[node1].insert(node2);
1935 }
1936 }
1937 }
1938 }
1939 }
1940 }
1941 }
1942 }
1943 }
1944 }
1945 }
1946 ✗ ++ielCurrentGeo;
1947 ✗ ++numElement[eltType];
1948 }
1949 }
1950 }
1951 }
1952
1953 // Store the pattern in CSR style
1954 ✗ std::vector<felInt> iCSR, jCSR;
1955 ✗ felInt dofSize = 0;
1956 ✗ felInt cptDof = 0;
1957 felInt pos;
1958
1959 ✗ iCSR.resize(dof().pattern().numRows() + 1, 0);
1960 ✗ for(std::size_t iNode=0; iNode<nodesNeighborhood.size(); ++iNode)
1961 ✗ dofSize += nodesNeighborhood[iNode].size();
1962
1963 ✗ jCSR.resize(dofSize, 0);
1964 ✗ for(std::size_t iNode=0; iNode<nodesNeighborhood.size(); ++iNode) {
1965 ✗ if(rankDof[iNode] == MpiInfo::rankProc()) {
1966 ✗ iCSR[cptDof + 1] = iCSR[cptDof] + nodesNeighborhood[iNode].size();
1967 ✗ pos = 0;
1968 ✗ for (auto it=nodesNeighborhood[iNode].begin(); it != nodesNeighborhood[iNode].end(); ++it) {
1969 ✗ jCSR[iCSR[cptDof] + pos] = *it;
1970 ✗ ++pos;
1971 }
1972 ✗ ++cptDof;
1973 }
1974 }
1975
1976 // Now, call merge pattern
1977 ✗ dof().mergeGlobalPattern(iCSR, jCSR);
1978 }
1979 54 }
1980
1981 /***********************************************************************************/
1982 /***********************************************************************************/
1983
1984 ✗ void LinearProblemNS::reviseSolution(PetscVector& alphai, std::vector<PetscVector>& stationnarySolutions)
1985 {
1986 // m_sol = u^n = \tilde{u}^n + \sum_i alpha_i^n u_i
1987
1988 // alpha_i^n = m_coefProblem->solution()[i]
1989 // u_i = m_stationnarySolutionsVec[i]
1990 // \tilde{u}^n = m_linearProblem[0]->solution() = m_sol already
1991
1992 PetscScalar coef ;
1993 ✗ for (PetscInt jj = 0 ; jj < (PetscInt) stationnarySolutions.size() ; jj++ ) {
1994 ✗ alphai.getValues(1, &jj, &coef);
1995 // m_sol = coef * statSol + m_sol
1996 ✗ solution().axpy(coef, stationnarySolutions[jj]);
1997 }
1998 }
1999
2000 /***********************************************************************************/
2001 /***********************************************************************************/
2002
2003 ✗ void LinearProblemNS::computeRHSThetaMethod(felInt& iel)
2004 {
2005 // Retrieve instance
2006 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
2007
2008 //==============================================
2009 // Provisional implementation
2010 //==============================================
2011
2012 ✗ m_elMatRHSThetaMethod[0]->mat().clear();
2013
2014 ✗ if ( ((m_fstransient->iteration == 1) && (r_instance.useALEformulation == 0))
2015 ✗ || ((m_fstransient->iteration > 0) && (r_instance.useALEformulation)) ) {
2016 std::size_t indexMatrix;
2017 ✗ if( r_instance.useALEformulation )
2018 ✗ indexMatrix = 0;
2019 else
2020 ✗ indexMatrix = 1;
2021 ✗ m_elMatRHSThetaMethod[indexMatrix]->mat().clear();
2022 ✗ if (m_useSymmetricStress) {
2023 ✗ m_elMatRHSThetaMethod[indexMatrix]->eps_phi_i_eps_phi_j((m_theta-1)*2*m_viscosity, *m_feVel, 0, 0, dimension());
2024 } else {
2025 ✗ m_elMatRHSThetaMethod[indexMatrix]->grad_phi_i_grad_phi_j((m_theta-1)*m_viscosity, *m_feVel, 0, 0, dimension());
2026 }
2027 }
2028
2029 ✗ if ( m_fstransient->iteration > 0) {
2030 //=======================
2031 // Convective Term
2032 //=======================
2033 ✗ switch(r_instance.NSequationFlag) {
2034 ✗ case 0:
2035 ✗ break;
2036 ✗ case 1:
2037 ✗ if(r_instance.useALEformulation) {
2038 ✗ m_elemFieldAdv.setValue(m_beta, *m_feVel, iel, m_iVelocity, m_ao, dof());
2039 } else {
2040 ✗ m_elemFieldAdv.setValue(m_solExtrapol, *m_feVel, iel, m_iVelocity, m_ao, dof());
2041 }
2042 ✗ m_elMatRHSThetaMethod[0]->u_grad_phi_j_phi_i((m_theta-1)*m_density,m_elemFieldAdv,*m_feVel,0,0,dimension());
2043 ✗ break;
2044 ✗ default:
2045 ✗ FEL_ERROR("Error: No theta-method for this kind of Convection Method ");
2046 }
2047 }
2048
2049 ✗ if( r_instance.useALEformulation == 0)
2050 ✗ m_elMatRHSThetaMethod[0]->mat() = m_elMatRHSThetaMethod[0]->mat() + m_elMatRHSThetaMethod[1]->mat();
2051
2052 ✗ m_elemFieldAdv.setValue(sequentialSolution(), *m_feVel, iel, m_iVelocity, m_ao, dof());
2053 ✗ m_elemFieldPres.setValue(sequentialSolution(), *m_fePres, iel, m_iPressure, m_ao, dof());
2054
2055 ✗ UBlasVector auxVec(m_elementVector[0]->vec().size());
2056 ✗ std::size_t entry=0;
2057
2058 ✗ for (std::size_t velComp=0; velComp <m_elemFieldAdv.get_val().size1() ; velComp++) {
2059 ✗ for (std::size_t velDof=0; velDof <m_elemFieldAdv.get_val().size2() ; velDof++) {
2060 ✗ auxVec(entry) = row(m_elemFieldAdv.get_val(),velComp)(velDof);
2061 ✗ entry++;
2062 }
2063 }
2064 ✗ for (std::size_t pressDof=0; pressDof <m_elemFieldPres.get_val().size2() ; pressDof++) {
2065 ✗ auxVec(entry) = row(m_elemFieldPres.get_val(),0)(pressDof);
2066 ✗ entry++;
2067 }
2068
2069 ✗ axpy_prod(m_elMatRHSThetaMethod[0]->mat(),auxVec,m_elementVector[0]->vec(), false);
2070
2071 }
2072
2073 /***********************************************************************************/
2074 /***********************************************************************************/
2075
2076 ✗ void LinearProblemNS::assembleFaceOrientedStabilization()
2077 {
2078 // Retrieve instance
2079 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
2080
2081 /*!
2082 * Assembly loop to compute the face-oriented stabilization
2083 * E. Burman, M. A. Fernandez and P. Hansbo, Continuous interior penalty finite element method
2084 * for Oseen's equations, SIAM J. Numer. Anal., 44(3), 1248-1274, 2006.
2085 *
2086 * This function works if the mesh has only one type of element !!!
2087 */
2088 // definition of variables
2089 ✗ std::pair<bool, GeometricMeshRegion::ElementType> adjElt;
2090 GeometricMeshRegion::ElementType eltType; // Type of element
2091
2092 ✗ felInt ielCurrentLocalGeo = 0; // local geometric id of the current element
2093 ✗ felInt ielCurrentGlobalGeo = 0; // global geometric id of the current element
2094 ✗ felInt ielOppositeGlobalGeo = 0; // global geometric id of the opposite element
2095 felInt numFacesPerType; // number of faces of an element of a given type
2096 felInt numPointPerElt; // number of vertices by element
2097 felInt ielCurrentGlobalSup; // global support element id of the current element
2098 felInt ielOppositeGlobalSup; // global support element id of the opposite element
2099
2100 ✗ std::vector<felInt> idOfFacesCurrent; // ids of the current element edges
2101 ✗ std::vector<felInt> idOfFacesOpposite; // ids of the opposite element edges
2102 ✗ std::vector<felInt> currentElemIdPoint; // ids of the vertices of the current element
2103 ✗ std::vector<felInt> oppositeElemIdPoint; // ids of the vertices of the opposite element
2104 ✗ std::vector<Point*> currentElemPoint; // point coordinates of the current element vertices
2105 ✗ std::vector<Point*> oppositeElemPoint; // point coordinates of the opposite element vertices
2106
2107 ✗ FlagMatrixRHS flag = FlagMatrixRHS::only_matrix; // flag to only assemble the matrix
2108 felInt idFaceToDo;
2109 ✗ bool allDone = false;
2110
2111 ✗ ElementField elemFieldAdvFace;
2112
2113 felInt rank;
2114 ✗ MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
2115
2116 ✗ const int iVel = m_listVariable.getVariableIdList(velocity); // TODO D.C. ALREADY DEFINED IN CALSS???
2117 ✗ const int iPres = m_listVariable.getVariableIdList(pressure);
2118 ✗ m_velocity = &m_listVariable[iVel ];
2119 ✗ m_pressure = &m_listVariable[iPres];
2120
2121 ✗ if ( m_velocity->idMesh() != m_pressure->idMesh() && m_pressure->idMesh() != static_cast<std::size_t>(m_currentMesh) )
2122 ✗ FEL_ERROR("Velocity and pressure defined on different meshes, or not defined on m_currentMesh");
2123
2124 // bag element type vector
2125 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
2126 ✗ eltType = bagElementTypeDomain[0];
2127
2128 // Finite Element With Bd for the opposite element
2129 CurrentFiniteElementWithBd* ptmp;
2130 CurrentFiniteElementWithBd* oppositeFEWithBdVel;
2131 CurrentFiniteElementWithBd* oppositeFEWithBdPre;
2132
2133 ✗ const GeoElement *geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
2134
2135 ✗ felInt typeOfFEVel = m_velocity->finiteElementType();
2136 ✗ const RefElement *refEleVel = geoEle->defineFiniteEle(eltType, typeOfFEVel, *m_mesh[m_currentMesh]);
2137 ✗ oppositeFEWithBdVel = new CurrentFiniteElementWithBd(*refEleVel, *geoEle, m_velocity->degreeOfExactness(), m_velocity->degreeOfExactness());
2138
2139 ✗ felInt typeOfFEPre = m_pressure->finiteElementType();
2140 ✗ const RefElement *refElePre = geoEle->defineFiniteEle(eltType, typeOfFEPre, *m_mesh[m_currentMesh]);
2141 ✗ oppositeFEWithBdPre = new CurrentFiniteElementWithBd(*refElePre, *geoEle, m_pressure->degreeOfExactness(), m_pressure->degreeOfExactness());
2142
2143
2144 // initializing variables
2145 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
2146 ✗ numFacesPerType = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numBdEle();
2147
2148 // resize of vectors
2149 ✗ idOfFacesCurrent.resize(numFacesPerType, -1);
2150 ✗ idOfFacesOpposite.resize(numFacesPerType, -1);
2151 ✗ currentElemIdPoint.resize(numPointPerElt, -1);
2152 ✗ oppositeElemIdPoint.resize(numPointPerElt, -1);
2153 ✗ currentElemPoint.resize(numPointPerElt, nullptr);
2154 ✗ oppositeElemPoint.resize(numPointPerElt, nullptr);
2155
2156 // define finite element
2157 ✗ defineFiniteElement(eltType);
2158 ✗ initElementArray();
2159 ✗ defineCurrentFiniteElementWithBd(eltType);
2160
2161 // allocate arrays for assembling the matrix
2162 ✗ allocateArrayForAssembleMatrixRHS(flag);
2163
2164 // init variables
2165 ✗ initPerElementType(eltType, flag);
2166 ✗ m_feVelWithBd = m_listCurrentFiniteElementWithBd[m_iVelocity];
2167 ✗ m_fePreWithBd = m_listCurrentFiniteElementWithBd[m_iPressure];
2168
2169
2170 ✗ CurrentFiniteElementWithBd* firstCurrentVel = m_feVelWithBd;
2171 ✗ CurrentFiniteElementWithBd* firstCurrentPre = m_fePreWithBd;
2172 ✗ CurrentFiniteElementWithBd* firstOppositeVel = oppositeFEWithBdVel;
2173 ✗ CurrentFiniteElementWithBd* firstOppositePre = oppositeFEWithBdPre;
2174
2175 // get informations on the current element
2176 ✗ setElemPoint(eltType, 0, currentElemPoint, currentElemIdPoint, &ielCurrentGlobalSup);
2177
2178 // update the finite elements
2179 ✗ m_feVelWithBd->updateFirstDeriv(0, currentElemPoint);
2180 ✗ m_feVelWithBd->updateBdMeasNormal();
2181 ✗ m_fePreWithBd->updateFirstDeriv(0, currentElemPoint);
2182 ✗ m_fePreWithBd->updateBdMeasNormal();
2183
2184 // get the global id of the first geometric element
2185 ✗ ISLocalToGlobalMappingApply(m_mappingElem[m_currentMesh], 1, &ielCurrentLocalGeo, &ielCurrentGlobalGeo);
2186
2187 // the std::unordered_map to remember what contribution have already been computed
2188 ✗ std::map<felInt, std::vector<bool> > contribDone;
2189 ✗ addNewFaceOrientedContributor(numFacesPerType, ielCurrentGlobalGeo, contribDone[ielCurrentGlobalGeo]);
2190
2191 // build the stack to know what is the next element
2192 ✗ std::stack<felInt> nextInStack;
2193 ✗ nextInStack.push(ielCurrentGlobalGeo);
2194
2195 // get all the faces of the element
2196 ✗ if(dimension() == 2)
2197 ✗ m_mesh[m_currentMesh]->getAllEdgeOfElement(ielCurrentGlobalGeo, idOfFacesCurrent);
2198 else
2199 ✗ m_mesh[m_currentMesh]->getAllFaceOfElement(ielCurrentGlobalGeo, idOfFacesCurrent);
2200
2201 // loop over all the element (snake style)
2202 ✗ while(!nextInStack.empty()) {
2203 // check the faces and use the first one that is an inner face and that has not been done yet.
2204 ✗ idFaceToDo = -1;
2205 ✗ allDone = true;
2206 ✗ for(std::size_t iface=0; iface<contribDone[ielCurrentGlobalGeo].size(); ++iface) {
2207 ✗ if(!contribDone[ielCurrentGlobalGeo][iface]) {
2208 // This is an inner face, check if the contribution has already been computed or not
2209 ✗ if(idFaceToDo == -1) {
2210 ✗ idFaceToDo = iface;
2211 } else {
2212 ✗ allDone = false;
2213 }
2214 }
2215 }
2216
2217 // update the stack
2218 ✗ if(!allDone)
2219 ✗ nextInStack.push(ielCurrentGlobalGeo);
2220
2221 ✗ if(nextInStack.top() == ielCurrentGlobalGeo && allDone)
2222 ✗ nextInStack.pop();
2223
2224
2225 // assemble terms
2226 ✗ if(idFaceToDo != -1) {
2227 // get the opposite id of the element
2228 ✗ ielOppositeGlobalGeo = ielCurrentGlobalGeo;
2229 ✗ adjElt = m_mesh[m_currentMesh]->getAdjElement(ielOppositeGlobalGeo, idFaceToDo);
2230
2231 // get the type of the opposite element
2232 // eltTypeOpp = adjElt.second;
2233
2234 // update the opposite finite element and std::set ielOppositeGlobalSup
2235 ✗ updateFaceOrientedFEWithBd(oppositeFEWithBdVel, oppositeFEWithBdPre, idOfFacesOpposite, numPointPerElt, ielOppositeGlobalGeo, ielOppositeGlobalSup);
2236
2237 // find the local id of the edge in the opposite element
2238 ✗ felInt localIdFaceOpposite = -1;
2239 ✗ for(std::size_t jface=0; jface<idOfFacesOpposite.size(); ++jface) {
2240 ✗ if(idOfFacesCurrent[idFaceToDo] == idOfFacesOpposite[jface]) {
2241 ✗ localIdFaceOpposite = jface;
2242 }
2243 }
2244
2245
2246 // compute coefficients
2247 ✗ elemFieldAdvFace.initialize(DOF_FIELD, m_feVelWithBd->bdEle(idFaceToDo), dimension());
2248 ✗ elemFieldAdvFace.setValue(m_solExtrapol, *m_feVelWithBd, idFaceToDo, ielCurrentGlobalSup, m_iVelocity, m_ao, dof());
2249 ✗ UBlasVector tmpval = prod(trans(elemFieldAdvFace.val), m_feVelWithBd->bdEle(idFaceToDo).normal[0]);
2250 ✗ double normVelNormalBd = norm_inf(tmpval);
2251 ✗ double normVelBd = 0.;
2252 ✗ for(std::size_t jcol=0; jcol<elemFieldAdvFace.val.size2(); ++jcol)
2253 ✗ normVelBd = std::max(normVelBd, norm_2(column(elemFieldAdvFace.val, jcol)));
2254
2255 ✗ double hK = 0.5 * (m_feVelWithBd->diameter() + oppositeFEWithBdVel->diameter());
2256 ✗ double ReK = m_density * normVelBd * hK / m_viscosity;
2257
2258 ✗ double coeffVel = r_instance.stabFOVel * std::min(1., ReK) * hK * hK;
2259 ✗ double gammaU = coeffVel * normVelNormalBd;
2260 ✗ double gammaP = 0;
2261 ✗ if(r_instance.NSequationFlag == 0) {
2262 ✗ gammaP = r_instance.stabFOPre * hK * hK * hK / m_viscosity;
2263 } else {
2264 ✗ double xiPre = ReK < 1 ? m_density * hK * hK * hK / m_viscosity : hK * hK / normVelBd;
2265 ✗ gammaP = r_instance.stabFOPre * xiPre;
2266 }
2267
2268 // We can now compute all the integrals.
2269 // There is a minus sign for gammaP because there is one for q div(u).
2270 // current element for phi_i and phi_j
2271 ✗ m_elementMat[0]->zero();
2272 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammaU, *m_feVelWithBd, *m_feVelWithBd, idFaceToDo, idFaceToDo, 0, 0, dimension());
2273 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(-gammaP, *m_fePreWithBd, *m_fePreWithBd, idFaceToDo, idFaceToDo, dimension(), dimension(), 1);
2274 ✗ setValueMatrixRHS(ielCurrentGlobalSup, ielCurrentGlobalSup, flag);
2275
2276 // current element for phi_i and opposite element for phi_j
2277 ✗ m_elementMat[0]->zero();
2278 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammaU, *oppositeFEWithBdVel, *m_feVelWithBd, localIdFaceOpposite, idFaceToDo, 0, 0, dimension());
2279 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(-gammaP, *oppositeFEWithBdPre, *m_fePreWithBd, localIdFaceOpposite, idFaceToDo, dimension(), dimension(), 1);
2280 ✗ setValueMatrixRHS(ielCurrentGlobalSup, ielOppositeGlobalSup, flag);
2281
2282 // opposite element for phi_i and current element for phi_j
2283 ✗ m_elementMat[0]->mat() = trans(m_elementMat[0]->mat());
2284 ✗ setValueMatrixRHS(ielOppositeGlobalSup, ielCurrentGlobalSup, flag);
2285
2286 // opposite element for phi_i and phi_j
2287 ✗ m_elementMat[0]->zero();
2288 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammaU, *oppositeFEWithBdVel, *oppositeFEWithBdVel, localIdFaceOpposite, localIdFaceOpposite, 0, 0, dimension());
2289 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(-gammaP, *oppositeFEWithBdPre, *oppositeFEWithBdPre, localIdFaceOpposite, localIdFaceOpposite, dimension(), dimension(), 1);
2290 ✗ setValueMatrixRHS(ielOppositeGlobalSup, ielOppositeGlobalSup, flag);
2291
2292 // update the std::unordered_map to say that this face is done
2293 ✗ contribDone[ielCurrentGlobalGeo][idFaceToDo] = true;
2294
2295 // check if the opposite element is on this proc
2296 ✗ if(m_eltPart[m_currentMesh][ielOppositeGlobalGeo] == rank) {
2297 ✗ if(contribDone.find(ielOppositeGlobalGeo) == contribDone.end()) {
2298 // not found in the std::unordered_map, add it and initialize it.
2299 ✗ addNewFaceOrientedContributor(numFacesPerType, ielOppositeGlobalGeo, contribDone[ielOppositeGlobalGeo]);
2300 }
2301 ✗ contribDone[ielOppositeGlobalGeo][localIdFaceOpposite] = true;
2302
2303 // update the finite element as the opposite finite element.
2304 ✗ ptmp = m_feVelWithBd;
2305 ✗ m_feVelWithBd = oppositeFEWithBdVel;
2306 ✗ oppositeFEWithBdVel = ptmp;
2307
2308 ✗ ptmp = m_fePreWithBd;
2309 ✗ m_fePreWithBd = oppositeFEWithBdPre;
2310 ✗ oppositeFEWithBdPre = ptmp;
2311
2312 ✗ ielCurrentGlobalGeo = ielOppositeGlobalGeo;
2313 ✗ ielCurrentGlobalSup = ielOppositeGlobalSup;
2314 ✗ idOfFacesCurrent = idOfFacesOpposite;
2315 } else {
2316 // not on this rank, take the next one in the stack
2317 ✗ if(!nextInStack.empty()) {
2318 ✗ ielCurrentGlobalGeo = nextInStack.top();
2319 ✗ updateFaceOrientedFEWithBd(m_feVelWithBd, m_fePreWithBd, idOfFacesCurrent, numPointPerElt, ielCurrentGlobalGeo, ielCurrentGlobalSup);
2320 }
2321 }
2322 ✗ } else {
2323 // All contribution have been already computed on this element
2324 // take the next element in the stack
2325 ✗ if(!nextInStack.empty()) {
2326 ✗ ielCurrentGlobalGeo = nextInStack.top();
2327 ✗ updateFaceOrientedFEWithBd(m_feVelWithBd, m_fePreWithBd, idOfFacesCurrent, numPointPerElt, ielCurrentGlobalGeo, ielCurrentGlobalSup);
2328 }
2329 }
2330 }
2331
2332 // desallocate array use for assemble with petsc.
2333 ✗ desallocateArrayForAssembleMatrixRHS(flag);
2334
2335 // desallocate opposite finite elements
2336 ✗ m_feVelWithBd = firstCurrentVel;
2337 ✗ m_fePreWithBd = firstCurrentPre;
2338 ✗ delete firstOppositeVel;
2339 ✗ delete firstOppositePre;
2340 }
2341
2342 ✗ void LinearProblemNS::addNewFaceOrientedContributor(felInt size, felInt idElt, std::vector<bool>& vec)
2343 {
2344 felInt ielTmp;
2345 ✗ std::pair<bool, GeometricMeshRegion::ElementType> adjElt;
2346
2347 ✗ vec.resize(size, false);
2348 ✗ for(felInt iface=0; iface<size; ++iface) {
2349 // check if this face is an inner face or a boundary
2350 ✗ ielTmp = idElt;
2351 ✗ adjElt = m_mesh[m_currentMesh]->getAdjElement(ielTmp, iface);
2352
2353 ✗ if(!adjElt.first) {
2354 // not an inner face, set the contribution to done.
2355 ✗ vec[iface] = true;
2356 }
2357 }
2358 }
2359
2360 ✗ void LinearProblemNS::updateFaceOrientedFEWithBd(CurrentFiniteElementWithBd* fevel, CurrentFiniteElementWithBd* fepre, std::vector<felInt>& idOfFaces, felInt numPoints, felInt idElt, felInt& idSup)
2361 {
2362 ✗ std::vector<felInt> elemIdPoint(numPoints, -1);
2363 ✗ std::vector<Point*> elemPoint(numPoints, nullptr);
2364
2365 // get information on the opposite element
2366 // same as the function "setElemPoint" but here ielOppositeGlobalGeo is global
2367 // we assume that the supportDofMesh is the same for all unknown (like in setElemPoint)
2368 ✗ m_supportDofUnknown[0].getIdElementSupport(idElt, idSup);
2369 ✗ m_mesh[m_currentMesh]->getOneElement(idElt, elemIdPoint);
2370 ✗ for (felInt iPoint=0; iPoint<numPoints; ++iPoint)
2371 ✗ elemPoint[iPoint] = &(m_mesh[m_currentMesh]->listPoints()[elemIdPoint[iPoint]]);
2372
2373 // update the finite elements
2374 ✗ fevel->updateFirstDeriv(0, elemPoint);
2375 ✗ fevel->updateBdMeasNormal();
2376 ✗ fepre->updateFirstDeriv(0, elemPoint);
2377 ✗ fepre->updateBdMeasNormal();
2378
2379 // get all the faces of the element
2380 ✗ if(dimension() == 2)
2381 ✗ m_mesh[m_currentMesh]->getAllEdgeOfElement(idElt, idOfFaces);
2382 else
2383 ✗ m_mesh[m_currentMesh]->getAllFaceOfElement(idElt, idOfFaces);
2384 }
2385
2386 /***********************************************************************************/
2387 /***********************************************************************************/
2388
2389 12 void LinearProblemNS::computeNormalField(const std::vector<int> & labels, felInt idVecVariable)
2390 {
2391 12 m_auxiliaryInts.resize(1);
2392 12 m_auxiliaryInts[0]=idVecVariable;//for instance m_iVelocity;
2393 /*!
2394 For each dof d of the surface the sum: \f$\sum_{K} m(K)\vec n(K)\f$ is computed.
2395 - K are all the boundary elements that share the dof d.
2396 - m(K) is the measure of the element
2397 - n(K) is the P0 normal defined on the element
2398 */
2399 12 assemblyLoopBoundaryGeneral(&LinearProblemNS::normalFieldComputer,labels,&LinearProblemNS::initPerETNormVel,&LinearProblemNS::updateFeNormVel);
2400
2401
3/6
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
 12 m_vecs.Get("normalField").assembly();
2402
3/6
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
 12 m_vecs.Get("measStar").assembly();
2403
2404 /// Than the vector is divided by \f$\sum_{K} m(K)\f$ to obtain the P1 field.
2405
7/14
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 12 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 12 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 12 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 12 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 12 times.
✗ Branch 23 not taken.
 12 pointwiseDivide(m_vecs.Get("normalField"), m_vecs.Get("normalField"), m_vecs.Get("measStar"));
2406
5/10
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 12 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 12 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 12 times.
✗ Branch 16 not taken.
 12 gatherVector( m_vecs.Get("normalField"), m_seqVecs.Get("normalField") );
2407
2408 /// using normalizeComputer we normalize the result.
2409 12 assemblyLoopBoundaryGeneral(&LinearProblemNS::normalizeComputer,labels,&LinearProblemNS::initPerETNormVel,&LinearProblemNS::updateFeNormVel);
2410
2411
3/6
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
 12 m_vecs.Get("normalField").assembly();
2412
5/10
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 12 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 12 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 12 times.
✗ Branch 16 not taken.
 12 gatherVector( m_vecs.Get("normalField"), m_seqVecs.Get("normalField") );
2413 12 }
2414
2415 /***********************************************************************************/
2416 /***********************************************************************************/
2417
2418 22499 void LinearProblemNS::normalFieldComputer(felInt ielSupportDof)
2419 {
2420 double value;
2421
2/2
✓ Branch 3 taken 67417 times.
✓ Branch 4 taken 22499 times.
 89916 for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numDof(); iSurfDof++) {
2422 67417 double feMeas(m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->measure());
2423
2/2
✓ Branch 3 taken 202091 times.
✓ Branch 4 taken 67417 times.
 269508 for ( std::size_t iCoor(0); iCoor < ( std::size_t ) m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numCoor(); iCoor++) {
2424 felInt iDofGlob;
2425
1/2
✓ Branch 3 taken 202091 times.
✗ Branch 4 not taken.
 202091 dof().loc2glob( ielSupportDof, iSurfDof, m_auxiliaryInts[0], iCoor, iDofGlob);
2426
1/2
✓ Branch 1 taken 202091 times.
✗ Branch 2 not taken.
 202091 AOApplicationToPetsc(m_ao,1,&iDofGlob);
2427 202091 value=0.0;
2428
2/2
✓ Branch 3 taken 1211586 times.
✓ Branch 4 taken 202091 times.
 1413677 for ( std::size_t iQuadBd(0.0) ; iQuadBd < ( std::size_t )m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numQuadraturePoint(); iQuadBd++) {
2429
2/4
✓ Branch 3 taken 1211586 times.
✗ Branch 4 not taken.
✓ Branch 9 taken 1211586 times.
✗ Branch 10 not taken.
 1211586 value += m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->weightMeas( iQuadBd ) * m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->normal[iQuadBd](iCoor);
2430 }
2431
3/6
✓ Branch 2 taken 202091 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 202091 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 202091 times.
✗ Branch 9 not taken.
 202091 m_vecs.Get("normalField").setValue(iDofGlob, value, ADD_VALUES);
2432
3/6
✓ Branch 2 taken 202091 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 202091 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 202091 times.
✗ Branch 9 not taken.
 202091 m_vecs.Get("measStar").setValue(iDofGlob, feMeas, ADD_VALUES);
2433 }
2434 }
2435 22499 }
2436
2437 /***********************************************************************************/
2438 /***********************************************************************************/
2439
2440 44998 void LinearProblemNS::updateFeNormVel(const std::vector<Point*>& elemPoint,const std::vector<int>& /*elemIdPoint*/)
2441 {
2442 44998 m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->updateMeasNormal(0,elemPoint);
2443 44998 }
2444
2445 /***********************************************************************************/
2446 /***********************************************************************************/
2447
2448 22499 void LinearProblemNS::normalizeComputer(felInt ielSupportDof)
2449 {
2450 22499 const int numComp = m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numCoor( );
2451 22499 const int numDof = m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numDof( );
2452
2453
1/2
✓ Branch 2 taken 22499 times.
✗ Branch 3 not taken.
 22499 std::vector< double > norm(numDof,0.0);
2454
2/4
✓ Branch 3 taken 22499 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 22499 times.
✗ Branch 7 not taken.
 44998 std::vector< std::vector< double > > NN(numDof, std::vector< double >(numComp, 0.0));
2455
2/2
✓ Branch 0 taken 67417 times.
✓ Branch 1 taken 22499 times.
 89916 for ( int idof = 0; idof < numDof; idof++ ) {
2456
2/2
✓ Branch 0 taken 202091 times.
✓ Branch 1 taken 67417 times.
 269508 for ( int icomp = 0; icomp < numComp; icomp++ ) {
2457 int iDofGlob;
2458 double value;
2459
1/2
✓ Branch 3 taken 202091 times.
✗ Branch 4 not taken.
 202091 dof().loc2glob( ielSupportDof, idof, m_auxiliaryInts[0], icomp, iDofGlob);
2460
1/2
✓ Branch 1 taken 202091 times.
✗ Branch 2 not taken.
 202091 AOApplicationToPetsc(m_ao, 1, &iDofGlob);
2461
3/6
✓ Branch 2 taken 202091 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 202091 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 202091 times.
✗ Branch 9 not taken.
 202091 m_seqVecs.Get("normalField").getValues(1, &iDofGlob, &value);
2462 202091 norm[idof] += value*value;
2463 202091 NN[idof][icomp] = value;
2464 }
2465 67417 norm[idof] = std::sqrt( norm[idof] );
2466
2/2
✓ Branch 0 taken 202091 times.
✓ Branch 1 taken 67417 times.
 269508 for ( int icomp = 0; icomp < numComp; icomp++ )
2467 202091 NN[idof][icomp] = NN[idof][icomp]/norm[idof];
2468 }
2469
2/2
✓ Branch 0 taken 67417 times.
✓ Branch 1 taken 22499 times.
 89916 for ( int idof = 0; idof < numDof; idof++ ) {
2470
2/2
✓ Branch 0 taken 202091 times.
✓ Branch 1 taken 67417 times.
 269508 for ( int icomp = 0; icomp < numComp; icomp++ ) {
2471 int iDofGlob;
2472
1/2
✓ Branch 3 taken 202091 times.
✗ Branch 4 not taken.
 202091 dof().loc2glob( ielSupportDof, idof, m_auxiliaryInts[0], icomp, iDofGlob);
2473
1/2
✓ Branch 1 taken 202091 times.
✗ Branch 2 not taken.
 202091 AOApplicationToPetsc(m_ao, 1, &iDofGlob);
2474
3/6
✓ Branch 2 taken 202091 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 202091 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 202091 times.
✗ Branch 11 not taken.
 202091 m_vecs.Get("normalField").setValue(iDofGlob, NN[idof][icomp], INSERT_VALUES);
2475 }
2476 }
2477 22499 }
2478
2479 /***********************************************************************************/
2480 /***********************************************************************************/
2481
2482 #ifdef FELISCE_WITH_CVGRAPH
2483 ✗ void LinearProblemNS::assembleMassBoundaryAndInitKSP( std::size_t iConn )
2484 {
2485 ✗ LinearProblem::assembleMassBoundaryAndInitKSP(&LinearProblemNS::massMatrixComputer,
2486 ✗ slave()->interfaceLabels(iConn),
2487 &LinearProblemNS::initPerETMass,
2488 &LinearProblemNS::updateFE,
2489 iConn);
2490 }
2491
2492 ✗ void LinearProblemNS::sendData()
2493 {
2494 ✗ gatherSolution();
2495 ✗ for ( std::size_t iConn(0); iConn<slave()->numConnections(); ++iConn ) {
2496 ✗ std::vector<PetscVector> vecs;
2497 ✗ for (std::size_t cVar(0); cVar<slave()->numVarToSend(iConn); ++cVar) {
2498 ✗ std::string varToSend=slave()->sendVariable(iConn,cVar);
2499 ✗ if ( varToSend== "STRESS" ) {
2500 ✗ prepareResidual(iConn);
2501 ✗ vecs.push_back(m_seqStressOut);
2502 ✗ } else if ( varToSend == "VELOCITY" ) {
2503 ✗ vecs.push_back(sequentialSolution());
2504 ✗ } else if ( varToSend == "DISPLACEMENT" ) {
2505 //TODO curDisp = oldDisp+dt*sequentialSolution();
2506 //vecs.push_back(curDisp);
2507 ✗ FEL_ERROR("TODO: implement it!");
2508 } else {
2509 ✗ std::cout<<"variable to send: "<<varToSend<<std::endl;
2510 ✗ FEL_ERROR("Error in variable to read");
2511 }
2512 }
2513 ✗ slave()->sendData(vecs,iConn);
2514 }
2515 }
2516
2517 ✗ void LinearProblemNS::readData()
2518 {
2519 ✗ bool usingDisp(false);
2520 ✗ bool usingVel(false);
2521 ✗ for ( std::size_t iConn(0); iConn<slave()->numConnections(); ++iConn ) {
2522 ✗ std::vector<PetscVector> vecs;
2523 ✗ for (std::size_t cVar(0); cVar<slave()->numVarToRead(iConn); ++cVar) {
2524 ✗ std::string varToRead=slave()->readVariable(iConn,cVar);
2525 ✗ if ( varToRead == "STRESS" ) {
2526 ✗ vecs.push_back(m_seqVecs.Get("cvgraphSTRESS"));
2527 ✗ } else if ( varToRead == "DISPLACEMENT" ) {
2528 ✗ usingDisp=true;
2529 ✗ vecs.push_back(m_seqVecs.Get("externalDisplacement"));
2530 ✗ } else if ( varToRead == "VELOCITY" ) {
2531 ✗ usingVel=true;
2532 ✗ vecs.push_back(m_seqVecs.Get("cvgraphVelocity"));
2533 } else {
2534 ✗ std::cout<<"Connection: "<<iConn<<"variable to read: "<<varToRead<<std::endl;
2535 ✗ FEL_ERROR("Error in variable to read");
2536 }
2537 }
2538 ✗ slave()->receiveData(vecs,iConn);
2539 }
2540 ✗ if ( usingDisp && usingVel ) {
2541 // It is not safe to mix velocity vars and displacement vars, even from different connections.
2542 // It is possible to do it of course, but the implementation should be checked carefully.
2543 ✗ FEL_ERROR("You can not mix velocity and displacement.");
2544 }
2545 ✗ if ( usingDisp ) {
2546 ✗ m_seqVecs.Get("cvgraphVelocity").axpbypcz(1./m_fstransient->timeStep, -1./m_fstransient->timeStep, 0, m_seqVecs.Get("externalDisplacement"), m_seqVecs.Get("externalDisplacementOld"));
2547 }
2548 }
2549
2550 /*! \brief Function to assemble the mass matrix
2551 * Function to be called in the assemblyLoopBoundaryGeneral
2552 */
2553 ✗ void LinearProblemNS::initPerETMass()
2554 {
2555 ✗ felInt numDofTotal = 0;
2556 //pay attention this numDofTotal is bigger than expected since it counts for all the VARIABLES
2557 ✗ for (std::size_t iFe = 0; iFe < m_listCurvilinearFiniteElement.size(); iFe++) {//this loop is on variables while it should be on unknown
2558 ✗ numDofTotal += m_listCurvilinearFiniteElement[iFe]->numDof()*m_listVariable[iFe].numComponent();
2559 }
2560 ✗ m_globPosColumn.resize(numDofTotal); m_globPosRow.resize(numDofTotal); m_matrixValues.resize(numDofTotal*numDofTotal);
2561 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[ m_iVelocity ];
2562 }
2563
2564 ✗ void LinearProblemNS::massMatrixComputer(felInt ielSupportDof)
2565 {
2566 ✗ m_elementMatBD[0]->zero();
2567 ✗ m_elementMatBD[0]->phi_i_phi_j(/*coef*/1.,*m_curvFeVel,/*iblock*/0,/*iblock*/0,dimension());
2568 ✗ setValueMatrixBD(ielSupportDof);
2569 }
2570
2571 ✗ void LinearProblemNS::updateFE(const std::vector<Point*>& elemPoint, const std::vector<int>&)
2572 {
2573 ✗ m_curvFeVel->updateMeasNormal(0, elemPoint);
2574 }
2575
2576 ✗ void LinearProblemNS::computeBoundaryStress(PetscVector& vec)
2577 {
2578 // Retrieve instance
2579 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
2580
2581 // This function replaces computeBoundaryStress implemented in linearProblem.cpp
2582 // I think this function is safer and more general (should work in parallel as well)
2583
2584 ✗ ElementField velocityDofValuesVol;
2585 ✗ ElementField pressureDofValuesBD;
2586 ✗ ElementField tractionVel; // velocity part of \sigma \vect n
2587 ✗ ElementField tractionPre; // pressure part of \sigma \vect n: -p\vect n
2588
2589 /*===========================================*/
2590 // initializations for the loop
2591 ✗ ElementType eltType,volEltType,previousVolEltType=GeometricMeshRegion::Nod;
2592 ✗ int numPointPerElt = 0;
2593 ✗ int currentLabel = 0;
2594 ✗ felInt numEltPerLabel = 0;
2595 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
2596 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
2597 ✗ eltType = (ElementType)ityp;
2598 ✗ numElement[eltType] = 0;
2599 }
2600 ✗ if(m_mesh[m_currentMesh]->statusFaces() == false) {
2601 ✗ m_mesh[m_currentMesh]->buildFaces();
2602 }
2603 ✗ std::vector<Point*> elemPoint;
2604 ✗ std::vector<felInt> elemIdPoint;
2605 ✗ std::vector<felInt> vectorIdSupport;
2606 felInt ielSupportDof;
2607 ✗ allocateVectorBoundaryConditionDerivedLinPb();
2608 ✗ const std::vector<ElementType>& bagElementTypeDomainBoundary = m_meshLocal[m_currentMesh]->bagElementTypeDomainBoundary();
2609 /*===========================================*/
2610
2611 ✗ for (std::size_t i=0; i < bagElementTypeDomainBoundary.size(); ++i) {// loop over the boundary element type
2612 ✗ eltType = bagElementTypeDomainBoundary[i];
2613 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
2614 ✗ elemPoint.resize(numPointPerElt, nullptr);
2615 ✗ elemIdPoint.resize(numPointPerElt, 0);
2616 ✗ defineCurvilinearFiniteElement(eltType); // m_listCurvilinearFiniteElement update
2617 ✗ initElementArrayBD();
2618 ✗ allocateArrayForAssembleMatrixRHSBD(FlagMatrixRHS::only_rhs); // works in combo with setValueCustomVectorBD
2619
2620 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
2621 ✗ m_curvFePre = m_listCurvilinearFiniteElement[m_iPressure];
2622 // dof values
2623 ✗ pressureDofValuesBD.initialize(DOF_FIELD, *m_curvFePre,1/*numComp*/);
2624 // for the results
2625 ✗ tractionVel.initialize(DOF_FIELD, *m_curvFeVel, dimension());
2626 ✗ tractionPre.initialize(DOF_FIELD, *m_curvFePre, dimension());
2627
2628 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin();
2629 ✗ itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {// loop over the boundary labels that have this ref
2630 ✗ currentLabel = itRef->first;
2631 ✗ numEltPerLabel = itRef->second.second;
2632 ✗ initPerDomainBoundaryCondition(elemPoint, elemIdPoint, currentLabel, numEltPerLabel, &ielSupportDof, eltType, numElement[eltType], FlagMatrixRHS::only_rhs); //it updates label
2633 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {// loop over the elements of this label
2634 // for each element:
2635 ✗ m_elementVectorBD[0]->zero();
2636 // - get the points and their ids and vectorIdSupport
2637 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, vectorIdSupport);
2638 ✗ ielSupportDof=vectorIdSupport[0];// We assume we are not working with duplicated support dof: TODO
2639 // - get the volume points and their ids
2640 ✗ std::vector<Point*> elemPointVol; // Points of the current domain element.
2641 ✗ std::vector<felInt> elemIdPointVol; // Id of points of the domain element.
2642 ✗ std::vector<felInt> bd2vol;
2643 ✗ int idLocFace = -1; // Local identity of the face of a domain element
2644 ✗ felInt idElemVol = -1; // ID of the domain element associated to the face with identity idLocFace
2645 ✗ m_mesh[m_currentMesh]->getElementFromBdElem(elemIdPoint,elemIdPointVol,elemPointVol,idLocFace,idElemVol);
2646 int dummy;
2647 ✗ m_mesh[m_currentMesh]->getTypeElemFromIdElem(idElemVol, volEltType, dummy);
2648 ✗ if (volEltType != previousVolEltType ) {
2649 ✗ defineFiniteElement(volEltType);
2650 ✗ initElementArray();
2651 ✗ previousVolEltType=volEltType;
2652 ✗ m_feVel = m_listCurrentFiniteElement[m_iVelocity];
2653 ✗ velocityDofValuesVol.initialize(DOF_FIELD, *m_feVel,dimension());
2654 ✗ m_feVel->allocateOnDofsDataStructures();
2655 }
2656 ✗ m_feVel->updateFirstDerivOnDofs(0, elemPointVol);
2657 ✗ m_curvFeVel->updateMeasNormal(0,elemPoint);
2658 ✗ m_curvFePre->updateMeasNormal(0,elemPoint);
2659 ✗ m_feVel->identifyLocBDDof(*m_curvFeVel,bd2vol);
2660
2661 // Extract the corresponding values of both pressure and velocity from the sequentialSolution.
2662 ✗ velocityDofValuesVol.setValue(sequentialSolution(),*m_feVel,idElemVol,m_iVelocity,m_ao,dof());
2663 ✗ pressureDofValuesBD.setValue(sequentialSolution(),*m_curvFePre,ielSupportDof,m_iPressure,m_ao,dof());
2664
2665 // Create an element field (dof) to store the traction \vect t = \sigma \vect n
2666 // compute \vect t = -p \vect n
2667
2668 //Attention: we have assumed the normal is constant on the element, and we have taken it on the first quadrature node.
2669 // if we want to generalize this we have to transform tractionPres into an element field defined on quadrature nodes
2670 // and we have to evalu
2671 ✗ tractionPre.val = -outer_prod(m_curvFePre->normal[0],pressureDofValuesBD.valAsVec()); // t_{comp,dof) = - n_comp * p_dof
2672
2673
2674 // compute \vect t = \sigma \vect n (only vel part) (depending on the flag about symmetric stres?)
2675 ✗ UBlasVector tmp(dimension()),tmp2(m_feVel->numDof());
2676
2677 //\nu \nabla u^T \cdot n
2678 ✗ for (int hbdDof(0); hbdDof<m_curvFeVel->numDof(); hbdDof++) {
2679 ✗ tmp2 = prod(trans(velocityDofValuesVol.val),m_curvFePre->normal[0]);
2680 ✗ tmp = prod(m_feVel->dPhiOnDofs[bd2vol.at(hbdDof)], tmp2);
2681 ✗ for (int jcomp(0); jcomp<dimension(); jcomp++) {
2682 ✗ tractionVel.val(jcomp,hbdDof) = r_instance.viscosity*tmp(jcomp);
2683 }
2684 }
2685 //\nu \nabla u \cdot n
2686 ✗ for (int hbdDof(0); hbdDof<m_curvFeVel->numDof(); hbdDof++) {
2687 ✗ tmp2 = prod(trans(m_feVel->dPhiOnDofs[bd2vol.at(hbdDof)]),m_curvFePre->normal[0]);
2688 ✗ tmp = prod(velocityDofValuesVol.val,tmp2);
2689 ✗ for (int jcomp(0); jcomp<dimension(); jcomp++) {
2690 ✗ tractionVel.val(jcomp,hbdDof) += r_instance.viscosity*tmp(jcomp);
2691 }
2692 }
2693
2694 ✗ m_elementVectorBD[0]->source( /* coef*/-1.0, *m_curvFeVel, tractionVel, /*iBlock*/0, dimension());
2695 ✗ m_elementVectorBD[0]->source( /* coef*/-1.0, *m_curvFePre, tractionPre, /*iBlock*/0, dimension());
2696
2697 // Place m_elementVectorBD[0] in the correct PetscVector
2698 ✗ setValueCustomVectorBD( /* bd element id */ielSupportDof, ADD_VALUES, vec );
2699 ✗ numElement[eltType]++;
2700 }
2701 }
2702 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_rhs);
2703 }
2704 ✗ vec.assembly();
2705 }
2706
2707 ✗ void LinearProblemNS::cvgraphNaturalBC(felInt iel)
2708 {
2709 // Retrieve instance
2710 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
2711
2712 BoundaryCondition* BC;
2713 // Loop over the cvgraph connections
2714 ✗ for ( std::size_t iConn(0); iConn<slave()->numConnections(); ++iConn ) {
2715 // Extract the labels of the current connection
2716 ✗ std::vector<int> labels = slave()->interfaceLabels(iConn);
2717 // Verify if the current label is part of the interface of this connection
2718 ✗ if ( std::find( labels.begin(), labels.end(), m_currentLabel ) != labels.end() ) {
2719 ✗ switch ( slave()->numVarToRead(iConn) ) {
2720
2721 ✗ case 1: // only one var to read: neumann interface BC
2722 // Verify that we have to apply the stress in its strong form.
2723 ✗ if ( slave()->readVariable(iConn,/*cVar*/0) == "STRESS" && slave()->residualInterpolationType(iConn) == 0 ) {
2724 // Look through the different neumann conditions to find the correct one.
2725 ✗ for (std::size_t iNeumann=0; iNeumann<m_boundaryConditionList.numNeumannBoundaryCondition(); iNeumann++ ) {
2726 ✗ BC = m_boundaryConditionList.Neumann(iNeumann);//get the the BC
2727 ✗ if ( std::find( BC->listLabel().begin(), BC->listLabel().end(),m_currentLabel) != BC->listLabel().end() ) {
2728 //compute the force term.
2729 ✗ m_elemFieldNeumann[iNeumann].setValue( m_seqVecs.Get("cvgraphSTRESS"), *m_curvFeVel , iel, m_iVelocity, m_ao, dof());
2730 }
2731 }
2732 }
2733 ✗ break;
2734 ✗ case 2: //two variables to read: robin interface BC
2735 // Look through the different robin conditions to find the correct one.
2736 ✗ for (std::size_t iRobin=0; iRobin<m_boundaryConditionList.numRobinBoundaryCondition(); iRobin++ ) {
2737 ✗ BC = m_boundaryConditionList.Robin(iRobin);
2738 ✗ if ( std::find( BC->listLabel().begin(), BC->listLabel().end(),m_currentLabel ) != BC->listLabel().end() ) {
2739 ✗ if ( slave()->residualInterpolationType(iConn) == 0 ) {
2740 ✗ m_elemFieldRobin[iRobin].setValue( m_seqVecs.Get("cvgraphSTRESS"), *m_curvFeVel , iel, m_iVelocity, m_ao, dof());
2741 } else {
2742 ✗ m_elemFieldRobin[iRobin].val *=0 ;
2743 }
2744 ✗ m_robinAux.setValue(m_seqVecs.Get("cvgraphVelocity"),*m_curvFeVel , iel, m_iVelocity, m_ao, dof());
2745 ✗ m_elemFieldRobin[iRobin].val += r_instance.alphaRobin[iRobin]*m_robinAux.val;
2746 }
2747 }
2748 ✗ break;
2749 ✗ default:
2750 ✗ FEL_ERROR("Three variables to read not yet implemented for this linear problem.")
2751 ✗ break;
2752 }
2753 }
2754 }
2755 }
2756 #endif
2757
2758 /***********************************************************************************/
2759 /***********************************************************************************/
2760
2761 ✗ void LinearProblemNS::localizeItfQpInBackMesh(std::size_t itfMsh, std::size_t /* bckMsh*/, felInt itfVar, std::vector<std::vector<felInt>>& mapItfQpInBackMesh)
2762 {
2763 ✗ const ElementType eltType = m_mesh[itfMsh]->bagElementTypeDomain()[0];
2764 ✗ const int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
2765 ✗ std::vector<Point*> elemPoint(numPointsPerElt, nullptr);
2766 ✗ std::vector<felInt> elemIdPoint(numPointsPerElt, 0);
2767
2768 ✗ const GeoElement *geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
2769 ✗ int typeOfFiniteElement = m_listVariable[itfVar].finiteElementType();
2770 ✗ const RefElement *refEle = geoEle->defineFiniteEle(eltType, typeOfFiniteElement, *m_mesh[itfMsh]);
2771 ✗ CurvilinearFiniteElement feStruc(*refEle, *geoEle, m_listVariable[itfVar].degreeOfExactness());
2772
2773 ✗ const felInt nbrQuaPnt = feStruc.numQuadraturePoint();
2774 ✗ const felInt numEltPerType = m_mesh[itfMsh]->numElements(eltType);
2775
2776 ✗ if ( mapItfQpInBackMesh.size() != static_cast<std::size_t>(numEltPerType) )
2777 ✗ mapItfQpInBackMesh.resize(numEltPerType, std::vector<felInt>(nbrQuaPnt,-1));
2778
2779
2780 // Parallel localization ( i prefer this way since the local mesh is usually owned by only one proc)
2781 ✗ const felInt numEltByProc = numEltPerType/MpiInfo::numProc();
2782 ✗ const felInt beginIdEltLocal = MpiInfo::rankProc()*numEltByProc;
2783 ✗ const felInt numEltProc = ( MpiInfo::rankProc() == MpiInfo::numProc()-1 ) ? numEltPerType-beginIdEltLocal : numEltByProc;
2784
2785 ✗ std::vector<felInt> recvcn(MpiInfo::numProc(), nbrQuaPnt*numEltByProc);
2786 ✗ recvcn.back() = nbrQuaPnt * ( numEltPerType - numEltByProc * ( MpiInfo::numProc() - 1 ) );
2787
2788 ✗ std::vector<felInt> displs(MpiInfo::numProc(), 0);
2789 ✗ std::transform( recvcn.begin(), recvcn.end()-1, displs.begin(), displs.begin()+1, [](auto i, auto j ){ return i + j; } );
2790
2791 ✗ std::vector<Point> listQuadPoint(nbrQuaPnt*numEltProc);
2792
2793 // Loop on element in the region with the type: eltType
2794 ✗ for (felInt iel=0; iel<numEltProc; ++iel) {
2795
2796 // Get elemPoint and elemIdPoint.
2797 ✗ m_mesh[itfMsh]->getOneElement(eltType, beginIdEltLocal+iel, elemIdPoint, elemPoint, 0);
2798
2799 // Update structure finite element
2800 ✗ feStruc.updatePoint(elemPoint);
2801 ✗ feStruc.computeCurrentQuadraturePoint();
2802
2803 // Add quadrature point to list
2804 ✗ for (int ig=0; ig<nbrQuaPnt; ig++)
2805 ✗ listQuadPoint[iel*nbrQuaPnt+ig] = feStruc.currentQuadPoint[ig];
2806 }
2807
2808 ✗ std::vector<felInt> indElt(listQuadPoint.size());
2809 ✗ std::vector<ElementType> typElt(listQuadPoint.size());
2810
2811 // Retrieve instance
2812 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
2813
2814 // Get the timer
2815 ✗ auto& r_timer = r_instance.timer;
2816
2817 ✗ r_timer.Start("Locator" + std::to_string(identifierProblem()) + "::findListPointsInMesh", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
2818 ✗ m_pLocator->findListPointsInMesh(indElt, typElt, listQuadPoint);
2819 ✗ r_timer.Stop("Locator" + std::to_string(identifierProblem()) + "::findListPointsInMesh", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
2820
2821 // Now i need to build the global indElt vector
2822 ✗ std::vector<felInt> indEltGlb(nbrQuaPnt*numEltPerType);
2823 ✗ MPI_Allgatherv(indElt.data(), recvcn[MpiInfo::rankProc()], MPI_INT, indEltGlb.data(), recvcn.data(), displs.data(), MPI_INT, MpiInfo::petscComm() );
2824
2825 ✗ for (felInt iel=0; iel<numEltPerType; iel++){
2826 ✗ for (felInt ig=0; ig<nbrQuaPnt; ig++)
2827 ✗ mapItfQpInBackMesh[iel][ig] = indEltGlb[iel*nbrQuaPnt+ig];
2828 }
2829 }
2830
2831 ✗ void LinearProblemNS::computeSolidFluidMaps()
2832 {
2833 // Retrieve instance
2834 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
2835
2836 ✗ const std::size_t velMesh = m_listVariable[m_iVelocity].idMesh();
2837 ✗ const std::size_t lagMesh = m_listVariable[m_iLagrange].idMesh();
2838
2839 ✗ const ElementType eltType = m_mesh[lagMesh]->bagElementTypeDomain()[0];
2840
2841 // Save previous localization map
2842 ✗ std::vector<std::vector<felInt> > mapLagQpInVelMeshOld(m_mapLagQpInVelMesh.begin(), m_mapLagQpInVelMesh.end());
2843
2844 // Localize structure mesh in fluid mesh
2845 ✗ localizeItfQpInBackMesh(lagMesh, velMesh, m_iLagrange, m_mapLagQpInVelMeshSeq);
2846
2847 // Resize m_mapLagQpInVelMesh if needed
2848 ✗ m_mapLagQpInVelMesh.resize(m_meshLocal[lagMesh]->numElements(eltType));
2849
2850 // Initialize parallel map
2851 ✗ const auto rank = MpiInfo::rankProc();
2852 ✗ auto it = m_eltPart[lagMesh].begin();
2853 ✗ std::copy_if(m_mapLagQpInVelMeshSeq.begin(), m_mapLagQpInVelMeshSeq.end(), m_mapLagQpInVelMesh.begin(), [&](auto&) { return rank == *it++; } );
2854
2855 // Check if localization has changed
2856 ✗ m_hasLocalizationChanged = false;
2857 ✗ if ( m_mapLagQpInVelMesh != mapLagQpInVelMeshOld )
2858 ✗ m_hasLocalizationChanged = true;
2859
2860 ✗ MPI_Allreduce(MPI_IN_PLACE, &m_hasLocalizationChanged, 1, MPI_CXX_BOOL, MPI_LOR, MpiInfo::petscComm());
2861
2862 ✗ if ( r_instance.massConstraint && r_instance.useFicItf ) {
2863
2864 ✗ const std::size_t japMesh = 2; //m_listVariable[m_iLagrange].idMesh(); TODO D.C.
2865
2866 ✗ const ElementType eltType = m_mesh[japMesh]->bagElementTypeDomain()[0];
2867
2868 // Localize structure mesh in fluid mesh
2869 ✗ localizeItfQpInBackMesh(japMesh, velMesh, m_iLagrange, m_mapJapQpInVelMeshSeq);
2870
2871 // Resize m_mapJapQpInVelMesh if needed
2872 ✗ m_mapJapQpInVelMesh.resize(m_meshLocal[japMesh]->numElements(eltType));
2873
2874 // Initialize parallel map
2875 ✗ auto it = m_eltPart[japMesh].begin();
2876 ✗ std::copy_if(m_mapJapQpInVelMeshSeq.begin(), m_mapJapQpInVelMeshSeq.end(), m_mapJapQpInVelMesh.begin(), [&](auto&) { return rank == *it++; } );
2877 }
2878 }
2879
2880 ✗ void LinearProblemNS::evaluateFluidStructurePattern()
2881 {
2882 // Retrieve instance
2883 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
2884
2885 // Get the timer
2886 ✗ auto& r_timer = r_instance.timer;
2887 ✗ r_timer.Start("LinearProblem" + std::to_string(identifierProblem()) + "::evaluateFluidStructurePattern", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
2888
2889 ✗ if ( !m_isCsrNoInterfaceStored ){
2890 ✗ m_csrNoInterface.set(dof().pattern().rowPointer(),dof().pattern().columnIndices());
2891 ✗ m_isCsrNoInterfaceStored = true;
2892 }
2893
2894 ✗ if ( m_hasLocalizationChanged ) {
2895
2896 // Define variable
2897 ✗ std::vector< std::set<felInt> > nodesNeighborhood(m_numDof);
2898 felInt ielSupportDofVel, ielSupportDofLag;
2899 felInt node1, node2, idElt;
2900
2901 ✗ const std::size_t velMesh = m_listVariable[m_iVelocity].idMesh();
2902 ✗ const std::size_t lagMesh = m_listVariable[m_iLagrange].idMesh();
2903
2904 // Get element type
2905 ✗ const ElementType eltType = m_mesh[lagMesh]->bagElementTypeDomain()[0];
2906 ✗ const ElementType eltTypeFl = m_mesh[velMesh]->bagElementTypeDomain()[0];
2907
2908 // Set static pattern
2909 ✗ dof().pattern().set(m_csrNoInterface.rowPointer(),m_csrNoInterface.columnIndices());
2910
2911 // Loop on element in the region with the type: eltType
2912 ✗ const felInt numEltPerType = m_mesh[lagMesh]->numElements(eltType);
2913 ✗ for (felInt iel = 0; iel < numEltPerType; ++iel) {
2914 ✗ m_supportDofUnknown[m_iUnknownLag].getIdElementSupport(eltType, iel, ielSupportDofLag);
2915
2916 // For every LM component
2917 ✗ for (std::size_t iComp = 0; iComp < m_listVariable[m_iLagrange].numComponent(); ++iComp) {
2918 ✗ const int iConnect = dof().getNumGlobComp(m_iUnknownLag, iComp);
2919
2920 // For every dof in support element
2921 ✗ for (int iSup = 0; iSup < m_supportDofUnknown[m_iUnknownLag].getNumSupportDof(ielSupportDofLag); ++iSup) {
2922
2923 // get the global id of the support dof
2924 ✗ dof().loc2glob(ielSupportDofLag, iSup, m_iLagrange, iComp, node1);
2925
2926 // Loop over qpoint
2927 ✗ for (std::size_t iq = 0; iq < m_mapLagQpInVelMeshSeq[iel].size(); ++iq){
2928
2929 ✗ for (std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); ++iUnknown) {
2930 ✗ const int iVariable = m_listUnknown.idVariable(iUnknown);
2931 ✗ if ( m_listVariable[iVariable].idMesh() != velMesh )
2932 ✗ continue;
2933
2934 // get fluid element id
2935 ✗ idElt = m_mapLagQpInVelMeshSeq[iel][iq];
2936
2937 // get id support dof
2938 ✗ m_supportDofUnknown[iUnknown].getIdElementSupport(eltTypeFl, idElt, ielSupportDofVel);
2939
2940 // For every velocity component
2941 ✗ for (std::size_t jComp = 0; jComp < m_listVariable[iVariable].numComponent(); ++jComp) {
2942 ✗ const int jConnect = dof().getNumGlobComp(iUnknown, jComp);
2943
2944 // For every dof in support element
2945 ✗ for (int jSup = 0; jSup < m_supportDofUnknown[iUnknown].getNumSupportDof(ielSupportDofVel); ++jSup) {
2946
2947 // get the global id of the support dof
2948 ✗ dof().loc2glob(ielSupportDofVel, jSup, iVariable, jComp, node2);
2949
2950 // if this node is on this process
2951 ✗ if ( m_listUnknown.mask()(iConnect, jConnect) ) // TODO D.C. not smart put this check here
2952 ✗ nodesNeighborhood[node1].insert(node2);
2953
2954 ✗ if ( m_listUnknown.mask()(jConnect, iConnect) ) // TODO D.C. not smart put this check here
2955 ✗ nodesNeighborhood[node2].insert(node1);
2956 }
2957 }
2958 }
2959 }
2960 }
2961 }
2962 }
2963
2964 // Store the pattern in CSR style
2965 ✗ felInt dofSize = 0;
2966 ✗ felInt cptDof = 0;
2967 felInt pos;
2968
2969 ✗ for(std::size_t iNode=0; iNode<nodesNeighborhood.size(); ++iNode)
2970 ✗ dofSize += nodesNeighborhood[iNode].size();
2971
2972 ✗ CSRMatrixPattern csr(dof().pattern().numRows(), dofSize);
2973 ✗ for(std::size_t iNode = 0; iNode < nodesNeighborhood.size(); ++iNode) {
2974 ✗ if(m_dofPart[iNode] == MpiInfo::rankProc()) {
2975 ✗ csr.rowPointer(cptDof + 1) = csr.rowPointer(cptDof) + nodesNeighborhood[iNode].size();
2976 ✗ pos = 0;
2977 ✗ for (auto it = nodesNeighborhood[iNode].begin(); it != nodesNeighborhood[iNode].end(); ++it) {
2978 ✗ csr.columnIndex(csr.rowPointer(cptDof) + pos) = *it;
2979 ✗ ++pos;
2980 }
2981 ✗ ++cptDof;
2982 }
2983 }
2984
2985 ✗ dof().mergeGlobalPattern(csr.rowPointer(),csr.columnIndices());
2986
2987 ///////////////////////////////////////////////////////////////
2988 // Matrix 0 reconstruction with new pattern
2989 ✗ const std::size_t numRows = dof().pattern().numRows();
2990 ✗ std::vector<felInt> nnz(numRows);
2991 ✗ for (size_t idof = 0; idof < numRows; idof++)
2992 ✗ nnz[idof] = dof().pattern().numNonzerosInRow(idof);
2993
2994 ✗ if ( MpiInfo::numProc() > 1 ) {
2995 ✗ std::vector<felInt> jCSR = dof().pattern().columnIndices();
2996 ✗ AOApplicationToPetsc(ao(), dof().pattern().numNonzeros(), jCSR.data());
2997
2998 ✗ m_matrices[0].destroy();
2999 ✗ m_matrices[0].createAIJ(MpiInfo::petscComm(), m_numDofLocal, m_numDofLocal, dof().numDof(), dof().numDof(), 0, FELISCE_PETSC_NULLPTR, 0, FELISCE_PETSC_NULLPTR);
3000 ✗ m_matrices[0].mpiAIJSetPreallocationCSR(dof().pattern().rowPointer().data(), jCSR.data(), NULL);
3001 ✗ m_matrices[0].setOption(MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
3002 ✗ m_matrices[0].setFromOptions();
3003 ✗ } else {
3004 ✗ m_matrices[0].destroy();
3005 ✗ m_matrices[0].createSeqAIJ(MpiInfo::petscComm(), numRows, numRows, 0, nnz.data());
3006 ✗ m_matrices[0].setOption(MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
3007 ✗ m_matrices[0].setFromOptions();
3008 }
3009 }
3010
3011 ✗ r_timer.Stop("LinearProblem" + std::to_string(identifierProblem()) + "::evaluateFluidStructurePattern", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3012 }
3013
3014 ✗ void LinearProblemNS::computeLagrangeMultiplierOnInterface(FlagMatrixRHS flagMatrixRHS)
3015 {
3016 // Retrieve instance
3017 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
3018
3019 // Get the timer
3020 ✗ auto& r_timer = r_instance.timer;
3021 ✗ r_timer.Start("LinearProblem" + std::to_string(identifierProblem()) + "::computeLagrangeMultiplierOnInterface", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3022
3023 // use to get element number in support dof mesh ordering.
3024 felInt ielSupportDofVel, ielSupportDofLag;
3025
3026 double val;
3027
3028 ✗ const felInt velMesh = m_listVariable[m_iVelocity].idMesh();
3029 ✗ const felInt lagMesh = m_listVariable[m_iLagrange].idMesh();
3030
3031 // Get element type
3032 ✗ std::vector<ElementType> eltTypeVec(m_listVariable.size());
3033 ✗ for (std::size_t iVar = 0; iVar < m_listVariable.size(); ++iVar)
3034 ✗ eltTypeVec[iVar] = m_mesh[m_listVariable[iVar].idMesh()]->bagElementTypeDomain()[0];
3035
3036 ✗ const ElementType eltType = m_mesh[lagMesh]->bagElementTypeDomain()[0];
3037 ✗ const ElementType eltTypeFl = m_mesh[velMesh]->bagElementTypeDomain()[0];
3038
3039 // Resize array.
3040 ✗ const int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
3041 ✗ std::vector<Point*> elemPoint(numPointsPerElt, nullptr);
3042 ✗ std::vector<felInt> elemIdPoint(numPointsPerElt, 0);
3043
3044 ✗ const int numPointsPerEltFl = GeometricMeshRegion::m_numPointsPerElt[eltTypeFl];
3045 ✗ std::vector<Point*> elemPointFl(numPointsPerEltFl, nullptr);
3046 ✗ std::vector<felInt> elemIdPointFl(numPointsPerEltFl, 0);
3047
3048 // Define all current finite element use in the problem from data
3049 ✗ defineCurvilinearFiniteElement(eltTypeVec);
3050 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
3051 ✗ m_curvFeLag = m_listCurvilinearFiniteElement[m_iLagrange];
3052
3053 ///////////////////////////////////////////////////////////////
3054 // Normal constraint
3055 ✗ m_listUnknown.setUnknownsRows({m_iUnknownLag});
3056 ✗ m_listUnknown.setUnknownsCols({m_iUnknownVel});
3057
3058 ✗ m_velBlock = 0;
3059 ✗ m_lagBlock = 0;
3060
3061 // Element matrix and vector initialisation
3062 ✗ initElementArrayBD(true);
3063
3064 // Allocate array use for assemble with petsc.
3065 ✗ allocateArrayForAssembleMatrixRHSBD(flagMatrixRHS);
3066
3067 // Initialize element field
3068 ✗ m_structVelQuad.initialize(QUAD_POINT_FIELD, *m_curvFeLag, m_curvFeLag->numCoor());
3069 ✗ m_structNormalQuad.initialize(QUAD_POINT_FIELD, *m_curvFeLag, m_curvFeLag->numCoor());
3070
3071 // Loop on element in the region with the type: eltType
3072 ✗ const felInt numEltPerType = m_meshLocal[lagMesh]->numElements(eltType);
3073 ✗ for (felInt iel = 0; iel < numEltPerType; ++iel) {
3074
3075 ✗ m_supportDofUnknownLocal[m_iUnknownLag].getIdElementSupport(eltType, iel, ielSupportDofLag);
3076 ✗ ISLocalToGlobalMappingApply(*m_mappingElemSupportPerUnknown[m_iUnknownLag], 1, &ielSupportDofLag, &ielSupportDofLag);
3077
3078 // Get structure element points
3079 ✗ m_meshLocal[lagMesh]->getOneElement(eltType, iel, elemIdPoint, elemPoint, 0);
3080
3081 // Update structure finite element
3082 ✗ m_curvFeLag->updateMeasQuadPt(0, elemPoint);
3083
3084 // **** Matrix ****************
3085 ✗ if ( flagMatrixRHS != FlagMatrixRHS::only_rhs ) {
3086 // Sum over quadrature points
3087 ✗ for (int ig = 0; ig < m_curvFeLag->numQuadraturePoint(); ++ig) {
3088
3089 // Quadrature point localization in fluid mesh
3090 ✗ felInt idElt = m_mapLagQpInVelMesh[iel][ig];
3091
3092 // Get id support dof
3093 ✗ m_supportDofUnknown[m_iUnknownVel].getIdElementSupport(eltTypeFl, idElt, ielSupportDofVel);
3094
3095 // Get fluid element points
3096 ✗ m_mesh[velMesh]->getOneElement(eltTypeFl, idElt, elemIdPointFl, elemPointFl, 0);
3097
3098 // Map structure quadrature point in fluid reference element
3099 Point qRefPoint;
3100 ✗ m_curvFeVel->mapPointInReferenceElement(0, elemPointFl, m_curvFeLag->currentQuadPoint[ig], qRefPoint);
3101
3102 // Clean matrices
3103 ✗ m_elementMatBD[0]->zero();
3104 ✗ m_elementMatTBD[0]->zero();
3105
3106 // Sum over components
3107 ✗ for (std::size_t ic = 0; ic <m_listVariable[m_iLagrange].numComponent(); ++ic){
3108 ✗ UBlasMatrixRange matrix = m_elementMatBD[0]->matBlock(m_lagBlock+ic,m_velBlock+ic);
3109 ✗ UBlasMatrixRange matrixT = m_elementMatTBD[0]->matBlock(m_velBlock+ic,m_lagBlock+ic);
3110
3111 // Sum over structure element nodes
3112 ✗ for(int idof = 0; idof < m_curvFeLag->numDof(); ++idof) {
3113
3114 // Sum over fluid element nodes
3115 ✗ for(int jdof = 0; jdof < m_curvFeVel->numDof(); ++jdof) {
3116
3117 // Compute integral
3118 ✗ val = m_curvFeLag->weightMeas(ig) * m_curvFeLag->phi[ig](idof) * m_curvFeVel->refEle().basisFunction().phi(jdof,qRefPoint);
3119
3120 ✗ matrix(idof,jdof) -= val;
3121 ✗ matrixT(jdof,idof) -= val;
3122 }
3123 }
3124 }
3125
3126 // Add values of elemMat in the global matrix: m_matrices[0].
3127 ✗ setValueCustomMatrixBD(ielSupportDofLag, ielSupportDofVel, m_matrices[0], true);
3128 }
3129 }
3130 // ****************************
3131
3132 // **** Vector ****************
3133 ✗ if ( flagMatrixRHS != FlagMatrixRHS::only_matrix ) {
3134 // Update structure velocity
3135 ✗ updateStructureVel(*m_curvFeLag, elemIdPoint);
3136
3137 // Clean vector
3138 ✗ m_elementVectorBD[0]->zero();
3139
3140 // Compute
3141 ✗ m_elementVectorBD[0]->source(-1., *m_curvFeLag, m_structVelQuad, m_lagBlock, m_listVariable[m_iLagrange].numComponent());
3142
3143 // Add values of elemVec in the global vector: m_vectors[0].
3144 ✗ setValueCustomVectorBD(ielSupportDofLag, ADD_VALUES, m_vectors[0]);
3145 }
3146 // ****************************
3147 }
3148
3149 // Deallocate array use for assemble with petsc.
3150 ✗ desallocateArrayForAssembleMatrixRHS(flagMatrixRHS);
3151
3152 ✗ m_listUnknown.setUnknownsRows({m_iUnknownVel,m_iUnknownPre});
3153 ✗ m_listUnknown.setUnknownsCols({m_iUnknownVel,m_iUnknownPre});
3154
3155 ✗ m_velBlock = m_listUnknown.getBlockPosition(m_iUnknownVel,0);
3156 ✗ m_lagBlock = m_listUnknown.getBlockPosition(m_iUnknownLag,0);
3157
3158 ✗ if( flagMatrixRHS == FlagMatrixRHS::only_rhs ) {
3159 // Real assembling of the right hand side (RHS).
3160 // for (std::size_t i = 0; i < m_vectors.size(); i++)
3161 ✗ m_vectors[0].assembly();
3162 }
3163
3164 ✗ r_timer.Stop("LinearProblem" + std::to_string(identifierProblem()) + "::computeLagrangeMultiplierOnInterface", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3165 }
3166
3167 ✗ void LinearProblemNS::computeBHStabLagLag(const double coeff, const int iq)
3168 {
3169 // \int_{\Sigma} \lambda \cdot \cxi
3170 ✗ UBlasMatrix tmpMat = coeff * m_curvFeLag->weightMeas(iq) * outer_prod(m_curvFeLag->phi[iq], m_curvFeLag->phi[iq]);
3171 ✗ for (std::size_t ic = 0; ic < m_listVariable[m_iLagrange].numComponent(); ++ic){
3172 ✗ UBlasMatrixRange matrix = m_elementMatBD[0]->matBlock(m_lagBlock+ic,m_lagBlock+ic);
3173 ✗ matrix -= tmpMat;
3174 }
3175 }
3176
3177 ✗ void LinearProblemNS::computeBHStabLagJap(const double coeff, const int iq)
3178 {
3179 double tmp;
3180
3181 // \int_{\Sigma} \lambda \cdot \cxi_j n and \int_{\Sigma} \lambda_j n \cdot \cxi
3182 ✗ for (std::size_t ic = 0; ic < m_listVariable[m_iLagrange].numComponent(); ++ic){
3183 ✗ UBlasMatrixRange matrix = m_elementMatBD[0]->matBlock(m_lagBlock+ic,m_japBlock);
3184 ✗ UBlasMatrixRange matrixT = m_elementMatBD[0]->matBlock(m_japBlock,m_lagBlock+ic);
3185
3186 ✗ for (int idof = 0; idof < m_curvFeLag->numDof(); ++idof){
3187 ✗ tmp = coeff * m_curvFeLag->weightMeas(iq) * m_curvFeLag->phi[iq](idof) * m_structNormalQuad.val(ic,iq);
3188 ✗ matrix(idof,0) -= tmp;
3189 ✗ matrixT(0,idof) -= tmp;
3190 }
3191 }
3192 }
3193
3194 ✗ void LinearProblemNS::computeBHStabJapJap(const double coeff, const int iq)
3195 {
3196 // \int_{\Sigma} \lambda_j n \cdot \cxi_j n
3197 ✗ UBlasMatrixRange matrix = m_elementMatBD[0]->matBlock(m_japBlock,m_japBlock);
3198 ✗ matrix(0,0) -= coeff * m_curvFeLag->weightMeas(iq);
3199 }
3200
3201 ✗ void LinearProblemNS::assembleBarbosaHughesStab()
3202 {
3203 // Retrieve instance
3204 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
3205
3206 // Get the timer
3207 ✗ auto& r_timer = r_instance.timer;
3208 ✗ r_timer.Start("LinearProblem" + std::to_string(identifierProblem()) + "::assembleBarbosaHughesStab", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3209
3210 ✗ std::vector<felInt> ielSupportDof;
3211 double coeff;
3212 ✗ const felInt lagMesh = m_listVariable[m_iLagrange].idMesh();
3213
3214 // Get element type
3215 ✗ std::vector<ElementType> eltTypeVec(m_listVariable.size());
3216 ✗ for (std::size_t iVar = 0; iVar < m_listVariable.size(); ++iVar)
3217 ✗ eltTypeVec[iVar] = m_mesh[m_listVariable[iVar].idMesh()]->bagElementTypeDomain()[0];
3218
3219 ✗ const ElementType eltType = m_mesh[lagMesh]->bagElementTypeDomain()[0];
3220
3221 // Resize array.
3222 ✗ const int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
3223 ✗ std::vector<Point*> elemPoint(numPointsPerElt, nullptr);
3224 ✗ std::vector<felInt> elemIdPoint(numPointsPerElt, 0);
3225
3226 // Define all current finite element use in the problem from data
3227 ✗ defineCurvilinearFiniteElement(eltTypeVec);
3228 ✗ m_curvFeLag = m_listCurvilinearFiniteElement[m_iLagrange];
3229
3230 ✗ const double inv_coeff = 1. / ( r_instance.stabCoeffLag * r_instance.viscosity );
3231
3232 felInt japMesh;
3233 ElementType eltTypeFl;
3234 ✗ if ( r_instance.massConstraint ){
3235 ✗ ielSupportDof.resize(2);
3236
3237 ✗ japMesh = m_listVariable[m_iJapanCon].idMesh();
3238 ✗ eltTypeFl = m_mesh[japMesh]->bagElementTypeDomain()[0];
3239
3240 ✗ m_listUnknown.setUnknownsRows({m_iUnknownLag, m_iUnknownJap});
3241 ✗ m_listUnknown.setUnknownsCols({m_iUnknownLag, m_iUnknownJap});
3242 ✗ m_lagBlock = 0;
3243 ✗ m_japBlock = m_listVariable[m_iLagrange].numComponent();
3244 }
3245 else {
3246 ✗ ielSupportDof.resize(1);
3247
3248 ✗ m_listUnknown.setUnknownsRows({m_iUnknownLag});
3249 ✗ m_listUnknown.setUnknownsCols({m_iUnknownLag});
3250 ✗ m_lagBlock = 0;
3251 }
3252
3253 // Initialize element field
3254 ✗ m_structNormalQuad.initialize(QUAD_POINT_FIELD, *m_curvFeLag, m_curvFeLag->numCoor());
3255
3256 // Element matrix and vector initialisation
3257 ✗ initElementArrayBD();
3258
3259 // Allocate array use for assemble with petsc.
3260 ✗ allocateArrayForAssembleMatrixRHSBD(FlagMatrixRHS::only_matrix);
3261
3262 // Loop on element in the region with the type: eltType
3263 ✗ const felInt numEltPerType = m_meshLocal[lagMesh]->numElements(eltType);
3264 ✗ for (felInt iel=0; iel<numEltPerType; iel++) {
3265
3266 ✗ m_supportDofUnknownLocal[m_iUnknownLag].getIdElementSupport(eltType, iel, ielSupportDof[0]);
3267 ✗ ISLocalToGlobalMappingApply(*m_mappingElemSupportPerUnknown[m_iUnknownLag], 1, &ielSupportDof[0], &ielSupportDof[0]);
3268
3269 // Get structure element points
3270 ✗ m_meshLocal[lagMesh]->getOneElement(eltType, iel, elemIdPoint, elemPoint, 0);
3271
3272 // Update structure finite element
3273 ✗ m_curvFeLag->updateMeas(0, elemPoint);
3274
3275 // Compute coefficiente
3276 ✗ coeff = m_curvFeLag->diameter() * inv_coeff;
3277
3278 // Update normals
3279 ✗ updateNormal(*m_curvFeLag, lagMesh, elemIdPoint);
3280
3281 // Sum over quadrature points
3282 ✗ for (int ig=0; ig<m_curvFeLag->numQuadraturePoint(); ig++) {
3283
3284 // Clear elementary matrix.
3285 ✗ m_elementMatBD[0]->zero();
3286
3287 // Assemble local matrix
3288 ✗ computeBHStabLagLag(coeff, ig);
3289
3290 ✗ if ( r_instance.massConstraint ){
3291 // Quadrature point localization in fluid mesh
3292 ✗ felInt idElt = m_mapLagQpInVelMesh[iel][ig];
3293
3294 // Get id support dof
3295 ✗ m_supportDofUnknown[m_iUnknownJap].getIdElementSupport(eltTypeFl, idElt, ielSupportDof[1]);
3296
3297 // Assemble local matrix
3298 ✗ computeBHStabLagJap(coeff, ig);
3299 ✗ computeBHStabJapJap(coeff, ig);
3300 }
3301
3302 // Add values of elemMat in the global matrix: m_matrices[0].
3303 ✗ setValueCustomMatrixBD(ielSupportDof, ielSupportDof, m_matrices[0]);
3304 }
3305 }
3306
3307 // Desallocate array use for assemble with petsc.
3308 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
3309
3310 // Reset
3311 ✗ m_listUnknown.setUnknownsRows({m_iUnknownVel,m_iUnknownPre}); // TODO write a method to reset to default
3312 ✗ m_listUnknown.setUnknownsCols({m_iUnknownVel,m_iUnknownPre});
3313
3314 ✗ m_lagBlock = m_listUnknown.getBlockPosition(m_iUnknownLag,0);
3315 ✗ if ( r_instance.massConstraint )
3316 ✗ m_japBlock = m_listUnknown.getBlockPosition(m_iUnknownJap,0);
3317
3318 ✗ r_timer.Stop("LinearProblem" + std::to_string(identifierProblem()) + "::assembleBarbosaHughesStab", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3319 }
3320
3321 ✗ void LinearProblemNS::assembleBrezziPitkarantaStab()
3322 {
3323 // Retrieve instance
3324 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
3325
3326 ElementType eltType;
3327 int numPointPerElt;
3328 felInt ielSupportDofLag;
3329
3330 ✗ std::vector<Point*> elemPoint;
3331 ✗ std::vector<felInt> elemIdPoint;
3332 ✗ std::vector<felInt> vectorIdSupport;
3333
3334 // Lagrange multiplier mesh
3335 ✗ const felInt lagMesh = m_listVariable[m_iLagrange].idMesh();
3336
3337 ✗ const double inv_coeff = 1. / ( r_instance.stabCoeffLag * r_instance.viscosity );
3338
3339 ✗ m_listUnknown.setUnknownsRows({m_iUnknownLag});
3340 ✗ m_listUnknown.setUnknownsCols({m_iUnknownLag});
3341
3342 ✗ m_lagBlock = 0;
3343
3344 // Assembly loop.
3345 // First loop on geometric type.
3346 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[lagMesh]->bagElementTypeDomain();
3347 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
3348 ✗ eltType = bagElementTypeDomain[i];
3349
3350 // Resize array.
3351 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
3352 ✗ elemPoint.resize(numPointPerElt, nullptr);
3353 ✗ elemIdPoint.resize(numPointPerElt, 0);
3354
3355 // define all current finite element use in the problem from data
3356 ✗ defineCurvilinearFiniteElement(eltType);
3357 ✗ m_curvFeLag = m_listCurvilinearFiniteElement[m_iLagrange];
3358
3359 // Element matrix and std::vector initialisation
3360 ✗ initElementArrayBD();
3361
3362 // Allocate array use for assemble with petsc.
3363 ✗ allocateArrayForAssembleMatrixRHSBD(FlagMatrixRHS::only_matrix);
3364
3365 // Loop on element in the region with the type: eltType
3366 ✗ const felInt numEltPerType = m_meshLocal[lagMesh]->numElements(eltType);
3367 ✗ for (felInt iel = 0; iel < numEltPerType; ++iel) {
3368
3369 // Get element info
3370 ✗ setElemPoint(eltType, iel, elemPoint, elemIdPoint, vectorIdSupport);
3371
3372 // Update structure finite element
3373 ✗ m_curvFeLag->updateMeasNormalContra(0, elemPoint);
3374
3375 // Set coefficient
3376 ✗ const double coeff = std::pow(m_curvFeLag->diameter(),3)*inv_coeff;
3377
3378 // Loop over all the support elements
3379 ✗ for (std::size_t it = 0; it < vectorIdSupport.size(); it++) {
3380
3381 // get the id of the support
3382 ✗ ielSupportDofLag = vectorIdSupport[it];
3383
3384 // clear elementary matrix.
3385 ✗ m_elementMatBD[0]->zero();
3386
3387 // Compute Brezzi-Pitkaranta stabilization
3388 ✗ m_elementMatBD[0]->grad_phi_i_grad_phi_j(-coeff, *m_curvFeLag, m_lagBlock, m_lagBlock, m_listVariable[m_iLagrange].numComponent());
3389
3390 // Add values of elemMat in the global matrix: m_matrices[0].
3391 ✗ setValueCustomMatrixBD(ielSupportDofLag, ielSupportDofLag, m_matrices[0]);
3392 }
3393 }
3394 // Deallocate array use for assemble with petsc.
3395 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
3396 }
3397
3398 // Reset
3399 ✗ m_listUnknown.setUnknownsRows({m_iUnknownVel,m_iUnknownPre}); // TODO write a method to reset to default
3400 ✗ m_listUnknown.setUnknownsCols({m_iUnknownVel,m_iUnknownPre});
3401
3402 ✗ m_lagBlock = m_listUnknown.getBlockPosition(m_iUnknownLag,0);
3403
3404 ✗ if ( r_instance.massConstraint ) {
3405 felInt numGlbDof;
3406 ✗ const double value = -1.e-7;
3407 ✗ m_dof.loc2glob(0, 0, m_iJapanCon, 0, numGlbDof);
3408 ✗ AOApplicationToPetsc(m_ao, 1, &numGlbDof);
3409 ✗ m_matrices[0].setValues(1, &numGlbDof, 1, &numGlbDof, &value, ADD_VALUES);
3410
3411 // m_elementMat[0]->phi_i_phi_j(-1.e-7, *m_feVel, m_lagBlock, m_lagBlock, m_listVariable[m_iLagrange].numComponent());
3412 }
3413 }
3414
3415 ✗ void LinearProblemNS::assembleJapanConstraintBoundary()
3416 {
3417 // Retrieve instance
3418 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
3419
3420 ✗ if ( r_instance.japanLabels.size() == 0 )
3421 ✗ return;
3422
3423 // Get the timer
3424 ✗ auto& r_timer = r_instance.timer;
3425 ✗ r_timer.Start("LinearProblem" + std::to_string(identifierProblem()) + "::assembleJapanConstraintBoundary", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3426
3427 ✗ m_listUnknown.setUnknownsRows({m_iUnknownVel});
3428 ✗ m_listUnknown.setUnknownsCols({m_iUnknownJap});
3429
3430 ✗ m_velBlock = 0;
3431 ✗ m_japBlock = 0;
3432
3433 ElementType eltType; //geometric element type in the mesh.
3434 ✗ int numPointPerElt = 0; //number of points per geometric element.
3435 ✗ int currentLabel = 0; //number of label domain.
3436 ✗ felInt numEltPerLabel = 0; //number of element for one label and one eltType.
3437
3438 // use to define a "global" numbering of element in the mesh.
3439 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
3440 ✗ for (int ityp=0; ityp< GeometricMeshRegion::m_numTypesOfElement; ityp++) {
3441 ✗ eltType = (ElementType)ityp;
3442 ✗ numElement[eltType] = 0;
3443 }
3444
3445 // contains points of the current element.
3446 ✗ std::vector<Point*> elemPoint;
3447
3448 // contains ids of point of the current element.
3449 ✗ std::vector<felInt> elemIdPoint;
3450
3451 // contains ids of support elements of a mesh element
3452 ✗ std::vector<felInt> vectorIdSupport;
3453
3454 // use to get element number in support dof mesh ordering.
3455 felInt ielSupportDofCurv;
3456
3457 //Assembly loop curvilinear.
3458 ✗ const std::vector<ElementType>& bagElementTypeDomainBoundary = m_meshLocal[m_currentMesh]->bagElementTypeDomainBoundary();
3459 ✗ for (std::size_t i = 0; i < bagElementTypeDomainBoundary.size(); ++i) {
3460 ✗ eltType = bagElementTypeDomainBoundary[i];
3461
3462 //resize array.
3463 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
3464 ✗ elemPoint.resize(numPointPerElt, nullptr);
3465 ✗ elemIdPoint.resize(numPointPerElt, 0);
3466
3467 //define all current finite element use in the problem from data
3468 //file configuration and allocate elemMat and elemVec (question: virtual ?).
3469 ✗ defineCurvilinearFiniteElement(eltType);
3470
3471 // Element matrix and vector initialisation
3472 ✗ initElementArrayBD(true);
3473
3474 //allocate array use for assemble with petsc.
3475 ✗ allocateArrayForAssembleMatrixRHSBD(FlagMatrixRHS::only_matrix);
3476
3477 // get velocity curvilinear fe
3478 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
3479
3480 // Second loop on region of the mesh.
3481 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
3482 ✗ currentLabel = itRef->first;
3483 ✗ numEltPerLabel = itRef->second.second;
3484
3485 ✗ auto bndRefIt = find(r_instance.japanLabels.begin(), r_instance.japanLabels.end(), currentLabel);
3486 ✗ if(bndRefIt == r_instance.japanLabels.end()){
3487 ✗ numElement[eltType] += numEltPerLabel;
3488 ✗ continue;
3489 }
3490
3491 // Third loop on element in the region with the type: eltType. ("real" loop on elements).
3492 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
3493 // return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
3494 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, vectorIdSupport);
3495
3496 // update velcity fe
3497 ✗ m_curvFeVel->updateMeasNormal(0, elemPoint);
3498
3499 // loop over all the support elements
3500 ✗ for (std::size_t it = 0; it < vectorIdSupport.size(); it++) {
3501 // get the id of the support
3502 ✗ ielSupportDofCurv = vectorIdSupport[it];
3503
3504 // clear elementary matrix.
3505 ✗ m_elementMatBD[0]->zero();
3506 ✗ m_elementMatTBD[0]->zero();
3507
3508 ✗ m_elementMatBD[0]->const_phi_i_dot_n(-1., *m_curvFeVel, m_velBlock, m_japBlock);
3509 ✗ m_elementMatTBD[0]->const_phi_j_dot_n(-1., *m_curvFeVel, m_japBlock, m_velBlock);
3510
3511 ✗ setValueCustomMatrixBD(ielSupportDofCurv, ielSupportDofCurv, m_matrices[0], true);
3512 }
3513
3514 ✗ numElement[eltType]++;
3515 }
3516 }
3517
3518 //allocate array use for assemble with petsc.
3519 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
3520 }
3521
3522 ✗ m_listUnknown.setUnknownsRows({m_iUnknownVel,m_iUnknownPre});
3523 ✗ m_listUnknown.setUnknownsCols({m_iUnknownVel,m_iUnknownPre});
3524
3525 ✗ m_velBlock = m_listUnknown.getBlockPosition(m_iUnknownVel,0);
3526 ✗ m_japBlock = m_listUnknown.getBlockPosition(m_iUnknownJap,0);
3527
3528 ✗ r_timer.Stop("LinearProblem" + std::to_string(identifierProblem()) + "::assembleJapanConstraintBoundary", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3529 }
3530
3531 ✗ void LinearProblemNS::assembleJapanConstraintInterface()
3532 {
3533 // Retrieve instance
3534 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
3535
3536 // Get the timer
3537 ✗ auto& r_timer = r_instance.timer;
3538 ✗ r_timer.Start("LinearProblem" + std::to_string(identifierProblem()) + "::assembleJapanConstraintInterface", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3539
3540 felInt ielSupportDofVel;
3541
3542 double val;
3543
3544 ✗ const felInt lagMesh = m_listVariable[m_iLagrange].idMesh();
3545 ✗ const felInt velMesh = m_listVariable[m_iVelocity].idMesh();
3546
3547 ✗ std::vector<felInt> listMeshes;
3548 ✗ std::vector< std::vector<std::vector<felInt> >* > listMapQuadPt;
3549 ✗ if ( r_instance.useFicItf ){
3550
3551 ✗ if ( r_instance.M1G.enable ){
3552 ✗ listMeshes = {lagMesh, 2}; // TODO D.C.
3553 ✗ listMapQuadPt = {&m_mapLagQpInVelMesh, &m_mapJapQpInVelMesh};
3554 }
3555 else {
3556 ✗ listMeshes = {2}; // TODO D.C.
3557 ✗ listMapQuadPt = {&m_mapJapQpInVelMesh};
3558 }
3559 } else {
3560 ✗ listMeshes = {lagMesh};
3561 ✗ listMapQuadPt = {&m_mapLagQpInVelMesh};
3562 }
3563
3564 ✗ m_listUnknown.setUnknownsRows({m_iUnknownVel});
3565 ✗ m_listUnknown.setUnknownsCols({m_iUnknownJap});
3566
3567 ✗ m_velBlock = 0;
3568 ✗ m_japBlock = 0;
3569
3570 ✗ for (auto iMsh = 0u; iMsh < listMeshes.size(); ++iMsh) {
3571 ✗ const felInt id_mesh = listMeshes[iMsh];
3572 ✗ const std::vector<std::vector<felInt> >& mapQuadPt = *listMapQuadPt[iMsh];
3573
3574 // Get element type
3575 ✗ std::vector<ElementType> eltTypeVec(m_listVariable.size());
3576 ✗ for (std::size_t iVar = 0; iVar < m_listVariable.size(); ++iVar)
3577 ✗ eltTypeVec[iVar] = m_mesh[m_listVariable[iVar].idMesh()]->bagElementTypeDomain()[0];
3578
3579 // Get element type
3580 ✗ const ElementType eltType = m_mesh[id_mesh]->bagElementTypeDomain()[0];
3581 ✗ const ElementType eltTypeFl = m_mesh[velMesh]->bagElementTypeDomain()[0];
3582
3583 // Resize array.
3584 ✗ const int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
3585 ✗ std::vector<Point*> elemPoint(numPointsPerElt, nullptr);
3586 ✗ std::vector<felInt> elemIdPoint(numPointsPerElt, 0);
3587
3588 ✗ const int numPointsPerEltFl = GeometricMeshRegion::m_numPointsPerElt[eltTypeFl];
3589 ✗ std::vector<Point*> elemPointFl(numPointsPerEltFl, nullptr);
3590 ✗ std::vector<felInt> elemIdPointFl(numPointsPerEltFl, 0);
3591
3592 // Define all current finite element use in the problem from data
3593 ✗ defineCurvilinearFiniteElement(eltTypeVec);
3594 ✗ m_curvFeLag = m_listCurvilinearFiniteElement[m_iLagrange];
3595 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
3596
3597 // Element matrix and vector initialisation
3598 ✗ initElementArrayBD(true);
3599
3600 // Allocate array use for assemble with petsc.
3601 ✗ allocateArrayForAssembleMatrixRHSBD(FlagMatrixRHS::only_matrix);
3602
3603 // Initialize element field
3604 ✗ m_structNormalQuad.initialize(QUAD_POINT_FIELD, *m_curvFeLag, m_curvFeLag->numCoor());
3605
3606 // Loop on element in the region with the type: eltType
3607 ✗ const felInt numEltPerType = m_meshLocal[id_mesh]->numElements(eltType);
3608 ✗ for (felInt iel = 0; iel < numEltPerType; ++iel) {
3609
3610 // return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
3611 ✗ m_meshLocal[id_mesh]->getOneElement(eltType, iel, elemIdPoint, elemPoint, 0);
3612
3613 // Update structure finite element
3614 ✗ m_curvFeLag->updateMeasQuadPt(0, elemPoint);
3615
3616 // Update normals
3617 ✗ updateNormal(*m_curvFeLag, id_mesh, elemIdPoint);
3618
3619 // Sum over quadrature points
3620 ✗ for (int ig = 0; ig < m_curvFeLag->numQuadraturePoint(); ++ig) {
3621
3622 // Quadrature point localization in fluid mesh
3623 ✗ felInt idElt = mapQuadPt[iel][ig];
3624
3625 // Get id support dof
3626 ✗ m_supportDofUnknown[m_iUnknownVel].getIdElementSupport(eltTypeFl, idElt, ielSupportDofVel);
3627
3628 // Get fluid element points
3629 ✗ m_mesh[velMesh]->getOneElement(eltTypeFl, idElt, elemIdPointFl, elemPointFl, 0);
3630
3631 // Map structure quadrature point in fluid reference element
3632 Point qRefPoint;
3633 ✗ m_curvFeVel->mapPointInReferenceElement(0, elemPointFl, m_curvFeLag->currentQuadPoint[ig], qRefPoint);
3634
3635 // Clear elementary matrices
3636 ✗ m_elementMatBD[0]->zero();
3637 ✗ m_elementMatTBD[0]->zero();
3638
3639 // Sum over components
3640 ✗ for (std::size_t ic = 0; ic < m_listVariable[m_iVelocity].numComponent(); ++ic){
3641
3642 // Get matrix block
3643 ✗ UBlasMatrixRange matrix = m_elementMatBD[0]->matBlock(m_velBlock+ic,m_japBlock);
3644 ✗ UBlasMatrixRange matrixT = m_elementMatTBD[0]->matBlock(m_japBlock,m_velBlock+ic);
3645
3646 // Sum over element nodes
3647 ✗ for(int idof = 0; idof < m_curvFeVel->numDof(); ++idof) {
3648
3649 // Compute integral
3650 ✗ val = m_curvFeLag->weightMeas(ig) * m_structNormalQuad.val(ic, ig) * m_curvFeVel->refEle().basisFunction().phi(idof,qRefPoint);
3651
3652 // Fill matrix block
3653 ✗ matrix(idof,0) -= val;
3654 ✗ matrixT(0,idof) -= val;
3655 }
3656 }
3657
3658 // add values of elemMat in the global matrix: m_matrices[0].
3659 ✗ setValueCustomMatrixBD(ielSupportDofVel, ielSupportDofVel, m_matrices[0], true);
3660 }
3661 }
3662
3663 // Deallocate array use for assemble with petsc.
3664 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
3665 }
3666
3667 ✗ m_listUnknown.setUnknownsRows({m_iUnknownVel,m_iUnknownPre});
3668 ✗ m_listUnknown.setUnknownsCols({m_iUnknownVel,m_iUnknownPre});
3669
3670 ✗ m_velBlock = m_listUnknown.getBlockPosition(m_iUnknownVel,0);
3671 ✗ m_japBlock = m_listUnknown.getBlockPosition(m_iUnknownJap,0);
3672
3673 ✗ r_timer.Stop("LinearProblem" + std::to_string(identifierProblem()) + "::assembleJapanConstraintInterface", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3674 }
3675
3676 ✗ void LinearProblemNS::updateStructureVel(CurBaseFiniteElement& fe, std::vector<felInt>& elemIdPoint)
3677 {
3678 // update on the structure dof the local structure velocity from the master at the struct element strucId
3679 // FEL_CHECK(m_structVel != 0, "Error: the pointer of the interface velocity has not been initialized yet.");
3680 // velocity comming from the FSI master ( The order is not the same! )
3681 // the local ordering is done by component while the global is by point
3682 ✗ for(int iq = 0; iq < fe.numQuadraturePoint(); ++iq) {
3683 ✗ for(int ic = 0; ic < fe.numCoor(); ++ic) {
3684 ✗ m_structVelQuad.val(ic, iq) = 0;
3685 ✗ for(int il = 0; il < fe.numDof(); ++il)
3686 ✗ m_structVelQuad.val(ic, iq) += (*m_structVel)[fe.numCoor()*elemIdPoint[il]+ic] * fe.phi[iq](il);
3687 }
3688 }
3689 }
3690
3691 ✗ void LinearProblemNS::computeItfDisp(std::vector<double>& itfDisp)
3692 {
3693 // Retrieve instance
3694 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
3695
3696 // Get number coor
3697 ✗ const int nbrCrd = dimension();
3698
3699 // Get id mesh
3700 ✗ const int idMshLag = m_listVariable[m_iLagrange].idMesh();
3701
3702 // Get number of real interface points
3703 ✗ const int nbrPntReal = m_mesh[idMshLag]->numPoints();
3704
3705 // Compute displacement real structure
3706 ✗ m_realDisp.assign( itfDisp.begin(), itfDisp.begin()+nbrCrd*nbrPntReal );
3707
3708 ✗ if ( r_instance.useFicItf ) {
3709
3710 // Get number of fictitious interface points
3711 ✗ const felInt nbrPntFict = m_mesh[2]->listPoints().size(); // TODO D.C.
3712
3713 // Resize fictitious displacement vector
3714 ✗ m_fictDisp.assign(nbrPntFict*nbrCrd,0.);
3715
3716 // User defined displacement for fictitious interface
3717 ✗ userComputeItfDisp();
3718 }
3719 }
3720
3721 ✗ void LinearProblemNS::moveItfMeshes()
3722 {
3723 // Retrieve instance
3724 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
3725
3726 ✗ const int idMshLag = listVariable()[m_iLagrange].idMesh();
3727
3728 ✗ mesh(idMshLag)->moveMesh(m_realDisp,1.0); // TODO D.C. mesh or meshLocal???
3729 ✗ meshLocal(idMshLag)->moveMesh(m_realDisp,1.0); // TODO D.C. mesh or meshLocal???
3730
3731 // Compute mesh normals
3732 ✗ mesh(idMshLag)->computeNormalTangent();
3733 ✗ auto& listNormals = mesh(idMshLag)->listNormals();
3734 ✗ if ( r_instance.flipItfNormal )
3735 ✗ std::for_each(listNormals.begin(), listNormals.end(), [](auto& p) { p *= -1.; } );
3736
3737 ✗ if ( r_instance.useFicItf ) {
3738
3739 // Move the fictitous mesh
3740 ✗ userMoveFictMesh();
3741
3742 // Compute mesh normals
3743 ✗ mesh(2)->computeNormalTangent();
3744 ✗ auto& listNormals = mesh(2)->listNormals();
3745 ✗ if ( r_instance.flipItfNormal )
3746 ✗ std::for_each(listNormals.begin(), listNormals.end(), [](auto& p) { p *= -1.; } );
3747 }
3748 }
3749
3750 ✗ void LinearProblemNS::updateNormal(CurvilinearFiniteElement& fe, int idMsh, std::vector<felInt>& elemIdPoint)
3751 {
3752 // Retrieve instance
3753 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
3754
3755 ✗ if ( r_instance.useAverageNormal ) {
3756 double norm;
3757 ✗ for(int iq=0; iq<fe.numQuadraturePoint(); ++iq) {
3758 ✗ for(int ic=0; ic<fe.numCoor(); ++ic) {
3759 ✗ m_structNormalQuad.val(ic, iq) = 0;
3760 ✗ for(int idof=0; idof<fe.numDof(); ++idof)
3761 ✗ m_structNormalQuad.val(ic, iq) += fe.phi[iq](idof) * m_mesh[idMsh]->listNormal(elemIdPoint[idof])[ic];
3762 }
3763
3764 ✗ norm = 0;
3765 ✗ for(int ic=0; ic<fe.numCoor(); ++ic)
3766 ✗ norm += m_structNormalQuad.val(ic, iq) * m_structNormalQuad.val(ic, iq);
3767
3768 ✗ norm = 1./std::sqrt(norm);
3769 ✗ for(int ic=0; ic<fe.numCoor(); ++ic)
3770 ✗ m_structNormalQuad.val(ic, iq) *= norm;
3771 }
3772 } else {
3773 ✗ fe.computeNormal();
3774
3775 ✗ const double coef = r_instance.flipItfNormal == true ? -1. : 1.;
3776
3777 ✗ for(int ic=0; ic<fe.numCoor(); ++ic)
3778 ✗ for(int iq=0; iq<fe.numQuadraturePoint(); ++iq)
3779 ✗ m_structNormalQuad.val(ic, iq) = coef * fe.normal[iq](ic);
3780 }
3781 }
3782
3783 ✗ void LinearProblemNS::setInterfaceExchangedData(std::vector<double>& itfVel, std::vector<double>& itfForc)
3784 {
3785 ✗ m_structVel = &itfVel;
3786 ✗ m_intfForc = &itfForc;
3787 }
3788
3789 ✗ void LinearProblemNS::computeInterfaceStress()
3790 {
3791 // Retrieve instance
3792 ✗ auto& r_instance = FelisceParam::instance(instanceIndex());
3793
3794 // Get the timer
3795 ✗ auto& r_timer = r_instance.timer;
3796 ✗ r_timer.Start("LinearProblem" + std::to_string(identifierProblem()) + "::computeInterfaceStress", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3797
3798 ✗ FEL_CHECK(m_intfForc != nullptr, "Error: the pointer of the interface force has not been initialized yet.");
3799
3800 ✗ std::fill(m_intfForc->begin(), m_intfForc->end(), 0);
3801
3802 // const felInt velMesh = m_listVariable[m_iVelocity].idMesh();
3803 ✗ const felInt lagMesh = m_listVariable[m_iLagrange].idMesh();
3804
3805 // Get element type
3806 ✗ std::vector<ElementType> eltTypeVec(m_listVariable.size());
3807 ✗ for (std::size_t iVar = 0; iVar < m_listVariable.size(); ++iVar)
3808 ✗ eltTypeVec[iVar] = m_mesh[m_listVariable[iVar].idMesh()]->bagElementTypeDomain()[0];
3809
3810 ✗ const ElementType eltType = m_mesh[lagMesh]->bagElementTypeDomain()[0];
3811 // const ElementType eltTypeFl = m_mesh[velMesh]->bagElementTypeDomain()[0];
3812
3813 // Resize array
3814 ✗ const int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
3815 ✗ std::vector<Point*> elemPoint(numPointsPerElt, nullptr);
3816 ✗ std::vector<felInt> elemIdPoint(numPointsPerElt, 0);
3817
3818 // Define structure finite element
3819 ✗ const GeoElement *geoEleStruc = m_mesh[lagMesh]->eltEnumToFelNameGeoEle[eltType].second;
3820 ✗ const RefElement *refEleStruc = geoEleStruc->defineFiniteEle(eltType, 0, *m_mesh[lagMesh]);
3821 ✗ m_curvFeStruc = new CurvilinearFiniteElement(*refEleStruc, *geoEleStruc, DegreeOfExactness_2); // TODO D.C.: Provide the value with the fluid or by data file
3822
3823 // Define all current finite element use in the problem from data
3824 ✗ defineCurvilinearFiniteElement(eltTypeVec);
3825 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
3826 ✗ m_curvFeLag = m_listCurvilinearFiniteElement[m_iLagrange];
3827
3828 // Initialize element fields
3829 ✗ m_structNormalQuad.initialize( QUAD_POINT_FIELD, *m_curvFeStruc, m_curvFeStruc->numCoor());
3830 ✗ m_structVelQuad.initialize(QUAD_POINT_FIELD, *m_curvFeStruc, m_curvFeStruc->numCoor());
3831 ✗ m_seqVelDofPts.initialize(DOF_FIELD, *m_curvFeVel, m_curvFeVel->numCoor());
3832 ✗ m_seqLagDofPts.initialize(DOF_FIELD, *m_curvFeLag, m_curvFeLag->numCoor());
3833
3834 ✗ ElementVector elemVec({m_curvFeStruc}, {static_cast<std::size_t>(m_curvFeStruc->numCoor())});
3835
3836 // Loop on element in the region with the type: eltType
3837 ✗ const felInt numEltPerType = m_meshLocal[lagMesh]->numElements(eltType);
3838 ✗ for (felInt iel = 0; iel < numEltPerType; ++iel) {
3839
3840 // return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
3841 ✗ m_meshLocal[lagMesh]->getOneElement(eltType, iel, elemIdPoint, elemPoint, 0);
3842
3843 // clear elementary vector.
3844 ✗ elemVec.zero();
3845
3846 // compute stress on interface element
3847 ✗ computeElementaryInterfaceStress(elemPoint, elemIdPoint, iel, elemVec);
3848
3849 // assembly m_intForc
3850 ✗ for (int ic = 0; ic < m_curvFeStruc->numCoor(); ++ic) {
3851 ✗ for (int ip = 0; ip < m_curvFeStruc->numPoint(); ++ip) {
3852 ✗ int gdof = elemIdPoint[ip]; // <- // TODO D.C.: don't like this, works only for P1
3853 ✗ (*m_intfForc)[m_curvFeStruc->numCoor()*gdof+ic] -= elemVec[ip+ic*m_curvFeStruc->numDof()];
3854 }
3855 }
3856 }
3857
3858 ✗ MPI_Allreduce(MPI_IN_PLACE, m_intfForc->data(), m_intfForc->size(), MPI_DOUBLE, MPI_SUM, MpiInfo::petscComm());
3859
3860 // if ( MpiInfo::rankProc() == 0 ){
3861 // std::cout << "Stress: ";
3862 // for (auto it = m_intfForc->begin(); it != m_intfForc->end(); ++it)
3863 // std::cout << *it << " ";
3864 // std::cout << std::endl;
3865 // }
3866
3867 ✗ delete m_curvFeStruc;
3868
3869 ✗ r_timer.Stop("LinearProblem" + std::to_string(identifierProblem()) + "::computeInterfaceStress", m_fstransient ? (m_fstransient->total_number_iterations > 0 ? m_fstransient->total_number_iterations : m_fstransient->iteration ) : 0);
3870 }
3871
3872 ✗ void LinearProblemNS::computeElementaryInterfaceStress(const std::vector<Point*>& elemPoint, std::vector<felInt>& /*elemIdPoint*/, felInt iel, ElementVector& elemVec)
3873 {
3874 double lag_ic;
3875 felInt ielGlb;
3876 ✗ const felInt lagMesh = m_listVariable[m_iLagrange].idMesh();
3877
3878 // update structure finite element
3879 ✗ m_curvFeStruc->updateMeasQuadPt(0, elemPoint);
3880
3881
3882 //--- Constraint on velocity
3883 ✗ ISLocalToGlobalMappingApply(mappingElem(lagMesh), 1, &iel, &ielGlb);
3884 ✗ m_seqLagDofPts.setValue(sequentialSolution(), *m_curvFeLag, ielGlb, m_iLagrange, m_ao, dof());
3885
3886 // for every component of the struct test function
3887 ✗ for (int ic = 0; ic < m_curvFeStruc->numCoor(); ++ic) {
3888 ✗ UBlasVectorRange vec = elemVec.vecBlock(ic);
3889
3890 // for every quadrature point
3891 ✗ for (int ig = 0; ig < m_curvFeStruc->numQuadraturePoint(); ig++) {
3892
3893 ✗ const Point& qRefPoint = m_curvFeStruc->quadratureRule().quadraturePoint(ig);
3894
3895 // evaluate \lambda_{ic}(x_{iq})
3896 ✗ lag_ic = 0;
3897 ✗ for(int idof = 0; idof < m_curvFeLag->numDof(); idof++)
3898 ✗ lag_ic += m_curvFeLag->refEle().basisFunction().phi(idof,qRefPoint) * m_seqLagDofPts.val(ic,idof);
3899
3900 // for every node
3901 ✗ for(int idof = 0; idof < m_curvFeStruc->numDof(); idof++)
3902 ✗ vec(idof) += m_curvFeStruc->weightMeas(ig) * lag_ic * m_curvFeStruc->phi[ig](idof);
3903 }
3904 }
3905 }
3906
3907 }
3908