GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemNitscheXFEM.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 1674 0.0%
Branches: 0 2671 0.0%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: M. Fernandez & F.M. Gerosa
13 //
14
15 /*!
16 \file linearProblemNitscheXFEM.cpp
17 \authors M. Fernandez & F.M. Gerosa
18 \date 04/2018
19 \brief Solver for a fictitious domain and Nitsche-XFEM fluid formulation.
20 */
21
22 // System includes
23 #include <numeric>
24
25 // External includes
26
27 // Project includes
28 #include "Core/felisceTools.hpp"
29 #include "FiniteElement/elementMatrix.hpp"
30 #include "FiniteElement/elementVector.hpp"
31 #include "Geometry/geo_utilities.hpp"
32 #include "Geometry/neighFacesOfEdges.hpp"
33 #include "linearProblemNitscheXFEM.hpp"
34 #include "Tools/math_utilities.hpp"
35
36 namespace felisce
37 {
38 ✗ LinearProblemNitscheXFEM::LinearProblemNitscheXFEM():
39 LinearProblem("NXFEM Fluid Solver",1,2),
40 ✗ m_viscosity(0.),
41 ✗ m_density(0.),
42 ✗ m_bdf(NULL)
43 {
44 ✗ m_isSeqBdfRHSAllocated = false;
45 ✗ m_isSeqVelExtrapolAllocated = false;
46
47 ✗ m_curvFeStruc = nullptr;
48
49 ✗ m_seqBdfRHS.resize(FelisceParam::instance().orderBdfNS);
50 ✗ for(std::size_t i=0; i<m_seqBdfRHS.size(); ++i)
51 ✗ m_seqBdfRHS[i] = PetscVector::null();
52
53 ✗ m_seqVelExtrapol.resize(FelisceParam::instance().orderBdfNS);
54 ✗ for(std::size_t i=0; i<m_seqVelExtrapol.size(); ++i)
55 ✗ m_seqVelExtrapol[i] = PetscVector::null();
56
57 ✗ m_intfVelo = nullptr;
58 ✗ m_intfForc = nullptr;
59 }
60
61 ✗ LinearProblemNitscheXFEM::~LinearProblemNitscheXFEM()
62 {
63
64 ✗ m_bdf = nullptr;
65 ✗ if (m_isSeqBdfRHSAllocated) {
66 ✗ for(std::size_t i=0; i<m_seqBdfRHS.size(); ++i)
67 ✗ m_seqBdfRHS[i].destroy();
68 ✗ m_seqBdfRHS.clear();
69 }
70
71 ✗ if (m_isSeqVelExtrapolAllocated) {
72 ✗ for(std::size_t i=0; i<m_seqVelExtrapol.size(); ++i)
73 ✗ m_seqVelExtrapol[i].destroy();
74 ✗ m_seqVelExtrapol.clear();
75 }
76
77 ✗ for(std::size_t iPt=0; iPt < m_mshPts.size(); ++iPt) {
78 ✗ delete m_mshPts[iPt];
79 }
80
81 ✗ for(std::size_t iPt=0; iPt < m_subItfPts.size(); ++iPt) {
82 ✗ delete m_subItfPts[iPt];
83 }
84
85 ✗ for(std::size_t i=0; i<m_solutionOld.size(); ++i)
86 ✗ m_solutionOld[i].destroy();
87
88 ✗ for(std::size_t i=0; i<m_sequentialSolutionOld.size(); ++i)
89 ✗ m_sequentialSolutionOld[i].destroy();
90
91 ✗ for(std::size_t i=0; i<m_aoOld.size(); ++i)
92 ✗ AODestroy(&m_aoOld[i]);
93 }
94
95 ✗ void LinearProblemNitscheXFEM::initialize(std::vector<GeometricMeshRegion::Pointer>& vec_mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES)
96 {
97 // calling initialize in LinearProblem
98 ✗ LinearProblem::initialize(vec_mesh, comm, doUseSNES);
99
100 // the problem
101 ✗ m_fstransient = fstransient;
102
103 // define variable of the linear system
104 ✗ std::vector<PhysicalVariable> listVariable(2);
105 ✗ std::vector<std::size_t> listNumComp(2);
106 ✗ listVariable[0] = velocity;
107 ✗ listNumComp[0] = dimension();
108 ✗ listVariable[1] = pressure;
109 ✗ listNumComp[1] = 1;
110
111 //define unknown of the linear system.
112 ✗ m_listUnknown.push_back(velocity);
113 ✗ m_listUnknown.push_back(pressure);
114
115 //adding Darcy pressure on a boundary
116 ✗ if ( FelisceParam::instance().contactDarcy > 0 ) {
117 ✗ listVariable.push_back(pressureDarcy); // this represents the constant LM
118 ✗ listNumComp.push_back(1); // to be used with fusionDof and AllToOne = true
119 ✗ m_listUnknown.push_back(pressureDarcy);
120 }
121
122 // link between variable and unknowns
123 ✗ definePhysicalVariable(listVariable, listNumComp);
124
125 // get id of the variables
126 ✗ m_iVelocity = m_listVariable.getVariableIdList(velocity);
127 ✗ m_iPressure = m_listVariable.getVariableIdList(pressure);
128
129 // get the variables
130 ✗ m_velocity = &m_listVariable[m_iVelocity];
131 ✗ m_pressure = &m_listVariable[m_iPressure];
132
133 // get id of the unknowns
134 ✗ m_iUnknownVel = m_listUnknown.getUnknownIdList(velocity);
135 ✗ m_iUnknownPre = m_listUnknown.getUnknownIdList(pressure);
136
137 // index block unknown
138 ✗ m_velBlock = m_listUnknown.getBlockPosition(m_iUnknownVel,0);
139 ✗ m_preBlock = m_listUnknown.getBlockPosition(m_iUnknownPre,0);
140
141 ✗ if ( FelisceParam::instance().contactDarcy > 0 ) {
142 ✗ m_iPreDarcy = m_listVariable.getVariableIdList(pressureDarcy); // get id of the variable
143 ✗ m_presDarcy = &m_listVariable[m_iPreDarcy]; // get the variable
144 ✗ m_iUnknownDar = m_listUnknown.getUnknownIdList(pressureDarcy); // get id of the unknown
145 ✗ m_darBlock = m_listUnknown.getBlockPosition(m_iUnknownDar,0); // index block unknown
146 }
147
148 // data from the data file
149 ✗ m_viscosity = FelisceParam::instance().viscosity;
150 ✗ m_density = FelisceParam::instance().density;
151 ✗ m_nitschePar = FelisceParam::instance().NitschePenaltyParam;
152
153 ✗ m_useSymmetricStress = FelisceParam::instance().useSymmetricStress;
154 ✗ m_useNSEquation = FelisceParam::instance().NSequationFlag;
155 ✗ m_useInterfaceStab = FelisceParam::instance().useInterfaceFlowStab;
156
157 // number of vertes for sub element
158 ✗ m_numVerPerFluidElt = mesh()->domainDim() + 1;
159 ✗ m_numVerPerSolidElt = mesh()->domainDim();
160
161 // create the vector of points for the fluid elements
162 ✗ m_mshPts.resize(m_numVerPerFluidElt);
163 ✗ for(std::size_t iver=0; iver < m_mshPts.size(); ++iver) {
164 ✗ m_mshPts[iver] = new Point();
165 }
166
167 ✗ m_itfPts.resize(m_numVerPerSolidElt);
168 ✗ m_subItfPts.resize(m_numVerPerSolidElt);
169 ✗ for(std::size_t iver=0; iver < m_subItfPts.size(); ++iver) {
170 ✗ m_subItfPts[iver] = new Point();
171 }
172
173 ✗ m_skipSmallVolume = m_skipSmallVolume && ( dimension() == 3 );
174 ✗ m_printSkipVolume = m_printSkipVolume && m_skipSmallVolume;
175 }
176
177 ✗ void LinearProblemNitscheXFEM::initSupportDofDerivedProblem()
178 {
179
180 ✗ m_duplicateSupportElements->duplicateIntersectedSupportElements(m_supportDofUnknown);
181 }
182
183 ✗ void LinearProblemNitscheXFEM::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS)
184 {
185 IGNORE_UNUSED_ELT_TYPE;
186 IGNORE_UNUSED_FLAG_MATRIX_RHS;
187
188 // Get the finite elements
189 ✗ m_feVel = m_listCurrentFiniteElement[m_iVelocity];
190 ✗ m_fePres = m_listCurrentFiniteElement[m_iPressure];
191
192 // Definition of the ElemField
193 ✗ m_elemFieldAdv.initialize( DOF_FIELD, *m_feVel, m_feVel->numCoor());
194 ✗ m_elemFieldAdvDup.initialize(DOF_FIELD, *m_feVel, m_feVel->numCoor());
195 ✗ m_elemFieldAdvTmp.initialize(DOF_FIELD, *m_feVel, m_feVel->numCoor());
196 ✗ m_elemFieldRHSbdf.initialize(DOF_FIELD, *m_feVel, m_velocity->numComponent());
197 ✗ m_elemFieldNormal.initialize(CONSTANT_FIELD, m_feVel->numCoor());
198
199 // solid
200 ✗ m_strucVelQuadPts.initialize(QUAD_POINT_FIELD, *m_feVel, m_velocity->numComponent()); // in computeElementArray
201 ✗ m_strucVelDofPts.initialize(DOF_FIELD, *m_curvFeStruc, m_curvFeStruc->numCoor()); // in computeElementArray
202 ✗ m_ddotComp.initialize(QUAD_POINT_FIELD, *m_fePres, m_pressure->numComponent()); // in computeElementArray
203 }
204
205 ✗ void LinearProblemNitscheXFEM::initPerElementTypeBoundaryCondition(ElementType& eltType, FlagMatrixRHS flagMatrixRHS)
206 {
207 IGNORE_UNUSED_ELT_TYPE;
208 IGNORE_UNUSED_FLAG_MATRIX_RHS;
209
210 // Get the curvilinear finite elements
211 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
212 ✗ m_curvFePress = m_listCurvilinearFiniteElement[m_iPressure];
213 }
214
215 ✗ void LinearProblemNitscheXFEM::userChangePattern(int numProc, int rankProc)
216 {
217 (void) numProc;
218 (void) rankProc;
219 // There may be integrals that involve the dof of the two duplicated support elements.
220 // The integrals for the ghost penalty modify also the pattern
221 // same for the face-oriented stabilization and the DG terms for the tip
222
223 // definition of variables
224 felInt currentIelGeo, oppositeIelGeo;
225 felInt node1, node2;
226 felInt idVar1, idVar2;
227 felInt ielSupport, jelSupport;
228 ✗ std::vector<felInt> vecSupportCurrent, vecSupportOpposite;
229 bool isAnInnerFace;
230
231 GeometricMeshRegion::ElementType eltType, eltTypeOpp;
232 felInt numEltPerLabel;
233 felInt numElement[GeometricMeshRegion::m_numTypesOfElement];
234 felInt numFacesPerElement;
235
236 ✗ const std::vector<felInt>& iSupportDof = m_supportDofUnknown[0].iSupportDof();
237 ✗ const std::vector<felInt>& iEle = m_supportDofUnknown[0].iEle();
238
239 // the initial repartition of dof
240 ✗ const felInt numDof = dof().numDof();
241 ✗ const felInt numDofByProc = numDof/MpiInfo::numProc();
242 ✗ const felInt beginIdDofLocal = MpiInfo::rankProc()*numDofByProc;
243 ✗ const felInt numDofProc = ( MpiInfo::rankProc() == MpiInfo::numProc()-1 ) ? numDof-beginIdDofLocal : numDofByProc;
244 ✗ const felInt endIdDofLocal = beginIdDofLocal + numDofProc;
245
246 ✗ std::vector<std::set<felInt> > nodesNeighborhood(numDofProc);
247
248 // build the edges or faces depending on the dimension
249 // In this function, "faces" refers to either edges or faces.
250 ✗ if( FelisceParam::instance().useGhostPenalty || FelisceParam::instance().NSStabType == 1 || !FelisceParam::instance().confIntersection ) {
251 ✗ if(dimension() == 2)
252 ✗ mesh()->buildEdges();
253 else
254 ✗ mesh()->buildFaces();
255 }
256
257 // ------------------------------------------------------------------------- //
258 // Complementary pattern duplicated elements //
259 // ------------------------------------------------------------------------- //
260 ✗ const auto& idInterElt = m_duplicateSupportElements->getIntersectedMshElt();
261 ✗ for(auto ielIt = idInterElt.begin(); ielIt != idInterElt.end(); ++ielIt) {
262 // get the geometric id of the element and the ids of its support elements
263 ✗ currentIelGeo = ielIt->first;
264
265 // loop over the unknowns of the problem
266 ✗ for (std::size_t iUnknown1 = 0; iUnknown1 < m_listUnknown.size(); ++iUnknown1) {
267 ✗ idVar1 = m_listUnknown.idVariable(iUnknown1);
268 // get the support element of the current element
269 ✗ m_supportDofUnknown[iUnknown1].getIdElementSupport(currentIelGeo, vecSupportCurrent);
270
271 // loop over the support elements
272 ✗ for(std::size_t ielSup = 0; ielSup < vecSupportCurrent.size(); ++ielSup) {
273 ✗ ielSupport = vecSupportCurrent[ielSup];
274 // for all support dof in this support element
275 ✗ for (felInt iSup = 0; iSup < m_supportDofUnknown[iUnknown1].getNumSupportDof(ielSupport); ++iSup) {
276 // for all component of the unknown
277 ✗ for (std::size_t iComp = 0; iComp < m_listVariable[idVar1].numComponent(); ++iComp) {
278 // get the global id of the support dof
279 ✗ dof().loc2glob(ielSupport, iSup, idVar1, iComp, node1);
280
281 // If this support dof is owned by this process
282 ✗ if(node1 >= beginIdDofLocal && node1 < endIdDofLocal) {
283 ✗ const felInt iConnect = dof().getNumGlobComp(iUnknown1, iComp);
284
285 // for all unknown in the linear problem
286 ✗ for (std::size_t iUnknown2 = 0; iUnknown2 < m_listUnknown.size(); ++iUnknown2) {
287 ✗ idVar2 = m_listUnknown.idVariable(iUnknown2);
288 // for all component of the second unknown
289 ✗ for (std::size_t jComp = 0; jComp < m_listVariable[idVar2].numComponent(); ++jComp) {
290 // If the two current components are connected
291 ✗ const felInt jConnect = dof().getNumGlobComp(iUnknown2, jComp);
292 ✗ if (m_listUnknown.mask()(iConnect, jConnect) > 0) {
293 // for all support element different than the ielSup
294 ✗ for(std::size_t jelSup=0; jelSup<vecSupportCurrent.size(); ++jelSup) {
295 ✗ if(jelSup != ielSup) {
296 ✗ jelSupport = vecSupportCurrent[jelSup];
297 // for all support dof in this support element
298 ✗ for (felInt jSup = 0; jSup < m_supportDofUnknown[iUnknown2].getNumSupportDof(jelSupport); ++jSup) {
299 // get the global id of the second support dof
300 ✗ dof().loc2glob(jelSupport, jSup, idVar2, jComp, node2);
301 ✗ if(node1 != node2)
302 ✗ nodesNeighborhood[node1-beginIdDofLocal].insert(node2);
303 }
304 }
305 }
306 }
307 }
308 }
309 }
310 }
311 }
312 }
313 }
314 }
315
316 // ----------------------------------------------------- //
317 // Complementary pattern for face-oriented stabilization //
318 // Complementary pattern for ghost penalty stabilization //
319 // ----------------------------------------------------- //
320 ✗ if( FelisceParam::instance().NSStabType == 1 || FelisceParam::instance().useGhostPenalty ) {
321
322 // first loop on element type
323 ✗ currentIelGeo = 0;
324 ✗ for (std::size_t i=0; i<mesh()->bagElementTypeDomain().size(); ++i) {
325 ✗ eltType = mesh()->bagElementTypeDomain()[i];
326 ✗ const GeoElement* geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
327 ✗ numElement[eltType] = 0;
328 ✗ numFacesPerElement = geoEle->numBdEle();
329 // second loop on region of the mesh.
330 ✗ for(GeometricMeshRegion::IntRefToBegEndIndex_type::const_iterator itRef = mesh()->intRefToBegEndMaps[eltType].begin(); itRef != mesh()->intRefToBegEndMaps[eltType].end(); itRef++) {
331 ✗ numEltPerLabel = itRef->second.second;
332 // third loop on element
333 ✗ for (felInt iElm = 0; iElm < numEltPerLabel; ++iElm) {
334
335 // loop over the boundaries of the element
336 ✗ for(felInt iFace = 0; iFace < numFacesPerElement; ++iFace) {
337 // check if this face is an inner face or a boundary
338 ✗ oppositeIelGeo = currentIelGeo;
339 ✗ eltTypeOpp = eltType;
340 ✗ isAnInnerFace = mesh()->getAdjElement(eltTypeOpp, oppositeIelGeo, iFace);
341 ✗ if(isAnInnerFace) {
342
343 // loop on the unknown of the linear problem
344 ✗ for (std::size_t iUnknown1 = 0; iUnknown1 < m_listUnknown.size(); ++iUnknown1) {
345 ✗ idVar1 = m_listUnknown.idVariable(iUnknown1);
346 // get the support element of the current element
347 ✗ m_supportDofUnknown[iUnknown1].getIdElementSupport(currentIelGeo, vecSupportCurrent);
348
349 // get number of dof in support element
350 ✗ const std::size_t sizeCurrent = m_supportDofUnknown[0].getNumSupportDof(vecSupportCurrent[0]);
351
352 // loop over the support elements
353 ✗ for(std::size_t ielSup = 0; ielSup < vecSupportCurrent.size(); ++ielSup) {
354 ✗ ielSupport = vecSupportCurrent[ielSup];
355
356 // for all support dof in this support element
357 ✗ for (felInt iSup = 0; iSup < m_supportDofUnknown[iUnknown1].getNumSupportDof(ielSupport); ++iSup) {
358 // for all component of the unknown
359 ✗ for (std::size_t iComp = 0; iComp < m_listVariable[idVar1].numComponent(); ++iComp) {
360 // get the global id of the support dof
361 ✗ dof().loc2glob(ielSupport, iSup, idVar1, iComp, node1);
362
363 // if this node is on this process
364 ✗ if(node1 >= beginIdDofLocal && node1 < endIdDofLocal) {
365 ✗ const felInt iConnect = dof().getNumGlobComp(iUnknown1, iComp);
366
367 // for all unknown
368 ✗ for (std::size_t iUnknown2 = 0; iUnknown2 < m_listUnknown.size(); ++iUnknown2) {
369 ✗ idVar2 = m_listUnknown.idVariable(iUnknown2);
370
371 // get the support element of the opposite element
372 ✗ m_supportDofUnknown[iUnknown2].getIdElementSupport(oppositeIelGeo, vecSupportOpposite);
373
374 // get number of dof in support element
375 ✗ const std::size_t sizeOpposite = m_supportDofUnknown[0].getNumSupportDof(vecSupportOpposite[0]);
376
377 ✗ for(std::size_t jelSup=0; jelSup < vecSupportOpposite.size(); ++jelSup) {
378 ✗ jelSupport = vecSupportOpposite[jelSup];
379
380 // compute dof in common
381 ✗ std::size_t countSupportDofInCommon = 0;
382 ✗ for(std::size_t iSup1=0; iSup1<sizeCurrent; ++iSup1) {
383 ✗ for(std::size_t iSup2=0; iSup2<sizeOpposite; ++iSup2)
384 ✗ if( iSupportDof[iEle[jelSupport] + iSup2] == iSupportDof[iEle[ielSupport] + iSup1] )
385 ✗ ++countSupportDofInCommon;
386 }
387
388 // if the two support elements share an edge/face
389 ✗ if( countSupportDofInCommon == sizeCurrent-1 ) {
390 // for all component of the second unknown
391 ✗ for (std::size_t jComp = 0; jComp < m_listVariable[idVar2].numComponent(); jComp++) {
392 // If the two current components are connected
393 ✗ const felInt jConnect = dof().getNumGlobComp(iUnknown2, jComp);
394 ✗ if (m_listUnknown.mask()(iConnect, jConnect) > 0) {
395 // for all support dof in this support element
396 ✗ for (felInt jSup = 0; jSup < m_supportDofUnknown[iUnknown2].getNumSupportDof(jelSupport); ++jSup) {
397 // get the global id of the second support dof
398 ✗ dof().loc2glob(jelSupport, jSup, idVar2, jComp, node2);
399 ✗ if(node1 != node2)
400 ✗ nodesNeighborhood[node1-beginIdDofLocal].insert(node2);
401 }
402 }
403 }
404 }
405 }
406 }
407 }
408 }
409 }
410 }
411 }
412 }
413 }
414 ✗ ++currentIelGeo;
415 ✗ ++numElement[eltType];
416 }
417 }
418 }
419 }
420
421 // ----------------------------------------------------- //
422 // Complementary pattern for tip-DG faces //
423 // ----------------------------------------------------- //
424 ✗ const felInt nbrTipElt = m_duplicateSupportElements->getNumTipElt();
425 ✗ if(nbrTipElt > 0 && !FelisceParam::instance().confIntersection ) {
426
427 // get set of tip elements
428 ✗ const auto& idTipElt = m_duplicateSupportElements->getIntersectedTipElt();
429
430 // copy of idTipElt in temporary structure
431 ✗ std::unordered_set<felInt> idTipAndNeiElt(idTipElt.begin(),idTipElt.end());
432
433 // get elttype and number of faces
434 ✗ eltType = mesh()->bagElementTypeDomain()[0];
435 ✗ numFacesPerElement = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numBdEle();
436
437 // loop on tip elements
438 ✗ for (auto it = idTipElt.begin(); it != idTipElt.end(); ++it) {
439 ✗ currentIelGeo = *it;
440
441 // looping on the faces of the element
442 ✗ for(felInt iFace=0; iFace < numFacesPerElement; ++iFace) {
443 // check if this face is an inner face or a boundary
444 ✗ oppositeIelGeo = currentIelGeo;
445 ✗ eltTypeOpp = eltType;
446 ✗ isAnInnerFace = mesh()->getAdjElement(eltTypeOpp, oppositeIelGeo, iFace);
447
448 ✗ if( isAnInnerFace )
449 ✗ idTipAndNeiElt.insert(oppositeIelGeo);
450 }
451 }
452
453 // loop on tip elements and neighbors
454 ✗ for (auto it = idTipAndNeiElt.begin(); it != idTipAndNeiElt.end(); ++it) {
455 ✗ currentIelGeo = *it;
456
457 ✗ bool isCurrentIelGeoATipElm = false;
458 ✗ if ( idTipElt.find(currentIelGeo) != idTipElt.end() )
459 ✗ isCurrentIelGeoATipElm = true;
460
461 // looping on the faces of the element
462 ✗ for(felInt iFace=0; iFace < numFacesPerElement; ++iFace) {
463 // check if this face is an inner face or a boundary
464 ✗ oppositeIelGeo = currentIelGeo;
465 ✗ eltTypeOpp = eltType;
466 ✗ isAnInnerFace = mesh()->getAdjElement(eltTypeOpp, oppositeIelGeo, iFace);
467
468 // if boundary face go to next one
469 ✗ if ( !isAnInnerFace )
470 ✗ continue;
471
472 // if current element is not a tip elt, check if the opposite is a tip elt
473 ✗ bool isOppositeIelGeoATipElm = false;
474 ✗ if ( !isCurrentIelGeoATipElm ) {
475 ✗ if ( idTipElt.find(oppositeIelGeo) != idTipElt.end() )
476 ✗ isOppositeIelGeoATipElm = true;
477 }
478
479 // if one of the two is a tip elt
480 ✗ if( isCurrentIelGeoATipElm || isOppositeIelGeoATipElm ) {
481 // loop on the unknown of the linear problem
482 ✗ for (std::size_t iUnknown1 = 0; iUnknown1 < m_listUnknown.size(); ++iUnknown1) {
483 ✗ idVar1 = m_listUnknown.idVariable(iUnknown1);
484 // get the support element of the current element
485 ✗ m_supportDofUnknown[iUnknown1].getIdElementSupport(currentIelGeo, vecSupportCurrent);
486
487 // loop over the support elements
488 ✗ for(std::size_t ielSup = 0; ielSup < vecSupportCurrent.size(); ++ielSup) {
489 ✗ ielSupport = vecSupportCurrent[ielSup];
490
491 // for all support dof in this support element
492 ✗ for (felInt iSup = 0; iSup < m_supportDofUnknown[iUnknown1].getNumSupportDof(ielSupport); ++iSup) {
493 // for all component of the unknown
494 ✗ for (std::size_t iComp = 0; iComp < m_listVariable[idVar1].numComponent(); ++iComp) {
495 // get the global id of the support dof
496 ✗ dof().loc2glob(ielSupport, iSup, idVar1, iComp, node1);
497
498 // if this node is on this process
499 ✗ if(node1 >= beginIdDofLocal && node1 < endIdDofLocal) {
500 ✗ const felInt iConnect = dof().getNumGlobComp(iUnknown1, iComp);
501
502 // for all unknown
503 ✗ for (std::size_t iUnknown2 = 0; iUnknown2 < m_listUnknown.size(); ++iUnknown2) {
504 ✗ idVar2 = m_listUnknown.idVariable(iUnknown2);
505 // get the support element of the opposite element
506 ✗ m_supportDofUnknown[iUnknown2].getIdElementSupport(oppositeIelGeo, vecSupportOpposite);
507
508 ✗ for(std::size_t jelSup=0; jelSup < vecSupportOpposite.size(); ++jelSup) {
509 ✗ jelSupport = vecSupportOpposite[jelSup];
510
511 // for all component of the second unknown
512 ✗ for (std::size_t jComp = 0; jComp < m_listVariable[idVar2].numComponent(); jComp++) {
513 // If the two current components are connected
514 ✗ const felInt jConnect = dof().getNumGlobComp(iUnknown2, jComp);
515 ✗ if (m_listUnknown.mask()(iConnect, jConnect) > 0) {
516 // for all support dof in this support element
517 ✗ for (felInt jSup = 0; jSup < m_supportDofUnknown[iUnknown2].getNumSupportDof(jelSupport); ++jSup) {
518 // get the global id of the second support dof
519 ✗ dof().loc2glob(jelSupport, jSup, idVar2, jComp, node2);
520 ✗ if(node1 != node2)
521 ✗ nodesNeighborhood[node1-beginIdDofLocal].insert(node2);
522 }
523 }
524 }
525 }
526 }
527 }
528 }
529 }
530 }
531 }
532 }
533 }
534 }
535 }
536
537 // ----------------------------------------------------------------- //
538 // build the complementary pattern and merge it with the current one //
539 // ----------------------------------------------------------------- //
540 ✗ std::size_t dofSize = 0;
541 ✗ for(std::size_t iNode=0; iNode<nodesNeighborhood.size(); ++iNode)
542 ✗ dofSize += nodesNeighborhood[iNode].size();
543
544 // the complementary pattern
545 ✗ if(dofSize > 0) {
546 ✗ std::vector<felInt> iCSR(numDofProc+1, 0);
547 ✗ std::vector<felInt> jCSR(dofSize , 0);
548 felInt pos;
549
550 ✗ for(felInt iNode = 0; iNode < numDofProc; ++iNode) {
551 ✗ iCSR[iNode + 1] = iCSR[iNode] + nodesNeighborhood[iNode].size();
552 ✗ pos = 0;
553 ✗ for(auto iCon = nodesNeighborhood[iNode].begin(); iCon != nodesNeighborhood[iNode].end(); ++iCon) {
554 ✗ jCSR[iCSR[iNode] + pos] = *iCon;
555 ✗ ++pos;
556 }
557 }
558
559 // Now, call merge pattern
560 ✗ dof().mergeGlobalPattern(iCSR, jCSR);
561 }
562 }
563
564 ✗ void LinearProblemNitscheXFEM::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel1, felInt& iel2, ElementType& eltType, felInt& ielGeo, FlagMatrixRHS flagMatrixRHS)
565 {
566 IGNORE_UNUSED_FLAG_MATRIX_RHS;
567 IGNORE_UNUSED_ELEM_ID_POINT;
568
569 // get the side of the support element looking at iel2
570 // here iel2 is refering to the solution because it will set the column in the global matrix
571 // (see the function setValueMatrixRHS in linearProblem.cpp)
572 felInt ielSupportDof;
573 ✗ m_supportDofUnknown[0].getIdElementSupport(eltType, ielGeo, ielSupportDof); // TODO here ielGeo is ID elt by type in the rest of the function is the global ID
574 ✗ const sideOfInterface sideOfSupportElm = ielSupportDof == iel2 ? LEFT : RIGHT;
575 ✗ const double side = sideOfSupportElm == LEFT ? 1. : -1.;
576
577 ✗ if( m_skipSmallVolume ){
578 ✗ m_totalVol = compAreaVolume(elemPoint);
579 ✗ m_skippedVol = 0;
580 ✗ m_totalArea = 0;
581 ✗ m_skippedArea = 0;
582 ✗ m_tol = m_eps * m_totalVol;
583 }
584
585 ✗ if( iel1 == iel2 ) {
586 //////////////////////////////////
587 // VOLUME ELEMENT - SUB ELEMENT //
588 //////////////////////////////////
589 ✗ felInt numSubElement = m_duplicateSupportElements->getNumSubMshEltPerMshElt(ielGeo);
590
591 // Get the support element of the old solutions (id of duplicated/non-dupl old element in the old support dof)
592 ✗ std::vector<felInt> ielOld(m_solutionOld.size());
593 ✗ std::vector< std::vector<felInt> > oldSupportElement(m_solutionOld.size());
594 ✗ for(std::size_t i = 0; i < oldSupportElement.size(); ++i) {
595 ✗ m_supportDofUnknownOld[i][0].getIdElementSupport(ielGeo, oldSupportElement[i]);
596
597 // set the id for the old solution
598 ✗ if( oldSupportElement[i].size() > 1 ) {
599 // this element is also duplicated for the old solution
600 ✗ ielOld[i] = sideOfSupportElm == LEFT ? oldSupportElement[i][0] : oldSupportElement[i][1];
601 } else {
602 // this element is NOT duplicated for the old solution
603 ✗ ielOld[i] = oldSupportElement[i][0];
604 }
605 }
606
607 // This element has been duplicted in the actual configuration
608 ✗ if( numSubElement > 0 ) {
609
610 // Volume sub elements
611 felInt subElemId;
612 ✗ std::vector<double> tmpPt(3);
613
614 // loop over the sub-elements
615 ✗ for(felInt iSubElt = 0; iSubElt < numSubElement; ++iSubElt) {
616 ✗ subElemId = m_duplicateSupportElements->getSubMshEltIdxMsh(ielGeo, iSubElt);
617
618 // if the side of the sub-element is the same as the element
619 ✗ if( m_duplicateSupportElements->getSubMshEltSide(ielGeo, iSubElt) == sideOfSupportElm ) {
620
621 // fill the vector with the points of the sub-element
622 ✗ for(std::size_t iPt = 0; iPt < m_numVerPerFluidElt; ++iPt) {
623 ✗ m_duplicateSupportElements->getSubMshEltVerCrd(subElemId, iPt, tmpPt);
624 ✗ *m_mshPts[iPt] = Point(tmpPt[0], tmpPt[1], tmpPt[2]);
625 }
626
627 ✗ if( m_skipSmallVolume ){
628 ✗ const double area = compAreaVolume(m_mshPts);
629 ✗ if( area <= m_tol ){
630 ✗ m_skippedVol += std::abs(area);
631 ✗ continue;
632 }
633 }
634
635 // Update finite element
636 ✗ m_feVel->updateSubElementFirstDeriv(0, elemPoint, m_mshPts);
637 ✗ m_fePres->updateSubElementFirstDeriv(0, elemPoint, m_mshPts);
638
639 // build classic NS matrix
640 ✗ computeElementMatrixRHS();
641 ✗ for(felInt i = 0; i < FelisceParam::instance().orderBdfNS; ++i)
642 ✗ computeElementMatrixRHSPartDependingOnOldTimeStep(iel1, ielOld[i], i);
643 }
644 }
645 ✗ } else {
646 // No volume sub elements
647 ✗ for(felInt i = 0; i < FelisceParam::instance().orderBdfNS; ++i) {
648
649 ✗ if( oldSupportElement[i].size() > 1 ) {
650 // old element is duplicated, get the sub elements
651 ✗ const felInt numSubElementOld = m_duplicateSupportElements->getNumSubMshEltMshOld(ielGeo, i);
652
653 felInt subElemIdOld, ielOldTmp;
654 ✗ std::vector<double> tmpPtOld(3);
655
656 ✗ for(felInt iSubElt = 0; iSubElt < numSubElementOld; ++iSubElt) {
657 ✗ subElemIdOld = m_duplicateSupportElements->getSubMshEltIdxMshOld(ielGeo, iSubElt, i);
658
659 // fill the vector with the points of the sub-element
660 ✗ for(std::size_t iPt = 0; iPt < m_numVerPerFluidElt; ++iPt) {
661 ✗ m_duplicateSupportElements->getSubMshEltVerCrdOld(subElemIdOld, iPt, tmpPtOld, i);
662 ✗ *m_mshPts[iPt] = Point(tmpPtOld[0], tmpPtOld[1], tmpPtOld[2]);
663 }
664
665 ✗ if( m_skipSmallVolume ){
666 ✗ const double area = compAreaVolume(m_mshPts);
667 ✗ if(area <= m_tol){
668 ✗ m_skippedVol += std::abs(area);
669 ✗ continue;
670 }
671 }
672
673 // Update finite element
674 ✗ m_feVel->updateSubElementFirstDeriv(0, elemPoint, m_mshPts);
675 ✗ m_fePres->updateSubElementFirstDeriv(0, elemPoint, m_mshPts);
676
677 // Get the id of the support element owning the sub element
678 ✗ ielOldTmp = m_duplicateSupportElements->getSubMshEltSideOld(subElemIdOld, i) == LEFT ? oldSupportElement[i][0] : oldSupportElement[i][1];
679
680 // build classic NS matrix
681 ✗ computeElementMatrixRHSPartDependingOnOldTimeStep(iel1, ielOldTmp, i);
682 }
683 }
684 }
685
686 // Update finite element if needed (needed in case we just modified the fem due to integration over subDomains)
687 ✗ if( !m_feVel->hasOriginalQuadPoint() || !m_fePres->hasOriginalQuadPoint() ) {
688 ✗ m_feVel->updateBasisAndFirstDeriv(0, elemPoint);
689 ✗ m_fePres->updateBasisAndFirstDeriv(0, elemPoint);
690 } else {
691 ✗ m_feVel->updateFirstDeriv(0, elemPoint);
692 ✗ m_fePres->updateFirstDeriv(0, elemPoint);
693 }
694
695 // build classic NS matrix
696 ✗ computeElementMatrixRHS();
697
698 // build matrix and rhs depending on old time step if this element is not intersected by the old solutions
699 ✗ for(felInt i = 0; i < FelisceParam::instance().orderBdfNS; ++i) {
700 ✗ if( oldSupportElement[i].size() == 1 )
701 ✗ computeElementMatrixRHSPartDependingOnOldTimeStep(iel1, ielOld[i], i);
702 }
703 }
704
705
706 ////////////////////////////
707 // SUB STRUCTURE ELEMENTS //
708 ////////////////////////////
709 ✗ numSubElement = m_duplicateSupportElements->getNumSubItfEltPerMshElt(ielGeo);
710 ✗ if( numSubElement > 0 ) {
711 felInt subElemId;
712 felInt elemId;
713 ✗ std::vector<double> tmpPt(3);
714 ✗ std::vector<double> itfNormal(dimension(),0);
715 ✗ std::vector<felInt> strucElemPtsId(m_numVerPerSolidElt);
716 ✗ GeometricMeshRegion::ElementType eltTypeStruc = m_interfMesh->bagElementTypeDomain()[0];
717
718 // new position for penaltyParam
719 ✗ m_feVel->updateMeas(0, elemPoint);
720 ✗ const double penaltyParam = m_viscosity * m_nitschePar / m_feVel->diameter();
721
722 // for every interface sub-element
723 ✗ for(felInt iSubElt = 0; iSubElt < numSubElement; ++iSubElt) {
724 // get the id of the sub element
725 ✗ subElemId = m_duplicateSupportElements->getSubItfEltIdxMsh(ielGeo, iSubElt);
726
727 // get the id of the element owing this sub-element
728 ✗ elemId = m_duplicateSupportElements->getItfEltIdxOfSubItfElt(subElemId);
729
730 // skip this subElement if belongs to the fictitious interface (this is a tip)
731 ✗ if( m_duplicateSupportElements->getItfEltRef(elemId) < 0 )
732 ✗ continue;
733
734 // fill the vector with the points of the sub-element
735 ✗ for(std::size_t iPt = 0; iPt < m_numVerPerSolidElt; ++iPt) {
736 ✗ m_duplicateSupportElements->getSubItfEltVerCrd(subElemId, iPt, tmpPt);
737 ✗ *m_subItfPts[iPt] = Point(tmpPt[0], tmpPt[1], tmpPt[2]);
738 }
739
740 // get the coordinates of the element owning this sub-element
741 ✗ m_interfMesh->getOneElement(eltTypeStruc, elemId, strucElemPtsId, m_itfPts, 0);
742
743 ✗ if( m_skipSmallVolume ){
744 ✗ const double area = compAreaVolume(m_subItfPts);
745 ✗ m_totalArea += area;
746 ✗ if ( area <= m_eps*compAreaVolume(m_itfPts) ) {
747 ✗ m_skippedArea += std::abs(area);
748 ✗ continue;
749 }
750 }
751
752 // update the finite element
753 ✗ m_feVel->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_curvFeStruc);
754 ✗ m_fePres->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_curvFeStruc);
755
756 // Inflow stabilization for NS
757 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab ) {
758
759 // get the normal to this sub element
760 ✗ m_interfMesh->listElementNormal(eltTypeStruc, elemId).getCoor(itfNormal);
761 ✗ m_elemFieldNormal.setValue(itfNormal);
762
763 // compute (\beta \cdot n)_- u_i v_i
764 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(0.5 * side * m_density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
765 }
766
767 // compute the elementary matrix and rhs
768 // update the value of the solid velocity
769 ✗ updateStructureVel(elemId);
770 ✗ updateSubStructureVel(m_itfPts, m_subItfPts);
771
772 // Nitsche's penalty term, \gamma \mu / h u_i . v_i
773 ✗ m_elementMat[0]->phi_i_phi_j(penaltyParam, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
774
775 // Nitsche's consistency terms: - \sigma(u_i, p_i)n_i . v_i
776 ✗ computeNitscheSigma_u_v(side * m_viscosity, -side);
777
778 // Nitsche's consistency terms, - \sigma(v_i, -q_i)n_i . u_i
779 // BE CAREFUL: the equation for q is multiply by minus !!!!!
780 ✗ computeNitscheSigma_v_u(0., -side);
781
782 // Nische's penalty term, \gamma \mu / h \dot{d}^n v_i
783 ✗ m_elementVector[0]->source(penaltyParam, *m_feVel, m_strucVelQuadPts, m_velBlock, m_velocity->numComponent());
784
785 // ddot \sigma(v_i, -q_i) n_i . ddot //
786 ✗ computeNitscheDpoint_Sigma_v(0., -side);
787 }
788 }
789 ✗ } else {
790
791 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab ) {
792
793 // Get the support element of the old solutions
794 ✗ std::vector<felInt> ielOld(m_solutionOld.size());
795 ✗ std::vector<felInt> ielOldDup(m_solutionOld.size());
796 ✗ std::vector< std::vector<felInt> > oldSupportElement(m_solutionOld.size());
797 ✗ for(std::size_t i = 0; i < oldSupportElement.size(); ++i) {
798 ✗ m_supportDofUnknownOld[i][0].getIdElementSupport(ielGeo, oldSupportElement[i]);
799
800 // set the id for the old solution
801 ✗ if( oldSupportElement[i].size() > 1 ) {
802 // this element is also duplicated for the old solution
803 ✗ if( sideOfSupportElm == LEFT ) {
804 ✗ ielOld[i] = oldSupportElement[i][0];
805 ✗ ielOldDup[i] = oldSupportElement[i][1];
806 } else {
807 ✗ ielOld[i] = oldSupportElement[i][1];
808 ✗ ielOldDup[i] = oldSupportElement[i][0];
809 }
810 } else {
811 // this element is NOT duplicated for the old solution
812 ✗ ielOld[i] = oldSupportElement[i][0];
813 ✗ ielOldDup[i] = oldSupportElement[i][0];
814 }
815 }
816
817 ////////////////////////////
818 // SUB STRUCTURE ELEMENTS //
819 ////////////////////////////
820 ✗ const felInt numSubElement = m_duplicateSupportElements->getNumSubItfEltPerMshElt(ielGeo);
821 ✗ if( numSubElement > 0 ) {
822 felInt subElemId, elemId;
823 ✗ std::vector<double> tmpPt(3);
824 ✗ std::vector<double> itfNormal(dimension());
825 ✗ std::vector<felInt> strucElemPtsId(m_numVerPerSolidElt);
826 ✗ GeometricMeshRegion::ElementType eltTypeStruc = m_interfMesh->bagElementTypeDomain()[0];
827
828 // set the extrapolation of velocity
829 ✗ m_elemFieldAdv.val.clear();
830 ✗ m_elemFieldAdvDup.val.clear();
831 ✗ if( FelisceParam::instance().useProjectionForNitscheXFEM ) {
832 ✗ for(felInt i = 0; i < FelisceParam::instance().orderBdfNS; ++i) {
833 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, iel1, m_iVelocity, m_ao, dof());
834 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
835 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, iel2, m_iVelocity, m_ao, dof());
836 ✗ m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
837 }
838 } else {
839 ✗ for(felInt i = 0; i < FelisceParam::instance().orderBdfNS; ++i) {
840 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOld[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
841 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
842 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOldDup[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
843 ✗ m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
844 }
845 }
846
847 // for every sub-element
848 ✗ for(felInt iSubElt = 0; iSubElt < numSubElement; ++iSubElt) {
849 // get the id of the sub element
850 ✗ subElemId = m_duplicateSupportElements->getSubItfEltIdxMsh(ielGeo, iSubElt);
851
852 // get the id of the element owing this sub element
853 ✗ elemId = m_duplicateSupportElements->getItfEltIdxOfSubItfElt(subElemId);
854
855 // skip this subElement if belongs to the fictitious interface (this is a tip)
856 ✗ if( m_duplicateSupportElements->getItfEltRef(elemId) < 0 )
857 ✗ continue;
858
859 // fill the vector with the points of the sub-element
860 ✗ for(std::size_t iPt = 0; iPt < m_numVerPerSolidElt; ++iPt) {
861 ✗ m_duplicateSupportElements->getSubItfEltVerCrd(subElemId, iPt, tmpPt);
862 ✗ *m_subItfPts[iPt] = Point(tmpPt[0], tmpPt[1], tmpPt[2]);
863 }
864
865 // get the coordinates of the element owning this sub-element
866 ✗ m_interfMesh->getOneElement(eltTypeStruc, elemId, strucElemPtsId, m_itfPts, 0);
867
868 ✗ if( m_skipSmallVolume ){
869 ✗ const double area = compAreaVolume(m_subItfPts);
870 ✗ m_totalArea += area;
871 ✗ if ( area <= m_eps*compAreaVolume(m_itfPts) ) {
872 ✗ m_skippedArea += std::abs(area);
873 ✗ continue;
874 }
875 }
876
877 // get the normal to this sub element
878 ✗ m_interfMesh->listElementNormal(eltTypeStruc, elemId).getCoor(itfNormal);
879 ✗ m_elemFieldNormal.setValue(itfNormal);
880
881 // update the finite element
882 ✗ m_feVel->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_curvFeStruc);
883 ✗ m_fePres->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_curvFeStruc);
884
885 // compute Inflow stabilization for NS equation
886 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(0.25 * side * m_density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
887 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(0.25 * side * m_density, *m_feVel, m_elemFieldAdvDup, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
888 }
889 }
890 }
891 }
892
893 ✗ if( m_printSkipVolume ) {
894 ✗ if ( std::abs(m_skippedVol) > m_eps * std::abs(m_totalVol) )
895 ✗ std::cout << "In computeElementArray, iel1 = " << iel1 << " iel2 = " << iel2 << " Skipped volume is : " << std::fixed << std::setprecision(16) << m_skippedVol / m_totalVol * 100. << " %" << std::defaultfloat << " m_skippedVol = " << m_skippedVol << std::endl;
896
897 ✗ if ( std::abs(m_skippedArea) > m_eps * std::abs(m_totalArea) )
898 ✗ std::cout << "In computeElementArray, iel1 = " << iel1 << " iel2 = " << iel2 << " Skipped surface is : " << std::fixed << std::setprecision(16) << m_skippedArea / m_totalArea * 100. << " %" << std::defaultfloat << " m_skippedArea = " << m_skippedArea << std::endl;
899 }
900 }
901
902 namespace
903 {
904 template<int DimT>
905 ✗ void normal(std::vector<Point*>& pnt, std::vector<double>& N) {
906
907 if ( DimT == 2 ){
908 ✗ N[0] = pnt[0]->y() - pnt[1]->y();
909 ✗ N[1] = pnt[1]->x() - pnt[0]->x();
910
911 ✗ const double invNorm = 1./ std::sqrt( N[0]*N[0] + N[1]*N[1]);
912
913 ✗ N[0] *= invNorm;
914 ✗ N[1] *= invNorm;
915 }
916 else if ( DimT == 3 ){
917
918 ✗ MathUtilities::CrossProduct(N, (*pnt[1] - *pnt[0]).getCoor(), (*pnt[2] - *pnt[0]).getCoor());
919
920 ✗ const double invNorm = 1./ std::sqrt( N[0]*N[0] + N[1]*N[1] + N[2]*N[2]);
921
922 ✗ N[0] *= invNorm;
923 ✗ N[1] *= invNorm;
924 ✗ N[2] *= invNorm;
925 }
926 else{
927 std::cout << "ERROR in normal()" << std::endl;
928 }
929 }
930 }
931
932 ✗ void LinearProblemNitscheXFEM::assembleMatrixFrontPoints()
933 {
934 // This is to handle the DG terms for the elements containing the tip of the solid
935 // Assemble the matrix to take the crack point into account
936 // Enforce continuity accros the fictitious interface
937 ✗ const felInt nbrVerFro = m_duplicateSupportElements->getNumFrontPoints();
938
939 ✗ if(nbrVerFro > 0) {
940 ✗ std::vector<felInt> verFro;
941 ✗ std::vector<felInt> verFro_temp; // generic front vertex vector for checking multi fronts within element
942 felInt currentGlobalIelGeo, currentIelGeo;
943 ✗ GeometricMeshRegion::ElementType eltType = mesh()->bagElementTypeDomain()[0];
944
945 ✗ const felInt numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
946 ✗ std::vector<Point*> elemPoint(numPointPerElt);
947 ✗ std::vector<felInt> elemIdPoint(numPointPerElt, 0);
948
949 ✗ std::vector<felInt> vectorIdSupport;
950 ✗ std::vector<double> itfNormal(dimension(), 0.);
951
952 PetscInt naux;
953
954 // Define finite element
955 ✗ defineFiniteElement(eltType);
956
957 // Element matrix and vector initialisation
958 ✗ initElementArray();
959
960 // allocate arrays for assembling the matrix
961 ✗ FlagMatrixRHS flag = FlagMatrixRHS::only_matrix;
962 ✗ allocateArrayForAssembleMatrixRHS(flag);
963
964 // init variables
965 ✗ initPerElementType(eltType, flag);
966
967 ✗ for(felInt iVer=0; iVer<nbrVerFro; ++iVer) {
968 // get the info on the front point
969 ✗ m_duplicateSupportElements->getFrontPointInfo(iVer, verFro);
970
971 // get the fluid element in which the fictitious interface is
972 ✗ currentGlobalIelGeo = verFro[1];
973 ✗ ISGlobalToLocalMappingApply(mappingElem(), IS_GTOLM_MASK, 1, &currentGlobalIelGeo, &naux, &currentIelGeo);
974 ✗ if ( currentIelGeo == -1 )
975 ✗ continue;
976
977 // get the support elements
978 ✗ setElemPoint(eltType, currentIelGeo, elemPoint, elemIdPoint, vectorIdSupport);
979
980 // fill the vector with the points of the fictitious interface
981 ✗ *m_subItfPts[0] = Point(m_interfMesh->listPoint(verFro[0]));
982 ✗ *m_subItfPts[1] = Point(mesh()->listPoint(verFro[2]));
983
984 // now check if two fronts are in same element--then opposite vertex is other front)
985 ✗ for(felInt iVer2=0; iVer2<nbrVerFro; ++iVer2) {
986 ✗ if(iVer2!=iVer){
987 ✗ m_duplicateSupportElements->getFrontPointInfo(iVer2, verFro_temp);
988 ✗ if(verFro[1]==verFro_temp[1]){
989 ✗ *m_subItfPts[1] = Point(m_interfMesh->listPoint(verFro[2]));
990 }
991 }
992 }
993
994 // update the finite element
995 ✗ m_feVel->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_curvFeStruc);
996 ✗ m_fePres->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_curvFeStruc);
997
998 // compute the normal
999 ✗ normal<2>(m_subItfPts, itfNormal);
1000
1001 // setting normal in element field
1002 ✗ m_elemFieldNormal.setValue(itfNormal);
1003
1004 // parameters
1005 ✗ const double penaltyParam = m_viscosity * m_nitschePar / m_feVel->diameter();
1006
1007 // loop over the support elements
1008 ✗ for(std::size_t iSup1=0; iSup1 < vectorIdSupport.size(); ++iSup1) {
1009 ✗ felInt iel1 = vectorIdSupport[iSup1];
1010 ✗ for(std::size_t iSup2=0; iSup2 < vectorIdSupport.size(); ++iSup2) {
1011 ✗ felInt iel2 = vectorIdSupport[iSup2];
1012
1013 // clear elementary matrix
1014 ✗ for (std::size_t j=0; j<numberOfMatrices(); ++j)
1015 ✗ m_elementMat[j]->zero();
1016
1017 // get the side
1018 ✗ const sideOfInterface sideOfSupportElm = iSup2 == 0 ? LEFT : RIGHT;
1019 ✗ const double side = sideOfSupportElm == LEFT ? 1. : -1.;
1020
1021 // set the extrapolation of velocity
1022 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab ) {
1023 // Get the support element of the old solutions
1024 ✗ std::vector<felInt> ielOld(m_solutionOld.size());
1025 ✗ std::vector<felInt> ielOldDup(m_solutionOld.size());
1026 ✗ std::vector< std::vector<felInt> > oldSupportElement(m_solutionOld.size());
1027 ✗ for(std::size_t i=0; i<oldSupportElement.size(); ++i) {
1028 ✗ m_supportDofUnknownOld[i][0].getIdElementSupport(currentGlobalIelGeo, oldSupportElement[i]);
1029
1030 // set the id for the old solution
1031 ✗ if(oldSupportElement[i].size() > 1) {
1032 // this element is also duplicated for the old solution
1033 ✗ if(sideOfSupportElm == LEFT) {
1034 ✗ ielOld[i] = oldSupportElement[i][0];
1035 ✗ ielOldDup[i] = oldSupportElement[i][1];
1036 } else {
1037 ✗ ielOld[i] = oldSupportElement[i][1];
1038 ✗ ielOldDup[i] = oldSupportElement[i][0];
1039 }
1040 } else {
1041 // this element is NOT duplicated for the old solution
1042 ✗ ielOld[i] = oldSupportElement[i][0];
1043 ✗ ielOldDup[i] = oldSupportElement[i][0];
1044 }
1045 }
1046
1047 ✗ m_elemFieldAdv.val.clear();
1048 ✗ m_elemFieldAdvDup.val.clear();
1049 ✗ if(FelisceParam::instance().useProjectionForNitscheXFEM) {
1050 ✗ for(felInt i=0; i<FelisceParam::instance().orderBdfNS; ++i) {
1051 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, iel1, m_iVelocity, m_ao, dof());
1052 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1053 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, iel2, m_iVelocity, m_ao, dof());
1054 ✗ m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
1055 }
1056 } else {
1057 ✗ for(felInt i=0; i<FelisceParam::instance().orderBdfNS; ++i) {
1058 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOld[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1059 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1060 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOldDup[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1061 ✗ m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
1062 }
1063 }
1064 }
1065
1066 // The normal is defined external to Omega_1. (Omega_1: right and Omega_2: left)
1067 // side = 1, on the left and -1 on the right.
1068 // The unknown is telling us in which subdomain we are and therefore which normal we have to consider
1069 ✗ if(iSup1 == iSup2) {
1070
1071 // Nitsche's penalty term, \gamma \nu / h u_i . v_i
1072 ✗ m_elementMat[0]->phi_i_phi_j(penaltyParam, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1073
1074 // Nitsche's consistency terms, -\sigma(u_i, p_i)n_i . v_i - \sigma(v_i, q_i)n_i . u_i
1075 ✗ computeNitscheSigma_u_v(0.5 * side * m_viscosity, -0.5 * side);
1076 ✗ computeNitscheSigma_v_u(0.5 * side * m_viscosity, -0.5 * side);
1077
1078 // stabilize inflow
1079 ✗ if(m_useNSEquation == 1 && m_useInterfaceStab ){
1080 // (0.5 \rho_f \beta_1 \cdot n_2) (u_1 \cdot v_1)-(0.5 \rho_f \beta_2 \cdot n_2) (u_2 \cdot v_2)
1081 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(0.5 * side * m_density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1082 }
1083 } else {
1084 ✗ m_elementMat[0]->phi_i_phi_j(-penaltyParam, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1085
1086 ✗ computeNitscheSigma_u_v(-0.5 * side * m_viscosity, 0.5 * side);
1087 ✗ computeNitscheSigma_v_u( 0.5 * side * m_viscosity, -0.5 * side);
1088
1089 // stabilize inflow
1090 ✗ if(m_useNSEquation == 1 && m_useInterfaceStab ) {
1091 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(0.25 * side * m_density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1092 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(0.25 * side * m_density, *m_feVel, m_elemFieldAdvDup, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1093 }
1094 }
1095
1096 // add values of elemMat in the global matrix: matrix(0).
1097 ✗ setValueMatrixRHS(vectorIdSupport[iSup1], vectorIdSupport[iSup2], flag);
1098 }
1099 }
1100 }
1101
1102 // desallocate array use for assemble with petsc.
1103 ✗ desallocateArrayForAssembleMatrixRHS(flag);
1104 }
1105 }
1106
1107 ✗ void LinearProblemNitscheXFEM::assembleMatrixTip()
1108 {
1109 // This is to handle the DG terms for the elements containing the tip edges/faces of the solid
1110 // Enforce continuity accros the fictitious interface
1111
1112 // get number of fluid tip elements
1113 ✗ const felInt nbrTipElm = m_duplicateSupportElements->getNumTipElt();
1114 ✗ if( nbrTipElm > 0 ) {
1115
1116 // we assume fluid elements to be only triangles or tetrahedra
1117 ✗ GeometricMeshRegion::ElementType eltTypeFluid = meshLocal()->bagElementTypeDomain()[0];
1118 ✗ std::vector<Point*> fluidElmPoint(m_numVerPerFluidElt, nullptr);
1119 ✗ std::vector<felInt> fluidElmIdPoint(m_numVerPerFluidElt);
1120 // we assume structure elements to be only segments or triangles
1121 ✗ GeometricMeshRegion::ElementType eltTypeStruc = m_interfMesh->bagElementTypeDomain()[0];
1122 ✗ std::vector<felInt> strucElmIdPoint(m_numVerPerSolidElt);
1123
1124 felInt currentIelGeo, currentGlobalIelGeo, numSubElement;
1125 felInt elemId, subElemId;
1126
1127 ✗ std::vector<felInt> vectorIdSupport;
1128 ✗ std::vector<double> itfNormal(dimension());
1129 ✗ std::vector<double> tmpPt(3);
1130 PetscInt naux;
1131
1132 // define finite element
1133 ✗ defineFiniteElement(eltTypeFluid);
1134
1135 // Element matrix and vector initialisation
1136 ✗ initElementArray();
1137
1138 // allocate arrays for assembling the matrix
1139 ✗ allocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
1140
1141 // init variables
1142 ✗ initPerElementType(eltTypeFluid, FlagMatrixRHS::only_matrix);
1143
1144 // loop over fluid tip elements
1145 ✗ for(felInt iTipElm = 0; iTipElm < nbrTipElm; ++iTipElm) {
1146 // get fluid tip element global and local id, if element doesn't belong to current proc skip it
1147 ✗ currentGlobalIelGeo = m_duplicateSupportElements->getTipEltIdx(iTipElm);
1148 ✗ ISGlobalToLocalMappingApply(mappingElem(), IS_GTOLM_MASK, 1, &currentGlobalIelGeo, &naux, &currentIelGeo);
1149 ✗ if ( currentIelGeo == -1 )
1150 ✗ continue;
1151
1152 ✗ if( m_skipSmallVolume ){
1153 ✗ m_totalArea = 0;
1154 ✗ m_skippedArea = 0;
1155 }
1156
1157 // get the number of structure sub-element in the current fluid element
1158 ✗ numSubElement = m_duplicateSupportElements->getNumSubItfEltPerMshElt(currentGlobalIelGeo);
1159 ✗ for(felInt iSubElt = 0; iSubElt < numSubElement; ++iSubElt) {
1160 // get the id of the strucure sub-element
1161 ✗ subElemId = m_duplicateSupportElements->getSubItfEltIdxMsh(currentGlobalIelGeo, iSubElt);
1162 // get the id of the structure element owing this sub-element
1163 ✗ elemId = m_duplicateSupportElements->getItfEltIdxOfSubItfElt(subElemId);
1164 // if elemId doesn't belong to fictitious structure skip it
1165 ✗ if( m_duplicateSupportElements->getItfEltRef(elemId) > 0 )
1166 ✗ continue;
1167
1168 // get the normal to this sub element
1169 ✗ m_interfMesh->listElementNormal(eltTypeStruc, elemId).getCoor(itfNormal);
1170
1171 // fill the vector with the points of the sub-element
1172 ✗ for(std::size_t iPt = 0; iPt < m_numVerPerSolidElt; ++iPt) {
1173 ✗ m_duplicateSupportElements->getSubItfEltVerCrd(subElemId, iPt, tmpPt);
1174 ✗ *m_subItfPts[iPt] = Point(tmpPt[0], tmpPt[1], tmpPt[2]);
1175 }
1176
1177 // get the coordinates of the structure element owning this sub-element
1178 ✗ m_interfMesh->getOneElement(eltTypeStruc, elemId, strucElmIdPoint, m_itfPts, 0);
1179
1180 ✗ if( m_skipSmallVolume ){
1181 ✗ const double area = compAreaVolume(m_subItfPts);
1182 ✗ m_totalArea += area;
1183 ✗ if ( area <= m_eps*compAreaVolume(m_itfPts) ) {
1184 ✗ m_skippedArea += std::abs(area);
1185 ✗ continue;
1186 }
1187 }
1188
1189 // get the support elements
1190 ✗ setElemPoint(eltTypeFluid, currentIelGeo, fluidElmPoint, fluidElmIdPoint, vectorIdSupport);
1191
1192 // update the finite element
1193 ✗ m_feVel->updateSubElementFirstDeriv(0, fluidElmPoint, m_subItfPts, m_curvFeStruc);
1194 ✗ m_fePres->updateSubElementFirstDeriv(0, fluidElmPoint, m_subItfPts, m_curvFeStruc);
1195
1196 ✗ m_elemFieldNormal.setValue(itfNormal);
1197
1198 // parameters
1199 ✗ const double penaltyParam = m_viscosity * m_nitschePar / m_feVel->diameter();
1200
1201 // loop over the support elements
1202 ✗ for(std::size_t iSup1 = 0; iSup1 < vectorIdSupport.size(); ++iSup1) {
1203 ✗ const felInt iel1 = vectorIdSupport[iSup1];
1204 ✗ for(std::size_t iSup2 = 0; iSup2 < vectorIdSupport.size(); ++iSup2) {
1205 ✗ const felInt iel2 = vectorIdSupport[iSup2];
1206
1207 // clear elementary matrix
1208 ✗ for(std::size_t j = 0; j < numberOfMatrices(); ++j)
1209 ✗ m_elementMat[j]->zero();
1210
1211 // get the side
1212 ✗ const sideOfInterface sideOfSupportElm = iSup2 == 0 ? LEFT : RIGHT;
1213 ✗ const double side = sideOfSupportElm == LEFT ? 1. : -1.;
1214
1215 // set the extrapolation of velocity
1216 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab ) {
1217 // Get the support element of the old solutions
1218 ✗ std::vector<felInt> ielOld(m_solutionOld.size());
1219 ✗ std::vector<felInt> ielOldDup(m_solutionOld.size());
1220 ✗ std::vector< std::vector<felInt> > oldSupportElement(m_solutionOld.size());
1221 ✗ for(std::size_t i = 0; i < oldSupportElement.size(); ++i) {
1222
1223 ✗ m_supportDofUnknownOld[i][0].getIdElementSupport(currentGlobalIelGeo, oldSupportElement[i]);
1224
1225 // set the id for the old solution
1226 ✗ if( oldSupportElement[i].size() > 1 ) {
1227 // this element is also duplicated for the old solution
1228 ✗ if( sideOfSupportElm == LEFT ) {
1229 ✗ ielOld[i] = oldSupportElement[i][0];
1230 ✗ ielOldDup[i] = oldSupportElement[i][1];
1231 } else {
1232 ✗ ielOld[i] = oldSupportElement[i][1];
1233 ✗ ielOldDup[i] = oldSupportElement[i][0];
1234 }
1235 } else {
1236 // this element is NOT duplicated for the old solution
1237 ✗ ielOld[i] = oldSupportElement[i][0];
1238 ✗ ielOldDup[i] = oldSupportElement[i][0];
1239 }
1240 }
1241
1242 ✗ m_elemFieldAdv.val.clear();
1243 ✗ m_elemFieldAdvDup.val.clear();
1244 ✗ if( FelisceParam::instance().useProjectionForNitscheXFEM ) {
1245 ✗ for(felInt i = 0; i < FelisceParam::instance().orderBdfNS; ++i) {
1246 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, iel1, m_iVelocity, m_ao, dof());
1247 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1248 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, iel2, m_iVelocity, m_ao, dof());
1249 ✗ m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
1250 }
1251 } else {
1252 ✗ for(felInt i = 0; i < FelisceParam::instance().orderBdfNS; ++i) {
1253 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOld[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1254 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1255 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOldDup[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1256 ✗ m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
1257 }
1258 }
1259 }
1260
1261 // The normal is defined external to Omega_1. (Omega_1: right and Omega_2: left)
1262 // side = 1, on the left and -1 on the right.
1263 // The unknown is telling us in which subdomain we are and therefore which normal we have to consider
1264 ✗ if( iSup1 == iSup2 ) {
1265
1266 // Nitsche's penalty term, \gamma \mu / h u_i . v_i
1267 ✗ m_elementMat[0]->phi_i_phi_j(penaltyParam, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1268
1269 // Nitsche's consistency terms, -\sigma(u_i, p_i)n_i . v_i - \sigma(v_i, q_i)n_i . u_i
1270 ✗ computeNitscheSigma_u_v(0.5 * side * m_viscosity, -0.5 * side);
1271 ✗ computeNitscheSigma_v_u(0.5 * side * m_viscosity, -0.5 * side);
1272
1273 // stabilize inflow
1274 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab ){
1275 // (0.5 \rho_f \beta_1 \cdot n_2) (u_1 \cdot v_1)-(0.5 \rho_f \beta_2 \cdot n_2) (u_2 \cdot v_2)
1276 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(0.5 * side * m_density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1277 }
1278 } else {
1279 ✗ m_elementMat[0]->phi_i_phi_j(-penaltyParam, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1280
1281 ✗ computeNitscheSigma_u_v(-0.5 * side * m_viscosity, 0.5 * side);
1282 ✗ computeNitscheSigma_v_u( 0.5 * side * m_viscosity, -0.5 * side);
1283
1284 // stabilize inflow
1285 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab ) {
1286 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(0.25 * side * m_density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1287 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(0.25 * side * m_density, *m_feVel, m_elemFieldAdvDup, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1288 }
1289 }
1290
1291 // add values of elemMat in the global matrix: matrix(0).
1292 ✗ setValueMatrixRHS(vectorIdSupport[iSup1], vectorIdSupport[iSup2], FlagMatrixRHS::only_matrix);
1293 }
1294 }
1295 } // end for on the tip triangle
1296
1297 ✗ if( m_printSkipVolume ) {
1298 ✗ if ( std::abs(m_skippedArea) > m_eps * std::abs(m_totalArea) )
1299 ✗ std::cout << "In assembleMatrixTip, Skipped surface is : " << std::fixed << std::setprecision(16) << m_skippedArea / m_totalArea * 100. << " %" << std::defaultfloat << std::endl;
1300 }
1301 }
1302
1303 // desallocate array use for assemble with petsc.
1304 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
1305 }
1306 }
1307
1308 ✗ void LinearProblemNitscheXFEM::assembleMatrixTipFaces()
1309 {
1310 // ------------------------------------- //
1311 // assemble matrix for tip faces element //
1312 // ------------------------------------- //
1313
1314 // get number of fluid tip elements
1315 ✗ const felInt nbrTipElm = m_duplicateSupportElements->getNumTipElt();
1316 ✗ if( nbrTipElm > 0 ) {
1317
1318 // get set of tip elements
1319 ✗ const auto& idTipElt = m_duplicateSupportElements->getIntersectedTipElt();
1320 bool isOppTipElt;
1321
1322 // get the types of elements
1323 ✗ GeometricMeshRegion::ElementType eltType = meshLocal()->bagElementTypeDomain()[0];
1324 GeometricMeshRegion::ElementType eltTypeOpp;
1325 ✗ const std::size_t numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
1326 ✗ const std::size_t numFacesPerType = dimension() == 2 ?
1327 ✗ GeometricMeshRegion::m_numEdgesPerElt[eltType] :
1328 ✗ GeometricMeshRegion::m_numFacesPerElt[eltType];
1329
1330 // define finite element
1331 ✗ defineFiniteElement(eltType);
1332
1333 // Element matrix and vector initialisation
1334 ✗ initElementArray();
1335
1336 // allocate arrays for assembling the matrix
1337 ✗ allocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
1338
1339 // init variables
1340 ✗ initPerElementType(eltType, FlagMatrixRHS::only_matrix);
1341
1342 // element informations
1343 ✗ std::vector<felInt> ielCurrentSupportElt, ielOppositeSupportElt;
1344 ✗ std::vector<Point*> elemCurrentPoint(numPointPerElt), elemOppositePoint(numPointPerElt);
1345 ✗ std::vector<felInt> elemCurrentIdPoint(numPointPerElt), elemOppositeIdPoint(numPointPerElt);
1346 ✗ std::vector<Point*> facPts(numPointPerElt-1);
1347 ✗ std::vector<double> faceNorm(dimension()), tmpPt(3);
1348 bool isAnInnerFace;
1349 felInt subFaceId;
1350 felInt currentIelGeo;
1351 felInt currentGlobalIelGeo, oppositeGlobalIelGeo;
1352 felInt idCurrentSupElt, idOppositeSupElt;
1353 sideOfInterface sideOfSupportElement, sideOfSupportOpposite;
1354 PetscInt naux;
1355
1356 // mapping from Felisce to Xfem for local ordering of faces or edges
1357 ✗ auto& mapFELtoXFM = m_duplicateSupportElements->getMapFaceFelisceToIntersector( dimension() );
1358
1359 std::array< std::array<short int,3>, 4> mapFACtoVER;
1360 ✗ if( dimension() == 2 )
1361 ✗ mapFACtoVER = { {{0,1,-1},{1,2,-1},{2,0,-1},{-1,-1,-1}} }; // TODO D.C. map from element faces to vertices (should be moved in a class for every geometric element)
1362 else
1363 ✗ mapFACtoVER = { {{0,1,2},{0,3,1},{1,3,2},{0,2,3}} };
1364
1365 // finite element for the opposite element
1366 ✗ const GeoElement *geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
1367 ✗ const felInt typeOfFiniteElement = m_velocity->finiteElementType();
1368 ✗ const RefElement *refEle = geoEle->defineFiniteEle(eltType, typeOfFiniteElement, *meshLocal());
1369 ✗ CurrentFiniteElement oppositeEltFE(*refEle, *geoEle, m_velocity->degreeOfExactness());
1370
1371 // loop over fluid tip elements
1372 ✗ for(felInt iTipElm = 0; iTipElm < nbrTipElm; ++iTipElm) {
1373 // get fluid tip element global and local id, if element doesn't belong to current proc skip it
1374 ✗ currentGlobalIelGeo = m_duplicateSupportElements->getTipEltIdx(iTipElm);
1375 ✗ ISGlobalToLocalMappingApply(mappingElem(), IS_GTOLM_MASK, 1, &currentGlobalIelGeo, &naux, &currentIelGeo);
1376 ✗ if ( currentIelGeo == -1 )
1377 ✗ continue;
1378
1379 ✗ if( m_skipSmallVolume ){
1380 ✗ m_totalArea = 0;
1381 ✗ m_skippedArea = 0;
1382 }
1383
1384 // get informations on the current element
1385 ✗ setElemPoint(eltType, currentIelGeo, elemCurrentPoint, elemCurrentIdPoint, ielCurrentSupportElt);
1386
1387 ✗ m_feVel->updateMeas(0, elemCurrentPoint);
1388
1389 // compute penalty parameter
1390 ✗ const double penaltyParam = m_viscosity * m_nitschePar / m_feVel->diameter();
1391
1392 // loop over all the edges or faces of the current support element
1393 ✗ for(std::size_t iFace = 0; iFace < numFacesPerType; ++iFace) {
1394 // get the neighboor. If it exists, this is an inner edge/face
1395 ✗ oppositeGlobalIelGeo = currentGlobalIelGeo;
1396 ✗ eltTypeOpp = eltType;
1397 ✗ isAnInnerFace = mesh()->getAdjElement(eltTypeOpp, oppositeGlobalIelGeo, iFace);
1398 ✗ if( isAnInnerFace ) {
1399
1400 // get information on the opposite element
1401 ✗ m_supportDofUnknown[0].getIdElementSupport(oppositeGlobalIelGeo, ielOppositeSupportElt);
1402
1403 ✗ mesh()->getOneElement(eltTypeOpp, oppositeGlobalIelGeo, elemOppositeIdPoint, elemOppositePoint, 0);
1404
1405 ✗ if ( idTipElt.find(oppositeGlobalIelGeo) != idTipElt.end() )
1406 ✗ isOppTipElt = true;
1407 else
1408 ✗ isOppTipElt = false;
1409
1410 // compute interior normal
1411 ✗ for (int iPt = 0; iPt < dimension(); ++iPt)
1412 ✗ facPts[iPt] = elemCurrentPoint[mapFACtoVER[iFace][iPt]];
1413
1414 ✗ if( dimension() == 2 )
1415 ✗ normal<2>(facPts, faceNorm);
1416 else
1417 ✗ normal<3>(facPts, faceNorm);
1418
1419 // setting the normal on the face
1420 ✗ m_elemFieldNormal.setValue(faceNorm);
1421
1422 // get number of sub faces
1423 ✗ const felInt nbrSubFaces = m_duplicateSupportElements->getNumSubFacPerFacTipElt(iTipElm, mapFELtoXFM[iFace]);
1424
1425 // loop on subfaces
1426 ✗ for(felInt iSubFace = 0; iSubFace < nbrSubFaces; ++iSubFace){
1427 ✗ subFaceId = m_duplicateSupportElements->getSubTipFacIdxMsh(iTipElm, mapFELtoXFM[iFace], iSubFace);
1428
1429 // get the sign of the subface
1430 ✗ sideOfSupportElement = m_duplicateSupportElements->getSgnSubFace(subFaceId);
1431 ✗ idCurrentSupElt = ielCurrentSupportElt[ sideOfSupportElement == RIGHT ? 1 : 0 ];
1432
1433 // select opposite element support
1434 ✗ sideOfSupportOpposite = m_duplicateSupportElements->getSgnSubOppFace(subFaceId);
1435 ✗ idOppositeSupElt = ielOppositeSupportElt[ sideOfSupportOpposite == RIGHT ? 1 : 0 ];
1436
1437 // fill the vector with the points of the sub-element
1438 ✗ for (int iPt = 0; iPt < dimension(); ++iPt) {
1439 ✗ m_duplicateSupportElements->getSubTipFacVerCrd(subFaceId, iPt, tmpPt);
1440 ✗ *m_subItfPts[iPt] = Point(tmpPt[0], tmpPt[1], tmpPt[2]);
1441 }
1442
1443 ✗ if( m_skipSmallVolume ){
1444 ✗ const double area = compAreaVolume(m_subItfPts);
1445 ✗ m_totalArea += area;
1446 ✗ if ( area <= m_eps*compAreaVolume(facPts) ) {
1447 ✗ m_skippedArea += std::abs(area);
1448 ✗ continue;
1449 }
1450 }
1451
1452 // update the weight and the test function at quadrature points given the subface
1453 ✗ m_feVel->updateSubElementFirstDeriv(0, elemCurrentPoint, m_subItfPts, m_curvFeStruc);
1454
1455 // update the weight and the test function at quadrature points given the subface
1456 ✗ oppositeEltFE.updateSubElementFirstDeriv(0, elemOppositePoint, m_subItfPts, m_curvFeStruc);
1457
1458 // set the extrapolation of velocity
1459 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab ) {
1460 // Get the support element of the old solutions
1461 ✗ std::vector<felInt> ielOld(m_solutionOld.size());
1462 ✗ std::vector<felInt> ielOldOpp(m_solutionOld.size());
1463 ✗ std::vector< std::vector<felInt> > oldSupportElement(m_solutionOld.size());
1464 ✗ for(std::size_t i = 0; i < oldSupportElement.size(); ++i) {
1465
1466 // setting z_1
1467 ✗ m_supportDofUnknownOld[i][0].getIdElementSupport(currentGlobalIelGeo, oldSupportElement[i]);
1468
1469 // set the id for the old solution
1470 ✗ if( oldSupportElement[i].size() > 1 ) {
1471 // this element is also duplicated for the old solution
1472 ✗ if( sideOfSupportElement == LEFT ) {
1473 ✗ ielOld[i] = oldSupportElement[i][0];
1474 } else {
1475 ✗ ielOld[i] = oldSupportElement[i][1];
1476 }
1477 } else {
1478 // this element is NOT duplicated for the old solution
1479 ✗ ielOld[i] = oldSupportElement[i][0];
1480 }
1481
1482 // setting z_2
1483 ✗ m_supportDofUnknownOld[i][0].getIdElementSupport(oppositeGlobalIelGeo, oldSupportElement[i]);
1484
1485 // set the id for the old solution
1486 ✗ if( oldSupportElement[i].size() > 1 ) {
1487 // this element is also duplicated for the old solution
1488 ✗ if( sideOfSupportOpposite == LEFT ) {
1489 ✗ ielOldOpp[i] = oldSupportElement[i][0];
1490 } else {
1491 ✗ ielOldOpp[i] = oldSupportElement[i][1];
1492 }
1493 } else {
1494 // this element is NOT duplicated for the old solution
1495 ✗ ielOldOpp[i] = oldSupportElement[i][0];
1496 }
1497 }
1498
1499 ✗ m_elemFieldAdv.val.clear();
1500 ✗ m_elemFieldAdvDup.val.clear();
1501 ✗ for(felInt i = 0; i < FelisceParam::instance().orderBdfNS; ++i) {
1502 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOld[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1503 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1504 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOldOpp[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1505 ✗ m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
1506 }
1507 }
1508
1509 // -------------------------------------------------------- //
1510 // current element for phi_i and psi_j //
1511 // -------------------------------------------------------- //
1512 ✗ m_elementMat[0]->zero();
1513
1514 // Nitsche's penalty term, \gamma \mu / h u . v
1515 ✗ m_elementMat[0]->psi_j_phi_i(penaltyParam, *m_feVel, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1516
1517 //Nitsche's consistency terms, -\sigma(u , p) n . v
1518 ✗ if( m_useSymmetricStress )
1519 ✗ m_elementMat[0]->eps_psi_j_dot_vec_phi_i(m_viscosity, m_elemFieldNormal, *m_feVel, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1520 else
1521 ✗ m_elementMat[0]->grad_psi_j_dot_vec_phi_i(0.5 * m_viscosity, m_elemFieldNormal, *m_feVel, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1522 // same finite element for the pressure and the velocity
1523 ✗ for(std::size_t idim = 0; idim < m_velocity->numComponent(); ++idim)
1524 ✗ m_elementMat[0]->psi_j_phi_i(-0.5 * m_elemFieldNormal.val(idim,0), *m_feVel, *m_feVel, m_velBlock+idim, m_preBlock, m_pressure->numComponent());
1525
1526 // Nitsche's symmetry terms, - u . \sigma(v , - q )n,
1527 // !! the equation for q is multiplied by - !!
1528 ✗ if( m_useSymmetricStress )
1529 ✗ m_elementMat[0]->psi_j_eps_phi_i_dot_vec(m_viscosity, m_elemFieldNormal, *m_feVel, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1530 else
1531 ✗ m_elementMat[0]->psi_j_grad_phi_i_dot_vec(0.5 * m_viscosity, m_elemFieldNormal, *m_feVel, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1532 // same finite element for the pressure and velocity
1533 ✗ for(std::size_t idim = 0; idim < m_velocity->numComponent(); ++idim)
1534 ✗ m_elementMat[0]->psi_j_phi_i(-0.5 * m_elemFieldNormal.val(idim,0), *m_feVel, *m_feVel, m_preBlock, m_velBlock+idim, m_pressure->numComponent());
1535
1536 // stabilize inflow (0.5 \rho_f \beta_1 \cdot n) (u_1 \cdot v_1)
1537 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab )
1538 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(0.5 * m_density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1539
1540 ✗ if ( isOppTipElt )
1541 ✗ *m_elementMat[0] *= 0.5;
1542
1543 ✗ setValueMatrixRHS(idCurrentSupElt, idCurrentSupElt, FlagMatrixRHS::only_matrix);
1544
1545
1546 // -------------------------------------------------------- //
1547 // current element for phi_i and opposite element for psi_j //
1548 // -------------------------------------------------------- //
1549 ✗ m_elementMat[0]->zero();
1550
1551 // Nitsche's penalty term, \gamma \mu / h u . v
1552 ✗ m_elementMat[0]->psi_j_phi_i(-penaltyParam, *m_feVel, oppositeEltFE, m_velBlock, m_velBlock, m_velocity->numComponent());
1553
1554 // Nitsche's consistency terms, -\sigma(u , p) n . v
1555 ✗ if( m_useSymmetricStress )
1556 ✗ m_elementMat[0]->eps_psi_j_dot_vec_phi_i(m_viscosity, m_elemFieldNormal, *m_feVel, oppositeEltFE, m_velBlock, m_velBlock, m_velocity->numComponent());
1557 else
1558 ✗ m_elementMat[0]->grad_psi_j_dot_vec_phi_i(0.5 * m_viscosity, m_elemFieldNormal, *m_feVel, oppositeEltFE, m_velBlock, m_velBlock, m_velocity->numComponent());
1559 // same finite element for the pressure and the velocity
1560 ✗ for(std::size_t idim = 0; idim < m_velocity->numComponent(); ++idim)
1561 ✗ m_elementMat[0]->psi_j_phi_i(-0.5 * m_elemFieldNormal.val(idim,0), *m_feVel, oppositeEltFE, m_velBlock+idim, m_preBlock, m_pressure->numComponent());
1562
1563 // Nitsche's symmetry terms, - u . \sigma(v , - q )n,
1564 // !! the equation for q is multiplied by - !!
1565 ✗ if( m_useSymmetricStress )
1566 ✗ m_elementMat[0]->psi_j_eps_phi_i_dot_vec(-m_viscosity, m_elemFieldNormal, *m_feVel, oppositeEltFE, m_velBlock, m_velBlock, m_velocity->numComponent());
1567 else
1568 ✗ m_elementMat[0]->psi_j_grad_phi_i_dot_vec(-0.5 * m_viscosity, m_elemFieldNormal, *m_feVel, oppositeEltFE, m_velBlock, m_velBlock, m_velocity->numComponent());
1569 // same finite element for the pressure and velocity
1570 ✗ for(std::size_t idim = 0; idim < m_velocity->numComponent(); ++idim)
1571 ✗ m_elementMat[0]->psi_j_phi_i(0.5 * m_elemFieldNormal.val(idim,0), *m_feVel, oppositeEltFE, m_preBlock, m_velBlock+idim, m_pressure->numComponent());
1572
1573 // stabilize inflow +(0.25 \rho_f \beta_1 \cdot n) (u_2 \cdot v_1) + (0.25 \rho_f \beta_2 \cdot n) (u_2 \cdot v_1)
1574 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab ) {
1575 ✗ m_elementMat[0]->f_dot_n_psi_j_phi_i(-0.25 * m_density, *m_feVel, oppositeEltFE, m_elemFieldAdv, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1576 ✗ m_elementMat[0]->f_dot_n_psi_j_phi_i(-0.25 * m_density, *m_feVel, oppositeEltFE, m_elemFieldAdvDup, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1577 }
1578
1579 ✗ if ( isOppTipElt )
1580 ✗ *m_elementMat[0] *= 0.5;
1581
1582 ✗ setValueMatrixRHS(idCurrentSupElt, idOppositeSupElt, FlagMatrixRHS::only_matrix);
1583
1584
1585 // -------------------------------------------------------- //
1586 // opposite element for phi_i and current element for psi_j //
1587 // -------------------------------------------------------- //
1588 ✗ m_elementMat[0]->zero();
1589
1590 // Nitsche's penalty term, \gamma \mu / h u . v
1591 ✗ m_elementMat[0]->psi_j_phi_i(-penaltyParam, oppositeEltFE, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1592
1593 // Nitsche's consistency terms, -\sigma(u , p) n . v
1594 ✗ if( m_useSymmetricStress )
1595 ✗ m_elementMat[0]->eps_psi_j_dot_vec_phi_i(-m_viscosity, m_elemFieldNormal, oppositeEltFE, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1596 else
1597 ✗ m_elementMat[0]->grad_psi_j_dot_vec_phi_i(-0.5 * m_viscosity, m_elemFieldNormal, oppositeEltFE, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1598 // same finite element for the pressure and the velocity
1599 ✗ for(std::size_t idim = 0; idim < m_velocity->numComponent(); ++idim)
1600 ✗ m_elementMat[0]->psi_j_phi_i(0.5 * m_elemFieldNormal.val(idim,0), oppositeEltFE, *m_feVel, m_velBlock+idim, m_preBlock, m_pressure->numComponent());
1601
1602 // Nitsche's symmetry terms, - u . \sigma(v , - q )n,
1603 // !! the equation for q is multiplied by - !!
1604 ✗ if( m_useSymmetricStress )
1605 ✗ m_elementMat[0]->psi_j_eps_phi_i_dot_vec(m_viscosity, m_elemFieldNormal, oppositeEltFE, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1606 else
1607 ✗ m_elementMat[0]->psi_j_grad_phi_i_dot_vec(0.5 * m_viscosity, m_elemFieldNormal, oppositeEltFE, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
1608 // same finite element for the pressure and velocity
1609 ✗ for(std::size_t idim = 0; idim < m_velocity->numComponent(); ++idim)
1610 ✗ m_elementMat[0]->psi_j_phi_i(-0.5 * m_elemFieldNormal.val(idim,0), oppositeEltFE, *m_feVel, m_preBlock, m_velBlock+idim, m_pressure->numComponent());
1611
1612 // stabilize inflow. - (0.25 \rho_f \beta_1 \cdot n) (u_1 \cdot v_2) - (0.25 \rho_f \beta_2 \cdot n) (u_1 \cdot v_2)
1613 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab ) {
1614 ✗ m_elementMat[0]->f_dot_n_psi_j_phi_i(0.25 * m_density, oppositeEltFE, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1615 ✗ m_elementMat[0]->f_dot_n_psi_j_phi_i(0.25 * m_density, oppositeEltFE, *m_feVel, m_elemFieldAdvDup, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1616 }
1617
1618 ✗ if ( isOppTipElt )
1619 ✗ *m_elementMat[0] *= 0.5;
1620
1621 ✗ setValueMatrixRHS(idOppositeSupElt, idCurrentSupElt, FlagMatrixRHS::only_matrix);
1622
1623
1624 // -------------------------------------------------------- //
1625 // opposite element for phi_i and psi_j //
1626 // -------------------------------------------------------- //
1627 ✗ m_elementMat[0]->zero();
1628
1629 // Nitsche's penalty term, \gamma \mu / h u . v
1630 ✗ m_elementMat[0]->psi_j_phi_i(penaltyParam, oppositeEltFE, oppositeEltFE, m_velBlock, m_velBlock, m_velocity->numComponent());
1631
1632 // Nitsche's consistency terms, -\sigma(u , p) n . v
1633 ✗ if( m_useSymmetricStress )
1634 ✗ m_elementMat[0]->eps_psi_j_dot_vec_phi_i(-m_viscosity, m_elemFieldNormal, oppositeEltFE, oppositeEltFE, m_velBlock, m_velBlock, m_velocity->numComponent());
1635 else
1636 ✗ m_elementMat[0]->grad_psi_j_dot_vec_phi_i(-0.5 * m_viscosity, m_elemFieldNormal, oppositeEltFE, oppositeEltFE, m_velBlock, m_velBlock, m_velocity->numComponent());
1637 // same finite element for the pressure and the velocity
1638 ✗ for(std::size_t idim = 0; idim < m_velocity->numComponent(); ++idim)
1639 ✗ m_elementMat[0]->psi_j_phi_i(0.5 * m_elemFieldNormal.val(idim,0), oppositeEltFE, oppositeEltFE, m_velBlock+idim, m_preBlock, m_pressure->numComponent());
1640
1641 // Nitsche's symmetry terms, - u . \sigma(v , - q )n,
1642 // !! the equation for q is multiplied by - !!
1643 ✗ if( m_useSymmetricStress )
1644 ✗ m_elementMat[0]->psi_j_eps_phi_i_dot_vec(-m_viscosity, m_elemFieldNormal, oppositeEltFE, oppositeEltFE, m_velBlock, m_velBlock, m_velocity->numComponent());
1645 else
1646 ✗ m_elementMat[0]->psi_j_grad_phi_i_dot_vec(-0.5 * m_viscosity, m_elemFieldNormal, oppositeEltFE, oppositeEltFE, m_velBlock, m_velBlock, m_velocity->numComponent());
1647 // same finite element for the pressure and velocity
1648 ✗ for(std::size_t idim = 0; idim < m_velocity->numComponent(); ++idim)
1649 ✗ m_elementMat[0]->psi_j_phi_i(0.5 * m_elemFieldNormal.val(idim,0), oppositeEltFE, oppositeEltFE, m_preBlock, m_velBlock+idim, m_pressure->numComponent());
1650
1651 // stabilize inflow + (0.5 \rho_f \beta_2 \cdot n) (u_2 \cdot v_2)
1652 ✗ if( m_useNSEquation == 1 && m_useInterfaceStab )
1653 ✗ m_elementMat[0]->f_dot_n_phi_i_phi_j(-0.5 * m_density, oppositeEltFE, m_elemFieldAdvDup, m_elemFieldNormal, m_velBlock, m_velBlock, m_velocity->numComponent(), 0);
1654
1655 ✗ if ( isOppTipElt )
1656 ✗ *m_elementMat[0] *= 0.5;
1657
1658 ✗ setValueMatrixRHS(idOppositeSupElt, idOppositeSupElt, FlagMatrixRHS::only_matrix);
1659 }
1660 }
1661 }
1662
1663 ✗ if( m_printSkipVolume ) {
1664 ✗ if ( std::abs(m_skippedArea) > m_eps * std::abs(m_totalArea) )
1665 ✗ std::cout << "In assembleMatrixTipFaces, Skipped surface is : " << std::fixed << std::setprecision(16) << m_skippedArea / m_totalArea * 100. << " %" << std::defaultfloat << std::endl;
1666 }
1667 }
1668
1669 // desallocate array use for assemble with petsc.
1670 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
1671 }
1672 }
1673
1674 ✗ void LinearProblemNitscheXFEM::assembleMatrixGhostPenalty()
1675 {
1676 // ------------- //
1677 // Ghost penalty //
1678 // ------------- //
1679 // get the types of elements
1680 ✗ GeometricMeshRegion::ElementType eltTypeOpp, eltType = meshLocal()->bagElementTypeDomain()[0]; // list of domain element types in the mesh
1681 ✗ const felInt numFacePerType = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numFace(); //if we are in 3d we integrate on the faces
1682 ✗ const felInt numEdgePerType = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numEdge(); //if we are in 2d we integrate on the edges
1683 ✗ std::vector<felInt> idOfFacesCurrent( dimension() == 2 ? numEdgePerType : numFacePerType, -1);
1684 ✗ std::vector<felInt> idOfFacesOpposite(dimension() == 2 ? numEdgePerType : numFacePerType, -1);
1685
1686 // define finite element
1687 ✗ defineFiniteElement(eltType);
1688 ✗ defineCurrentFiniteElementWithBd(eltType);
1689
1690 // Element matrix and vector initialisation
1691 ✗ initElementArray();
1692
1693 // allocate arrays for assembling the matrix
1694 ✗ allocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
1695
1696 // init variables
1697 ✗ initPerElementType(eltType, FlagMatrixRHS::only_matrix);
1698 ✗ m_feVelWithBd = m_listCurrentFiniteElementWithBd[m_iVelocity];
1699
1700 // get global supportDofMesh informations
1701 ✗ std::vector<felInt>& iEle = m_supportDofUnknown[0].iEle();
1702 ✗ std::vector<felInt>& iSupportDof = m_supportDofUnknown[0].iSupportDof();
1703
1704 // element informations
1705 bool isAnInnerFace;
1706 ✗ std::vector<felInt> ielCurrentSupportElt, ielOppositeSupportElt;
1707 ✗ std::vector<Point*> elemCurrentPoint(m_numVerPerFluidElt), elemOppositePoint(m_numVerPerFluidElt);
1708 ✗ std::vector<felInt> elemCurrentIdPoint(m_numVerPerFluidElt), elemOppositeIdPoint(m_numVerPerFluidElt);
1709 felInt currentIelGeo;
1710 felInt currentGlobalIelGeo, oppositeGlobalIelGeo;
1711 felInt idCurrentSupElt, idOppositeSupElt;
1712 PetscInt naux;
1713
1714 // finite element with boundary for the opposite element
1715 ✗ const GeoElement *geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
1716 ✗ const felInt typeOfFiniteElement = m_velocity->finiteElementType();
1717 ✗ const RefElement *refEle = geoEle->defineFiniteEle(eltType, typeOfFiniteElement, *meshLocal());
1718 ✗ CurrentFiniteElementWithBd oppositeEltFEWithBd(*refEle, *geoEle, m_velocity->degreeOfExactness(), m_velocity->degreeOfExactness());
1719
1720 // loop over all the intersected elements
1721 ✗ const auto& idInterElt = m_duplicateSupportElements->getIntersectedMshElt();
1722 ✗ for(auto ielIt = idInterElt.begin(); ielIt != idInterElt.end() && ielIt->first < mesh()->numElements(eltType); ++ielIt) {
1723 // get the geometric id of the element and the ids of its support elements
1724 ✗ currentGlobalIelGeo = ielIt->first;
1725 ✗ ISGlobalToLocalMappingApply(mappingElem(), IS_GTOLM_MASK, 1, &currentGlobalIelGeo, &naux, &currentIelGeo);
1726
1727 // jump to the next intersected element if this element doesn't belong to this processor
1728 ✗ if( currentIelGeo == -1)
1729 ✗ continue;
1730
1731 // get all the faces of the element
1732 ✗ if( dimension() == 2 )
1733 ✗ mesh()->getAllEdgeOfElement(currentGlobalIelGeo, idOfFacesCurrent);
1734 else
1735 ✗ mesh()->getAllFaceOfElement(currentGlobalIelGeo, idOfFacesCurrent);
1736
1737 // get informations on the current element
1738 ✗ setElemPoint(eltType, currentIelGeo, elemCurrentPoint, elemCurrentIdPoint, ielCurrentSupportElt);
1739
1740 // loop over all the edges or faces of the current support element
1741 ✗ for(std::size_t iFace = 0; iFace < idOfFacesCurrent.size(); ++iFace) {
1742
1743 // get the neighboor. If it exists, this is an inner edge
1744 ✗ oppositeGlobalIelGeo = currentGlobalIelGeo;
1745 ✗ eltTypeOpp = eltType;
1746 ✗ isAnInnerFace = mesh()->getAdjElement(eltTypeOpp, oppositeGlobalIelGeo, iFace);
1747
1748 ✗ if( isAnInnerFace ) {
1749
1750 // get information on the opposite element
1751 ✗ m_supportDofUnknown[0].getIdElementSupport(oppositeGlobalIelGeo, ielOppositeSupportElt);
1752
1753 ✗ mesh()->getOneElement(eltTypeOpp, oppositeGlobalIelGeo, elemOppositeIdPoint, elemOppositePoint, 0);
1754
1755 // get number of dof in support element
1756 ✗ const std::size_t sizeCurrent = m_supportDofUnknown[0].getNumSupportDof(ielCurrentSupportElt[0]);
1757 ✗ const std::size_t sizeOpposite = m_supportDofUnknown[0].getNumSupportDof(ielOppositeSupportElt[0]);
1758
1759 // for each support elements
1760 ✗ for(std::size_t iSup = 0; iSup < ielCurrentSupportElt.size(); ++iSup) {
1761
1762 // get current support element
1763 ✗ idCurrentSupElt = ielCurrentSupportElt[iSup];
1764
1765 // the other element is also duplicated, take the one on the same side
1766 ✗ if( ielOppositeSupportElt.size() > 1 )
1767 ✗ idOppositeSupElt = ielOppositeSupportElt[iSup];
1768 else
1769 ✗ idOppositeSupElt = ielOppositeSupportElt[0];
1770
1771 // compute dof in common
1772 ✗ std::size_t countSupportDofInCommon = 0;
1773 ✗ for(std::size_t iSup1=0; iSup1<sizeCurrent; ++iSup1) {
1774 ✗ for(std::size_t iSup2=0; iSup2<sizeOpposite; ++iSup2)
1775 ✗ if( iSupportDof[iEle[idOppositeSupElt] + iSup2] == iSupportDof[iEle[idCurrentSupElt] + iSup1] )
1776 ✗ ++countSupportDofInCommon;
1777 }
1778
1779 // if the two support elements share an edge/face
1780 ✗ bool computeGhostPenalty = false;
1781 ✗ if( countSupportDofInCommon == sizeCurrent-1)
1782 ✗ computeGhostPenalty = true;
1783
1784
1785 // add the contribution of this edge/face
1786 ✗ if( computeGhostPenalty ) {
1787
1788
1789 // if opposite element is not duplicated compute stabilization both side
1790 ✗ bool computeBothSide = false;
1791 ✗ if( ielOppositeSupportElt.size() > 1 )
1792 ✗ computeBothSide = true;
1793
1794 // update the finite elements
1795 ✗ m_feVelWithBd->updateFirstDeriv(0, elemCurrentPoint);
1796 ✗ m_feVelWithBd->updateBdMeasNormal();
1797 ✗ oppositeEltFEWithBd.updateFirstDeriv(0, elemOppositePoint);
1798 ✗ oppositeEltFEWithBd.updateBdMeasNormal();
1799
1800 // find the local id of the edge in the opposite element
1801 ✗ if( dimension() == 2 )
1802 ✗ mesh()->getAllEdgeOfElement(eltType, oppositeGlobalIelGeo, idOfFacesOpposite);
1803 else
1804 ✗ mesh()->getAllFaceOfElement(eltType, oppositeGlobalIelGeo, idOfFacesOpposite);
1805
1806 ✗ felInt localIdFaceOpposite = -1;
1807 ✗ for(std::size_t jface = 0; jface < idOfFacesOpposite.size(); ++jface) {
1808 ✗ if( idOfFacesCurrent[iFace] == idOfFacesOpposite[jface] ) {
1809 ✗ localIdFaceOpposite = jface;
1810 }
1811 }
1812
1813 // set coefficient of the ghost penalty
1814 ✗ const double gammag = FelisceParam::instance().coefGhostPenalty * m_viscosity * m_feVelWithBd->bdEle(iFace).diameter();
1815
1816 // ----------------------------------- //
1817 // current element for phi_i and phi_j //
1818 // ----------------------------------- //
1819 ✗ m_elementMat[0]->zero();
1820 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammag, *m_feVelWithBd, *m_feVelWithBd, iFace, iFace, m_velBlock, m_velBlock, m_velocity->numComponent());
1821 ✗ setValueMatrixRHS(idCurrentSupElt, idCurrentSupElt, FlagMatrixRHS::only_matrix);
1822
1823 // -------------------------------------------------------- //
1824 // current element for phi_i and opposite element for phi_j //
1825 // -------------------------------------------------------- //
1826 ✗ m_elementMat[0]->zero();
1827 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammag, oppositeEltFEWithBd, *m_feVelWithBd, localIdFaceOpposite, iFace, m_velBlock, m_velBlock, m_velocity->numComponent());
1828 ✗ setValueMatrixRHS(idCurrentSupElt, idOppositeSupElt, FlagMatrixRHS::only_matrix);
1829
1830 // if the opposite element is not duplicated, it will not be handle in this loop
1831 // compute all the term in this case
1832 ✗ if( computeBothSide ) {
1833 // -------------------------------------------------------- //
1834 // opposite element for phi_i and current element for phi_j //
1835 // -------------------------------------------------------- //
1836 ✗ m_elementMat[0]->zero();
1837 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammag, *m_feVelWithBd, oppositeEltFEWithBd, iFace, localIdFaceOpposite, m_velBlock, m_velBlock, m_velocity->numComponent());
1838 ✗ setValueMatrixRHS(idOppositeSupElt, idCurrentSupElt, FlagMatrixRHS::only_matrix);
1839
1840 // ------------------------------------ //
1841 // opposite element for phi_i and phi_j //
1842 // ------------------------------------ //
1843 ✗ m_elementMat[0]->zero();
1844 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammag, oppositeEltFEWithBd, oppositeEltFEWithBd, localIdFaceOpposite, localIdFaceOpposite, m_velBlock, m_velBlock, m_velocity->numComponent());
1845 ✗ setValueMatrixRHS(idOppositeSupElt, idOppositeSupElt, FlagMatrixRHS::only_matrix);
1846 }
1847 }
1848 }
1849 }
1850 }
1851 }
1852
1853 // desallocate array use for assemble with petsc.
1854 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
1855 }
1856
1857 ✗ void LinearProblemNitscheXFEM::assembleMatrixFaceOriented()
1858 {
1859 // Use Face-oriented stabilization
1860 // E. Burman, M. A. Fernandez and P. Hansbo, Continuous interior penalty finite element method for
1861 // Oseen's equations, 2006.
1862 ✗ std::pair<bool, GeometricMeshRegion::ElementType> adjElt;
1863
1864 GeometricMeshRegion::ElementType eltType; // Type of element
1865 GeometricMeshRegion::ElementType eltTypeOpp; // Type of element
1866
1867 bool computeFaceOriented; // true if we need to compute the stabilization on this face
1868
1869 ✗ felInt ielCurrentLocalGeo = 0; // local geometric id of the current element
1870 ✗ felInt ielCurrentGlobalGeo = 0; // global geometric id of the current element
1871 ✗ felInt ielOppositeGlobalGeo = 0; // global geometric id of the opposite element
1872 felInt numFacesPerType; // number of faces of an element of a given type
1873 felInt numEltPerLabel; // number of element by reference
1874 felInt numPointPerElt; // number of vertices by element
1875 felInt currentSupport; // id of the current support element
1876 felInt oppositeSupport; // id of the opposite support element
1877 ✗ std::vector<felInt> ielCurrentGlobalSup; // local support element id of the current element
1878 ✗ std::vector<felInt> ielOppositeGlobalSup; // local support element id of the opposite element
1879
1880 ✗ std::vector<felInt> idOfFacesCurrent; // ids of the current element edges
1881 ✗ std::vector<felInt> idOfFacesOpposite; // ids of the opposite element edges
1882 ✗ std::vector<felInt> currentElemIdPoint; // ids of the vertices of the current element
1883 ✗ std::vector<felInt> oppositeElemIdPoint; // ids of the vertices of the opposite element
1884 ✗ std::vector<Point*> currentElemPoint; // point coordinates of the current element vertices
1885 ✗ std::vector<Point*> oppositeElemPoint; // point coordinates of the opposite element vertices
1886 ✗ std::vector<felInt> countElement; // count the element by type
1887
1888 ✗ ElementField elemFieldAdvFace;
1889 ✗ ElementField elemFieldAdvFaceTmp;
1890
1891 ✗ const std::vector<felInt>& iSupportDof = m_supportDofUnknown[0].iSupportDof();
1892 ✗ const std::vector<felInt>& iEle = m_supportDofUnknown[0].iEle();
1893
1894 // resize vector
1895 ✗ countElement.resize(GeometricMeshRegion::m_numTypesOfElement, 0);
1896
1897 // volume bag
1898 ✗ const std::vector<ElementType>& bagElementTypeDomain = meshLocal()->bagElementTypeDomain();
1899
1900 // finite element with bd for the opposite element
1901 ✗ std::map<GeometricMeshRegion::ElementType, CurrentFiniteElementWithBd*> oppositeEltFEWithBdPre;
1902 ✗ std::map<GeometricMeshRegion::ElementType, CurrentFiniteElementWithBd*> oppositeEltFEWithBdVel;
1903
1904 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
1905 ✗ eltType = bagElementTypeDomain[i];
1906 ✗ const GeoElement *geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
1907
1908 ✗ felInt typeOfFEVel = m_velocity->finiteElementType();
1909 ✗ const RefElement *refEleVel = geoEle->defineFiniteEle(eltType, typeOfFEVel, *mesh());
1910 ✗ oppositeEltFEWithBdVel[eltType] = new CurrentFiniteElementWithBd(*refEleVel, *geoEle, m_velocity->degreeOfExactness(), m_velocity->degreeOfExactness());
1911
1912 ✗ felInt typeOfFEPre = m_pressure->finiteElementType();
1913 ✗ const RefElement *refElePre = geoEle->defineFiniteEle(eltType, typeOfFEPre, *mesh());
1914 ✗ oppositeEltFEWithBdPre[eltType] = new CurrentFiniteElementWithBd(*refElePre, *geoEle, m_pressure->degreeOfExactness(), m_pressure->degreeOfExactness());
1915 }
1916
1917 // First loop on geometric type
1918 ✗ for (std::size_t itype = 0; itype < bagElementTypeDomain.size(); ++itype) {
1919 // get element type and number of vertices and "faces"
1920 ✗ eltType = bagElementTypeDomain[itype];
1921 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
1922 ✗ numFacesPerType = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numBdEle();
1923
1924 // resize of vectors
1925 ✗ idOfFacesCurrent.resize(numFacesPerType, -1);
1926 ✗ idOfFacesOpposite.resize(numFacesPerType, -1);
1927 ✗ currentElemIdPoint.resize(numPointPerElt, -1);
1928 ✗ oppositeElemIdPoint.resize(numPointPerElt, -1);
1929 ✗ currentElemPoint.resize(numPointPerElt, nullptr);
1930 ✗ oppositeElemPoint.resize(numPointPerElt, nullptr);
1931
1932 // define finite element
1933 ✗ defineFiniteElement(eltType);
1934 ✗ defineCurrentFiniteElementWithBd(eltType);
1935
1936 // Element matrix and vector initialisation
1937 ✗ initElementArray();
1938
1939 // allocate arrays for assembling the matrix
1940 ✗ allocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
1941
1942 // init variables
1943 ✗ initPerElementType(eltType, FlagMatrixRHS::only_matrix);
1944 ✗ m_feVelWithBd = m_listCurrentFiniteElementWithBd[m_iVelocity];
1945 ✗ m_fePreWithBd = m_listCurrentFiniteElementWithBd[m_iPressure];
1946
1947 // second loop on region of the mesh
1948 ✗ for(auto itRef = meshLocal()->intRefToBegEndMaps[eltType].begin(); itRef != meshLocal()->intRefToBegEndMaps[eltType].end(); itRef++) {
1949 // get the number of element with this reference
1950 ✗ numEltPerLabel = itRef->second.second;
1951
1952 // set the current label of LinearProblem to the current one of this loop
1953 ✗ initPerDomain(itRef->first, FlagMatrixRHS::only_matrix);
1954
1955 // third loop on element
1956 ✗ for(felInt iel = 0; iel < numEltPerLabel; ++iel) {
1957 // get the global id of the current geometric element
1958 ✗ ISLocalToGlobalMappingApply(mappingElem(), 1, &ielCurrentLocalGeo, &ielCurrentGlobalGeo);
1959
1960 // get all the faces of the element
1961 ✗ if( dimension() == 2 )
1962 ✗ mesh()->getAllEdgeOfElement(ielCurrentGlobalGeo, idOfFacesCurrent);
1963 else
1964 ✗ mesh()->getAllFaceOfElement(ielCurrentGlobalGeo, idOfFacesCurrent);
1965
1966 // get informations on the current element
1967 ✗ setElemPoint(eltType, countElement[eltType], currentElemPoint, currentElemIdPoint, ielCurrentGlobalSup);
1968
1969 // update the finite elements
1970 ✗ m_feVelWithBd->updateFirstDeriv(0, currentElemPoint);
1971 ✗ m_feVelWithBd->updateBdMeasNormal();
1972 ✗ m_fePreWithBd->updateFirstDeriv(0, currentElemPoint);
1973 ✗ m_fePreWithBd->updateBdMeasNormal();
1974
1975 // loop over the support element of the current geometric element
1976 ✗ for(std::size_t ielSup = 0; ielSup < ielCurrentGlobalSup.size(); ++ielSup) {
1977 ✗ currentSupport = ielCurrentGlobalSup[ielSup];
1978
1979 // get the old support dof
1980 ✗ std::vector< std::vector<felInt> > oldSupportElement(FelisceParam::instance().orderBdfNS);
1981 ✗ std::vector<felInt> ielOld(FelisceParam::instance().orderBdfNS);
1982
1983 ✗ for(felInt iextrap = 0; iextrap < FelisceParam::instance().orderBdfNS; ++iextrap) {
1984 // get the old support element
1985 ✗ m_supportDofUnknownOld[iextrap][0].getIdElementSupport(ielCurrentGlobalGeo, oldSupportElement[iextrap]);
1986
1987 // set the id for the old solution
1988 ✗ if( oldSupportElement[iextrap].size() > 1 )
1989 ✗ ielOld[iextrap] = oldSupportElement[iextrap][ielSup];
1990 else
1991 ✗ ielOld[iextrap] = oldSupportElement[iextrap][0];
1992 }
1993
1994 // compute coefficients
1995 ✗ m_elemFieldAdv.val.clear();
1996 ✗ if( FelisceParam::instance().useProjectionForNitscheXFEM ) {
1997 ✗ for(felInt iextrap = 0; iextrap < FelisceParam::instance().orderBdfNS; ++iextrap) {
1998 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[iextrap], *m_feVel, currentSupport, m_iVelocity, m_ao, dof());
1999 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
2000 }
2001 } else {
2002 ✗ for(felInt iextrap = 0; iextrap < FelisceParam::instance().orderBdfNS; ++iextrap) {
2003 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[iextrap], *m_feVel, ielOld[iextrap], m_iVelocity, m_aoOld[iextrap], m_dofOld[iextrap]);
2004 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
2005 }
2006 }
2007
2008 ✗ double normVelExtrap = 0;
2009 ✗ for(std::size_t jcol = 0; jcol < m_elemFieldAdv.val.size2(); ++jcol)
2010 ✗ normVelExtrap = std::max(normVelExtrap, norm_2(column(m_elemFieldAdv.val, jcol)));
2011
2012 ✗ const double hK = m_feVelWithBd->diameter();
2013 ✗ const double ReK = m_density * normVelExtrap * hK / m_viscosity;
2014 ✗ const double coeffVel = FelisceParam::instance().stabFOVel * std::min(1., ReK) * hK * hK;
2015
2016 ✗ double gammaP = 0.;
2017 ✗ if( m_useNSEquation == 0 ) {
2018 ✗ gammaP = FelisceParam::instance().stabFOPre * hK * hK * hK / m_viscosity;
2019 } else {
2020 ✗ double xiPre = ReK < 1 ? m_density * hK * hK * hK / m_viscosity : hK * hK / normVelExtrap;
2021 ✗ gammaP = FelisceParam::instance().stabFOPre * xiPre;
2022 }
2023
2024 // loop over all the faces of the current geometric element
2025 ✗ for(std::size_t iFace = 0; iFace < idOfFacesCurrent.size(); ++iFace) {
2026 // check if this face is an inner face or a boundary
2027 ✗ ielOppositeGlobalGeo = ielCurrentGlobalGeo;
2028
2029 // first = true if there is an adjacent element
2030 // second = the ElementType of the adjacent element if first is true
2031 ✗ adjElt = mesh()->getAdjElement(ielOppositeGlobalGeo, iFace);
2032
2033 // compute the stabilization only on inner faces
2034 ✗ if( adjElt.first ) {
2035 // get the type of the opposite element
2036 ✗ eltTypeOpp = adjElt.second;
2037
2038 // find the local id of the edge in the opposite element
2039 ✗ if( dimension() == 2 )
2040 ✗ mesh()->getAllEdgeOfElement(ielOppositeGlobalGeo, idOfFacesOpposite);
2041 else
2042 ✗ mesh()->getAllFaceOfElement(ielOppositeGlobalGeo, idOfFacesOpposite);
2043
2044 ✗ felInt localIdFaceOpposite = -1;
2045 ✗ for(std::size_t jface = 0; jface < idOfFacesOpposite.size(); ++jface) {
2046 ✗ if( idOfFacesCurrent[iFace] == idOfFacesOpposite[jface] )
2047 ✗ localIdFaceOpposite = jface;
2048 }
2049
2050 // get information on the opposite element
2051 // same as the function "setElemPoint" but here ielOppositeGlobalGeo is global
2052 // we assume that the supportDofMesh is the same for all unknown (like in setElemPoint)
2053 ✗ m_supportDofUnknown[0].getIdElementSupport(ielOppositeGlobalGeo, ielOppositeGlobalSup);
2054
2055 ✗ mesh()->getOneElement(eltTypeOpp, ielOppositeGlobalGeo, oppositeElemIdPoint, oppositeElemPoint, 0);
2056
2057 // get number of dof in support element
2058 ✗ const std::size_t sizeCurrent = m_supportDofUnknown[0].getNumSupportDof(ielCurrentGlobalSup[0]);
2059 ✗ const std::size_t sizeOpposite = m_supportDofUnknown[0].getNumSupportDof(ielOppositeGlobalSup[0]);
2060
2061 // compare the number of support element
2062 ✗ computeFaceOriented = false;
2063 ✗ if( ielCurrentGlobalSup.size() == ielOppositeGlobalSup.size() ) {
2064
2065 // both are duplicated or both are not. Take the element on the same side
2066 ✗ oppositeSupport = ielOppositeGlobalSup[ielSup];
2067
2068 // compute dof in common
2069 ✗ std::size_t countSupportDofInCommon = 0;
2070 ✗ for(std::size_t iSup1=0; iSup1<sizeCurrent; ++iSup1) {
2071 ✗ for(std::size_t iSup2=0; iSup2<sizeOpposite; ++iSup2)
2072 ✗ if( iSupportDof[iEle[oppositeSupport] + iSup2] == iSupportDof[iEle[currentSupport] + iSup1] )
2073 ✗ ++countSupportDofInCommon;
2074 }
2075
2076 // if the two support elements share an edge/face
2077 ✗ if( countSupportDofInCommon == sizeCurrent-1 )
2078 ✗ computeFaceOriented = true;
2079 }
2080 else {
2081
2082 // One is duplicated, the other is not
2083 ✗ for(std::size_t jelSup=0; jelSup < ielOppositeGlobalSup.size() && !computeFaceOriented; ++jelSup) {
2084 ✗ oppositeSupport = ielOppositeGlobalSup[jelSup];
2085
2086 // compute dof in common
2087 ✗ std::size_t countSupportDofInCommon = 0;
2088 ✗ for(std::size_t iSup1=0; iSup1<sizeCurrent; ++iSup1) {
2089 ✗ for(std::size_t iSup2=0; iSup2<sizeOpposite; ++iSup2)
2090 ✗ if( iSupportDof[iEle[oppositeSupport] + iSup2] == iSupportDof[iEle[currentSupport] + iSup1] )
2091 ✗ ++countSupportDofInCommon;
2092 }
2093
2094 // if the two support elements share an edge/face
2095 ✗ if( countSupportDofInCommon == sizeCurrent-1 )
2096 ✗ computeFaceOriented = true;
2097 }
2098 }
2099
2100 // if there really is an opposite support element
2101 // Even for an inner edge, it is possible that the current support element has no opposite
2102 // support element (if it is duplicated and the opposite element is not).
2103 ✗ if( computeFaceOriented ) {
2104 // update the opposite finite elements
2105 ✗ oppositeEltFEWithBdVel[eltTypeOpp]->updateFirstDeriv(0, oppositeElemPoint);
2106 ✗ oppositeEltFEWithBdVel[eltTypeOpp]->updateBdMeasNormal();
2107 ✗ oppositeEltFEWithBdPre[eltTypeOpp]->updateFirstDeriv(0, oppositeElemPoint);
2108 ✗ oppositeEltFEWithBdPre[eltTypeOpp]->updateBdMeasNormal();
2109
2110 // compute coefficient
2111 ✗ elemFieldAdvFace.initialize(DOF_FIELD, m_feVelWithBd->bdEle(iFace), m_velocity->numComponent());
2112 ✗ elemFieldAdvFaceTmp.initialize(DOF_FIELD, m_feVelWithBd->bdEle(iFace), m_velocity->numComponent());
2113 ✗ if( FelisceParam::instance().useProjectionForNitscheXFEM ) {
2114 ✗ for(felInt iextrap = 0; iextrap < FelisceParam::instance().orderBdfNS; ++iextrap) {
2115 ✗ elemFieldAdvFaceTmp.setValue(m_seqVelExtrapol[iextrap], *m_feVelWithBd, iFace, currentSupport, m_iVelocity, m_ao, dof());
2116 ✗ elemFieldAdvFace.val += elemFieldAdvFaceTmp.val;
2117 }
2118 } else {
2119 ✗ for(felInt iextrap = 0; iextrap < FelisceParam::instance().orderBdfNS; ++iextrap) {
2120 ✗ elemFieldAdvFaceTmp.setValue(m_seqVelExtrapol[iextrap], *m_feVelWithBd, iFace, ielOld[iextrap], m_iVelocity, m_aoOld[iextrap], m_dofOld[iextrap]);
2121 ✗ elemFieldAdvFace.val += elemFieldAdvFaceTmp.val;
2122 }
2123 }
2124
2125 ✗ UBlasVector normalVel = prod(trans(elemFieldAdvFace.val), m_feVelWithBd->bdEle(iFace).normal[0]);
2126 ✗ double normVelBd = norm_inf(normalVel);
2127 ✗ double gammaU = coeffVel * normVelBd;
2128
2129 // We can now compute all the integrals.
2130 // There is a minus sign for gammaP because there is one for q div(u).
2131 // current element for phi_i and phi_j
2132 ✗ m_elementMat[0]->zero();
2133 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n( gammaU, *m_feVelWithBd, *m_feVelWithBd, iFace, iFace, m_velBlock, m_velBlock, m_velocity->numComponent());
2134 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(-gammaP, *m_fePreWithBd, *m_fePreWithBd, iFace, iFace, m_preBlock, m_preBlock, m_pressure->numComponent());
2135 ✗ setValueMatrixRHS(currentSupport, currentSupport, FlagMatrixRHS::only_matrix);
2136
2137 // current element for phi_i and opposite element for phi_j
2138 ✗ m_elementMat[0]->zero();
2139 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n( gammaU, *oppositeEltFEWithBdVel[eltTypeOpp], *m_feVelWithBd, localIdFaceOpposite, iFace, m_velBlock, m_velBlock, m_velocity->numComponent());
2140 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(-gammaP, *oppositeEltFEWithBdPre[eltTypeOpp], *m_fePreWithBd, localIdFaceOpposite, iFace, m_preBlock, m_preBlock, m_pressure->numComponent());
2141 ✗ setValueMatrixRHS(currentSupport, oppositeSupport, FlagMatrixRHS::only_matrix);
2142
2143 // current element for phi_i and opposite element for phi_j
2144 ✗ m_elementMat[0]->mat() = trans(m_elementMat[0]->mat());
2145 ✗ setValueMatrixRHS(oppositeSupport, currentSupport, FlagMatrixRHS::only_matrix);
2146
2147 // current element for phi_i and phi_j
2148 ✗ m_elementMat[0]->zero();
2149 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n( gammaU, *oppositeEltFEWithBdVel[eltTypeOpp], *oppositeEltFEWithBdVel[eltTypeOpp], localIdFaceOpposite, localIdFaceOpposite, m_velBlock, m_velBlock, m_velocity->numComponent());
2150 ✗ m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(-gammaP, *oppositeEltFEWithBdPre[eltTypeOpp], *oppositeEltFEWithBdPre[eltTypeOpp], localIdFaceOpposite, localIdFaceOpposite, m_preBlock, m_preBlock, m_pressure->numComponent());
2151 ✗ setValueMatrixRHS(oppositeSupport, oppositeSupport, FlagMatrixRHS::only_matrix);
2152 }
2153 }
2154 }
2155 }
2156 ✗ ++countElement[eltType];
2157 ✗ ++ielCurrentLocalGeo;
2158 }
2159 }
2160 }
2161
2162 // desallocate array use for assemble with petsc.
2163 ✗ desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_matrix);
2164
2165 // desallocate opposite finite elements
2166 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
2167 ✗ eltType = bagElementTypeDomain[i];
2168 ✗ delete oppositeEltFEWithBdVel[eltType];
2169 ✗ delete oppositeEltFEWithBdPre[eltType];
2170 }
2171 }
2172
2173 ✗ void LinearProblemNitscheXFEM::computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*ielSupportDof*/, FlagMatrixRHS /*flagMatrixRHS*/)
2174 {
2175
2176 // call to the function in linear problem that updates the curvilinear finite element
2177 ✗ if( m_curvFeVel->hasOriginalQuadPoint() ) {
2178 ✗ m_curvFeVel->updateMeasNormalContra(0, elemPoint);
2179 ✗ m_curvFePress->updateMeasNormalContra(0, elemPoint);
2180 } else {
2181 ✗ m_curvFeVel->updateBasisAndNormalContra(0, elemPoint);
2182 ✗ m_curvFePress->updateBasisAndNormalContra(0, elemPoint);
2183 }
2184 }
2185
2186 ✗ void LinearProblemNitscheXFEM::computeAndApplyElementNaturalBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& ielSupportDof1, felInt& ielSupportDof2, felInt& ielGeoGlobal, FlagMatrixRHS flagMatrixRHS)
2187 {
2188
2189 ✗ if( ielSupportDof1 == ielSupportDof2 ) {
2190
2191 // get the side of the sub elements
2192 ✗ std::vector<sideOfInterface> bndElemSide;
2193 ✗ m_duplicateSupportElements->getSubBndEltSide(ielGeoGlobal, bndElemSide);
2194
2195 // check if the face belongs to the intersected mesh
2196 ✗ if( bndElemSide.size() > 0 ) {
2197
2198 // update finite elements
2199 ✗ computeElementArrayBoundaryCondition(elemPoint, elemIdPoint, ielSupportDof1, flagMatrixRHS);
2200
2201 // detect current side
2202 felInt ielSupport;
2203 ✗ m_supportDofUnknown[0].getIdElementSupport(ielGeoGlobal, ielSupport);
2204 ✗ const sideOfInterface side = ielSupport == ielSupportDof2 ? LEFT : RIGHT;
2205
2206 // get sub-faces coordinates
2207 ✗ std::vector< std::vector< Point > > bndElemCrd;
2208 ✗ m_duplicateSupportElements->getSubBndElt(ielGeoGlobal, bndElemCrd);
2209
2210 // loop over the sub-faces
2211 ✗ for(std::size_t iSubElt = 0; iSubElt < bndElemSide.size(); ++iSubElt) {
2212
2213 // check if sub-face is the good side
2214 ✗ if( bndElemSide[iSubElt] != side )
2215 ✗ continue;
2216
2217 // get sub-face coordinates
2218 ✗ for(std::size_t iPt = 0; iPt < m_numVerPerSolidElt; ++iPt)
2219 ✗ *m_subItfPts[iPt] = bndElemCrd[iSubElt][iPt];
2220
2221 // update finite elements
2222 ✗ m_curvFeVel->updateSubElementMeasNormal(0, elemPoint, m_subItfPts);
2223 ✗ m_curvFePress->updateSubElementMeasNormal(0, elemPoint, m_subItfPts);
2224
2225 // compute specific term of users
2226 ✗ userElementComputeNaturalBoundaryCondition(elemPoint, elemIdPoint, ielSupportDof1, ielSupportDof2, ielGeoGlobal, m_currentLabel);
2227
2228 // apply natural boundary condition
2229 ✗ applyNaturalBoundaryCondition(elemPoint, elemIdPoint, ielSupportDof1, flagMatrixRHS);
2230 }
2231 ✗ } else {
2232 ✗ computeElementArrayBoundaryCondition(elemPoint, elemIdPoint, ielSupportDof1, flagMatrixRHS);
2233
2234 ✗ userElementComputeNaturalBoundaryCondition(elemPoint, elemIdPoint, ielSupportDof1, ielSupportDof2, ielGeoGlobal, m_currentLabel);
2235
2236 ✗ applyNaturalBoundaryCondition(elemPoint, elemIdPoint, ielSupportDof1, flagMatrixRHS);
2237 }
2238 }
2239 }
2240
2241 ✗ void LinearProblemNitscheXFEM::assembleMatrixBCDarcy()
2242 {
2243
2244 ✗ const std::vector<int>& r_fusionlabel = FelisceParam::instance().FusionLabel;
2245 int currentLabel;
2246 ✗ FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::only_matrix;
2247 ✗ felInt numEltPerLabel, numElement = 0;
2248 ✗ felInt ielGeoGlobal, ielGeo = 0;
2249
2250 // contains ids of all the support elements associated to a mesh element
2251 ✗ std::vector<felInt> vectorIdSupport;
2252
2253 // use to get element number in support dof mesh ordering.
2254 felInt ielSupportDof;
2255
2256 // first loop on geometric type. (only one...)
2257 ✗ ElementType eltType = meshLocal()->bagElementTypeDomainBoundary()[0];
2258
2259 // number points per element
2260 ✗ const int numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
2261
2262 // get number domain elements
2263 ✗ const felInt numDomElt = meshLocal()->getNumDomainElement();
2264
2265 // contains points of the current element.
2266 ✗ std::vector<Point*> elemPoint(numPointPerElt, 0);
2267
2268 // contains ids of point of the current element.
2269 ✗ std::vector<felInt> elemIdPoint(numPointPerElt, 0);
2270
2271 // define all current finite element use in the problem from data
2272 // file configuration and allocate elemMat and elemVec
2273 ✗ defineCurvilinearFiniteElement(eltType);
2274
2275 // Element matrix and vector initialisation
2276 ✗ initElementArrayBD();
2277
2278 // allocate array use for assemble with petsc.
2279 ✗ allocateArrayForAssembleMatrixRHSBD(flagMatrixRHS);
2280
2281 // initialize the elementField for pressureDarcy
2282 ✗ initPerElementTypeDarcy(eltType, flagMatrixRHS);
2283
2284 // second loop on region of the mesh.
2285 ✗ for (auto itRef = meshLocal()->intRefToBegEndMaps[eltType].begin(); itRef != meshLocal()->intRefToBegEndMaps[eltType].end(); itRef++) {
2286 ✗ currentLabel = itRef->first;
2287 ✗ numEltPerLabel = itRef->second.second;
2288
2289 ✗ auto it_fus = std::find(r_fusionlabel.begin(), r_fusionlabel.end(), currentLabel);
2290 ✗ if ( it_fus != r_fusionlabel.end() ) {
2291
2292 // third loop on element in the region with the type: eltType. ("real" loop on elements).
2293 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
2294
2295 ✗ m_supportDofUnknownLocal[m_iPreDarcy].getIdElementSupport(eltType, numElement, vectorIdSupport);
2296 ✗ ISLocalToGlobalMappingApply(*m_mappingElemSupportPerUnknown[m_iPreDarcy], vectorIdSupport.size(), vectorIdSupport.data(), vectorIdSupport.data());
2297
2298 ✗ meshLocal()->getOneElement(eltType, numElement, elemIdPoint, elemPoint, 0);
2299
2300 ✗ ielGeo = numElement + numDomElt;
2301 ✗ ISLocalToGlobalMappingApply(mappingElem(), 1, &ielGeo, &ielGeoGlobal);
2302
2303 // fourth loop over all the support elements
2304 ✗ for (std::size_t it = 0; it < vectorIdSupport.size(); it++) {
2305
2306 // get the id of the support
2307 ✗ ielSupportDof = vectorIdSupport[it];
2308
2309 // clear elementary matrices
2310 ✗ for (std::size_t j=0; j<m_elementMatBD.size(); j++)
2311 ✗ m_elementMatBD[j]->zero();
2312
2313 // compute Darcy term on element
2314 ✗ computeAndApplyElementDarcy(elemPoint, ielSupportDof, ielSupportDof, ielGeoGlobal, flagMatrixRHS);
2315
2316 // add values of elemMat in the global matrix: m_matrices[0].
2317 ✗ setValueMatrixRHSBD(ielSupportDof, ielSupportDof, flagMatrixRHS);
2318 }
2319
2320 ✗ numElement++;
2321 }
2322 } else {
2323 ✗ numElement += numEltPerLabel;
2324 }
2325 }
2326
2327 //desallocate array use for assemble with petsc.
2328 ✗ desallocateArrayForAssembleMatrixRHS(flagMatrixRHS);
2329 }
2330
2331 ✗ void LinearProblemNitscheXFEM::computeAndApplyElementDarcy(const std::vector<Point*>& elemPoint, felInt ielSupportDof1, felInt ielSupportDof2, felInt ielGeoGlobal, FlagMatrixRHS flagMatrixRHS)
2332 {
2333 IGNORE_UNUSED_ARGUMENT(flagMatrixRHS);
2334 ✗ if(ielSupportDof1 == ielSupportDof2) {
2335
2336 // get the side of the sub elements
2337 ✗ std::vector<sideOfInterface> bndElemSide;
2338 ✗ m_duplicateSupportElements->getSubBndEltSide(ielGeoGlobal, bndElemSide);
2339
2340 // check if the face belongs to the intersected mesh
2341 ✗ if( bndElemSide.size() > 0 ) {
2342
2343 // update finite elements
2344 ✗ computeElementArrayDarcy(elemPoint);
2345
2346 // detect current side
2347 felInt ielSupport;
2348 ✗ m_supportDofUnknown[0].getIdElementSupport(ielGeoGlobal, ielSupport);
2349 ✗ const sideOfInterface side = ielSupport == ielSupportDof2 ? LEFT : RIGHT;
2350
2351 // get sub-faces coordinates
2352 ✗ std::vector< std::vector< Point > > bndElemCrd;
2353 ✗ m_duplicateSupportElements->getSubBndElt(ielGeoGlobal, bndElemCrd);
2354
2355 // loop over the sub-faces
2356 ✗ for(std::size_t iSubElt=0; iSubElt<bndElemSide.size(); ++iSubElt) {
2357
2358 // check if sub-face is the good side
2359 ✗ if(bndElemSide[iSubElt] != side)
2360 ✗ continue;
2361
2362 // get sub-face coordinates
2363 ✗ for(std::size_t iPt=0; iPt < m_numVerPerSolidElt; ++iPt)
2364 ✗ *m_subItfPts[iPt] = bndElemCrd[iSubElt][iPt];
2365
2366 // update finite elements
2367 ✗ m_curvFeVel->updateSubElementMeasNormal(0, elemPoint, m_subItfPts);
2368 ✗ m_curvFeDarcy->updateSubElementMeasNormal(0, elemPoint, m_subItfPts);
2369
2370 // compute Darcy term
2371 ✗ elementComputeDarcyOnBoundary(elemPoint);
2372 }
2373 ✗ } else {
2374 // update finite elements
2375 ✗ computeElementArrayDarcy(elemPoint);
2376
2377 // compute Darcy term
2378 ✗ elementComputeDarcyOnBoundary(elemPoint);
2379 }
2380 }
2381 }
2382
2383 ✗ void LinearProblemNitscheXFEM::setStrucFiniteElement(CurvilinearFiniteElement* strucFE)
2384 {
2385
2386 ✗ m_curvFeStruc = strucFE;
2387 }
2388
2389 ✗ void LinearProblemNitscheXFEM::addMatrixRHS()
2390 {
2391
2392 ✗ this->vector().axpy(1.0,vector(1));
2393 }
2394
2395 ✗ void LinearProblemNitscheXFEM::setDuplicateSupportObject(DuplicateSupportDof* object)
2396 {
2397
2398 ✗ m_duplicateSupportElements = object;
2399 }
2400
2401 ✗ void LinearProblemNitscheXFEM::gatherSeqVelExtrapol(std::vector<PetscVector>& parallelVec)
2402 {
2403 ✗ for(std::size_t i=0; i<parallelVec.size(); ++i)
2404 ✗ gatherVector(parallelVec[i], m_seqVelExtrapol[i]);
2405 ✗ m_isSeqVelExtrapolAllocated = true;
2406 }
2407
2408 ✗ void LinearProblemNitscheXFEM::gatherSeqBdfRHS(std::vector<PetscVector>& parallelVec)
2409 {
2410 ✗ for(std::size_t i=0; i<parallelVec.size(); ++i)
2411 ✗ gatherVector(parallelVec[i], m_seqBdfRHS[i]);
2412 ✗ m_isSeqBdfRHSAllocated = true;
2413 }
2414
2415 ✗ void LinearProblemNitscheXFEM::deleteDynamicData()
2416 {
2417 // here, we delete all the matrices, vectors and mapping of the linear problem.
2418 // everything will be rebuilt at the next iteration
2419
2420 // sequential solution
2421 ✗ bool& isSeqSolAlloc = isSequentialSolutionAllocated();
2422 ✗ if(isSeqSolAlloc) {
2423 ✗ sequentialSolution().destroy();
2424 ✗ isSeqSolAlloc = false;
2425 }
2426
2427 // element matrices and matrices
2428 ✗ for (unsigned int i = 0; i < numberOfMatrices(); i++)
2429 ✗ matrix(i).destroy();
2430
2431 ✗ m_elementMat.clear();
2432 ✗ m_elementVector.clear();
2433 ✗ m_elementMatBD.clear();
2434 ✗ m_elementVectorBD.clear();
2435
2436 // solution and rhs
2437 ✗ bool& areSolAndRhsAlloc = areSolutionAndRHSAllocated();
2438 ✗ if (areSolAndRhsAlloc) {
2439 ✗ for (unsigned int i = 0; i < numberOfVectors(); i++)
2440 ✗ vector(i).destroy();
2441
2442 ✗ solution().destroy();
2443 ✗ areSolAndRhsAlloc = false;
2444 }
2445
2446 // finite elements
2447 ✗ for (std::size_t i = 0; i < m_listCurrentFiniteElement.size(); i++)
2448 ✗ delete m_listCurrentFiniteElement[i];
2449 ✗ m_listCurrentFiniteElement.listCurrentFiniteElement().clear();
2450
2451 ✗ for (std::size_t i = 0; i < m_listCurvilinearFiniteElement.size(); i++)
2452 ✗ delete m_listCurvilinearFiniteElement[i];
2453 ✗ m_listCurvilinearFiniteElement.listCurvilinearFiniteElement().clear();
2454
2455 // mappings
2456 ✗ bool& isMapElemAlloc = isMappingElemAllocated();
2457 ✗ if(isMapElemAlloc) {
2458 ✗ ISLocalToGlobalMappingDestroy(&mappingElem());
2459 ✗ isMapElemAlloc = false;
2460 }
2461
2462 ✗ bool& isMapElemSupportAlloc = isMappingElemSupportAllocated();
2463 ✗ if(isMapElemSupportAlloc) {
2464 ✗ for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++ )
2465 ✗ ISLocalToGlobalMappingDestroy(mappingElemSupportPerUnknown(iUnknown));
2466
2467 ✗ isMapElemSupportAlloc = false;
2468 }
2469
2470 ✗ bool& isMapNodesAlloc = isMappingNodesAllocated();
2471 ✗ if(isMapNodesAlloc) {
2472 ✗ ISLocalToGlobalMappingDestroy(&mappingNodes());
2473 ✗ isMapNodesAlloc = false;
2474 }
2475
2476 ✗ bool& isMapFelToPetscAlloc = isMappingLocalFelisceToGlobalPetscAllocated();
2477 ✗ if(isMapFelToPetscAlloc) {
2478 ✗ for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++ )
2479 ✗ ISLocalToGlobalMappingDestroy(&mappingIdFelisceToPetsc(iUnknown));
2480
2481 ✗ isMapFelToPetscAlloc = false;
2482 }
2483
2484 // AO (copy the AO, it will be destroy with the old one)
2485 ✗ bool& isAOAlloc = isAOAllocated();
2486 ✗ if(isAOAlloc) {
2487 ✗ AODestroy(&ao());
2488 ✗ isAOAlloc = false;
2489 }
2490
2491 // dofs (the pattern is useless for the old dof)
2492 ✗ dof().numDof() = 0;
2493 ✗ dof().numDofPerUnknown().clear();
2494 ✗ dof().clearPattern();
2495
2496 // repartition of the dof
2497 ✗ m_dofPart.clear();
2498 ✗ m_eltPart.clear();
2499
2500 // supportDofMesh (copy the old one)
2501 ✗ m_supportDofUnknown.clear();
2502 ✗ m_supportDofUnknownLocal.clear();
2503
2504 // local mesh
2505 ✗ meshLocal()->domainDim() = GeometricMeshRegion::GeoMeshUndefined;
2506 ✗ meshLocal()->deleteElements();
2507
2508 // extrapolation and time rhs of this class
2509 ✗ if(m_isSeqVelExtrapolAllocated) {
2510 ✗ for(std::size_t i=0; i<m_seqVelExtrapol.size(); ++i) {
2511 ✗ m_seqVelExtrapol[i].destroy();
2512 ✗ m_seqVelExtrapol[i] = PetscVector::null();
2513 }
2514 ✗ m_isSeqVelExtrapolAllocated = false;
2515 }
2516
2517 ✗ if(m_isSeqBdfRHSAllocated) {
2518 ✗ for(std::size_t i=0; i<m_seqBdfRHS.size(); ++i) {
2519 ✗ m_seqBdfRHS[i].destroy();
2520 ✗ m_seqBdfRHS[i] = PetscVector::null();
2521 }
2522 ✗ m_isSeqBdfRHSAllocated = false;
2523 }
2524
2525 // KSP
2526 ✗ m_KSPInterface->destroy();
2527 }
2528
2529 ✗ void LinearProblemNitscheXFEM::initOldObjects()
2530 {
2531 // compute how many old solution we need to store
2532 ✗ felInt order = FelisceParam::instance().orderBdfNS;
2533
2534 // resize vector
2535 ✗ m_sequentialSolutionOld.resize(order);
2536 ✗ m_solutionOld.resize(order);
2537 ✗ m_aoOld.resize(order);
2538 ✗ m_supportDofUnknownOld.resize(order);
2539 ✗ m_dofOld.resize(order);
2540
2541 // init old solution
2542 ✗ for(std::size_t i=0; i<m_sequentialSolutionOld.size(); ++i) {
2543 ✗ m_sequentialSolutionOld[i].duplicateFrom(sequentialSolution());
2544 ✗ m_sequentialSolutionOld[i].copyFrom(sequentialSolution());
2545 }
2546
2547 ✗ for(std::size_t i=0; i<m_solutionOld.size(); ++i) {
2548 ✗ m_solutionOld[i].duplicateFrom(solution());
2549 ✗ m_solutionOld[i].copyFrom(solution());
2550 }
2551
2552 // init old AO
2553 ✗ for(std::size_t i=0; i<m_aoOld.size(); ++i) {
2554 ✗ m_aoOld[i] = ao();
2555 ✗ PetscObjectReference((PetscObject)ao());
2556 }
2557
2558 // init old support dof
2559 ✗ for(std::size_t i=0; i<m_supportDofUnknownOld.size(); ++i)
2560 ✗ m_supportDofUnknownOld[i] = m_supportDofUnknown;
2561
2562 // init old dof
2563 ✗ std::vector<SupportDofMesh*> pSupportDofUnknown(m_supportDofUnknown.size(), NULL);
2564 ✗ for(std::size_t i=0; i<m_dofOld.size(); ++i) {
2565 ✗ for(std::size_t iUnknown=0; iUnknown < m_supportDofUnknown.size(); ++iUnknown)
2566 ✗ pSupportDofUnknown[iUnknown] = &m_supportDofUnknownOld[i][iUnknown];
2567
2568 ✗ m_dofOld[i].setDof(m_listUnknown, m_listVariable, pSupportDofUnknown);
2569 ✗ m_dofOld[i].numDof() = dof().numDof();
2570 ✗ m_dofOld[i].numDofPerUnknown() = dof().numDofPerUnknown();
2571 }
2572 }
2573
2574 ✗ void LinearProblemNitscheXFEM::updateOldSolution(felInt countInitOld)
2575 {
2576
2577 ✗ if(countInitOld > 0) {
2578
2579 // sequential solution
2580 ✗ m_sequentialSolutionOld[m_sequentialSolutionOld.size() - 1].destroy();
2581 ✗ for(std::size_t i=m_sequentialSolutionOld.size() - 1; i>0; --i)
2582 ✗ m_sequentialSolutionOld[i] = m_sequentialSolutionOld[i-1];
2583
2584 ✗ m_sequentialSolutionOld[0].duplicateFrom(sequentialSolution());
2585 ✗ m_sequentialSolutionOld[0].copyFrom(sequentialSolution());
2586
2587 // parallel solution
2588 ✗ m_solutionOld[m_solutionOld.size() - 1].destroy();
2589 ✗ for(std::size_t i=m_solutionOld.size() - 1; i>0; --i)
2590 ✗ m_solutionOld[i] = m_solutionOld[i-1];
2591
2592 ✗ m_solutionOld[0].duplicateFrom(solution());
2593 ✗ m_solutionOld[0].copyFrom(solution());
2594
2595 // AO
2596 ✗ AODestroy(&m_aoOld[m_aoOld.size() - 1]);
2597 ✗ for(std::size_t i=m_aoOld.size() - 1; i>0; --i)
2598 ✗ m_aoOld[i] = m_aoOld[i-1];
2599 ✗ m_aoOld[0] = ao();
2600 ✗ PetscObjectReference((PetscObject)ao());
2601
2602 // supportDofMesh (copy the old one)
2603 ✗ for(std::size_t i=m_supportDofUnknownOld.size() - 1; i>0; --i)
2604 ✗ m_supportDofUnknownOld[i] = m_supportDofUnknownOld[i-1];
2605 ✗ m_supportDofUnknownOld[0] = m_supportDofUnknown;
2606
2607 // dofs
2608 ✗ for(std::size_t i=m_dofOld.size() - 1; i>0; --i) {
2609 ✗ m_dofOld[i].numDof() = m_dofOld[i-1].numDof();
2610 ✗ m_dofOld[i].numDofPerUnknown() = m_dofOld[i-1].numDofPerUnknown();
2611 }
2612
2613 ✗ m_dofOld[0].numDof() = dof().numDof();
2614 ✗ m_dofOld[0].numDofPerUnknown() = dof().numDofPerUnknown();
2615
2616 // delete extrapolations
2617 ✗ for(std::size_t i=0; i<m_seqVelExtrapol.size(); ++i) {
2618 ✗ m_seqVelExtrapol[i].destroy();
2619 ✗ m_seqVelExtrapol[i] = PetscVector::null();
2620 }
2621 ✗ m_isSeqVelExtrapolAllocated = false;
2622
2623 // delete rhs time
2624 ✗ for(std::size_t i=0; i<m_seqBdfRHS.size(); ++i) {
2625 ✗ m_seqBdfRHS[i].destroy();
2626 ✗ m_seqBdfRHS[i] = PetscVector::null();
2627 }
2628 ✗ m_isSeqBdfRHSAllocated = false;
2629 } else {
2630
2631 // copy the solutions
2632 ✗ for(std::size_t i=m_sequentialSolutionOld.size() - 1; i>0; --i)
2633 ✗ m_sequentialSolutionOld[i].copyFrom(m_sequentialSolutionOld[i-1]);
2634 ✗ m_sequentialSolutionOld[0].copyFrom(sequentialSolution());
2635
2636 ✗ for(std::size_t i=m_solutionOld.size() - 1; i>0; --i)
2637 ✗ m_solutionOld[i].copyFrom(m_solutionOld[i-1]);
2638 ✗ m_solutionOld[0].copyFrom(solution());
2639 }
2640 }
2641
2642 ✗ void LinearProblemNitscheXFEM::setInterfaceExchangedData(std::vector<double>& intfVelo, std::vector<double>& intfForc)
2643 {
2644 ✗ m_intfVelo = &intfVelo;
2645 ✗ m_intfForc = &intfForc;
2646 }
2647
2648 ✗ void LinearProblemNitscheXFEM::computeInterfaceStress()
2649 {
2650
2651 ✗ FEL_CHECK(m_intfForc != nullptr, "Error: the pointer of the interface force has not been initialized yet.");
2652
2653 ✗ std::fill(m_intfForc->begin(), m_intfForc->end(), 0);
2654
2655 felInt numEltPerLabel;
2656 int currentLabel;
2657
2658 // Get fictitious interface labels
2659 ✗ const std::vector<int>& fictitious = FelisceParam::instance().fictitious;
2660
2661 // Get element type
2662 ✗ const ElementType eltType = m_interfMesh->bagElementTypeDomain()[0];
2663 ✗ const ElementType eltTypeFl = meshLocal()->bagElementTypeDomain()[0];
2664
2665 // Resize array
2666 ✗ const int numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
2667 ✗ std::vector<Point*> elemPoint(numPointPerElt, 0);
2668 ✗ std::vector<felInt> elemIdPoint(numPointPerElt, 0);
2669
2670 // Define finite element
2671 ✗ defineFiniteElement(eltTypeFl);
2672
2673 // Element matrix and vector initialisation
2674 ✗ initElementArray();
2675
2676 // Get the finite elements
2677 ✗ m_feVel = m_listCurrentFiniteElement[m_iVelocity];
2678 ✗ m_fePres = m_listCurrentFiniteElement[m_iPressure];
2679
2680 // Initialize element fields
2681 ✗ m_elemFieldSolidRHS.initialize(QUAD_POINT_FIELD, *m_curvFeStruc, m_curvFeStruc->numCoor());
2682 ✗ m_seqVelDofPts.initialize(DOF_FIELD, *m_feVel, m_velocity->numComponent());
2683 ✗ m_seqPreDofPts.initialize(DOF_FIELD, *m_fePres, m_pressure->numComponent());
2684
2685 ✗ ElementVector elemVec( {m_curvFeStruc}, {static_cast<std::size_t>(m_curvFeStruc->numCoor())} );
2686
2687 // Used to define a "global" numbering of element in the mesh.
2688 ✗ felInt numElement = 0;
2689
2690 // Loop on element in the region with the type: eltType
2691 ✗ for(auto itRef = m_interfMesh->intRefToBegEndMaps[eltType].begin(); itRef != m_interfMesh->intRefToBegEndMaps[eltType].end(); itRef++) {
2692 ✗ currentLabel = itRef->first;
2693 ✗ numEltPerLabel = itRef->second.second;
2694
2695 // if fictitious interface skip
2696 ✗ auto it_fic = find (fictitious.begin(), fictitious.end(), currentLabel);
2697 ✗ if ( it_fic == fictitious.end() ) {
2698
2699 ✗ for(felInt iel = 0; iel < numEltPerLabel; iel++) {
2700
2701 // return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
2702 ✗ m_interfMesh->getOneElement(eltType, numElement, elemIdPoint, elemPoint, 0);
2703
2704 // clear elementary vector.
2705 ✗ elemVec.zero();
2706
2707 // compute stress on interface element
2708 ✗ m_computeElementaryInterfaceStress(elemPoint, elemIdPoint, eltType, numElement, elemVec);
2709
2710 // assembly m_intForc
2711 ✗ for(int ic = 0; ic < m_curvFeStruc->numCoor(); ++ic) {
2712 ✗ for(int ip = 0; ip < m_curvFeStruc->numPoint(); ++ip) { // TODO should be numDof but if P1 --> numPoint ok
2713 ✗ int gdof = elemIdPoint[ip];
2714 ✗ (*m_intfForc)[m_curvFeStruc->numCoor()*gdof+ic] += elemVec.vec()[ip+ic*m_curvFeStruc->numDof()];
2715 }
2716 }
2717
2718 ✗ numElement++;
2719 }
2720 } else {
2721
2722 ✗ numElement += numEltPerLabel;
2723 }
2724 }
2725 }
2726
2727 ✗ void LinearProblemNitscheXFEM::printStructStress(int indexTime, double dt)
2728 {
2729
2730 ✗ if ( MpiInfo::rankProc() != 0)
2731 ✗ return;
2732
2733 ✗ int nbp = m_interfMesh->numPoints();
2734
2735 // write case
2736 FILE * pFile;
2737 ✗ std::string fileName, fileNameDir, fileNameCase;
2738
2739 ✗ fileNameDir = FelisceParam::instance().resultDirStresses;
2740 ✗ fileNameCase = fileNameDir+"stresses.case";
2741
2742 ✗ pFile = fopen (fileNameCase.c_str(),"w");
2743 ✗ if (pFile==NULL) {
2744 ✗ std::string command = "mkdir -p " + fileNameDir;
2745 ✗ int ierr = system(command.c_str());
2746 ✗ if( ierr > 0)
2747 ✗ FEL_ERROR("Impossible to write "+fileNameCase+" case file");
2748 ✗ pFile = fopen (fileNameCase.c_str(),"w");
2749 }
2750 ✗ fprintf( pFile,"FORMAT\n");
2751 ✗ fprintf( pFile,"type: ensight\n");
2752 ✗ fprintf( pFile,"GEOMETRY\n");
2753 ✗ fprintf( pFile,"model: solid.geo\n");
2754 ✗ fprintf( pFile,"VARIABLE\n");
2755 ✗ fprintf( pFile,"vector per node: 1 stresses stresses.*****.vct\n");
2756 ✗ fprintf( pFile,"TIME\n");
2757 ✗ fprintf( pFile, "time set: %d\n", 1);
2758 ✗ fprintf( pFile,"number of steps: %d \n", indexTime+1);
2759 ✗ fprintf( pFile, "filename start number: %d\n", 0);
2760 ✗ fprintf( pFile, "filename increment: %d\n", 1);
2761 ✗ fprintf( pFile, "time values:\n");
2762 ✗ int count=0;
2763
2764 ✗ for (int i=0; i < indexTime+1; ++i) {
2765 ✗ fprintf( pFile, "%12.5e",i*dt);
2766 ✗ ++count;
2767 ✗ if ( count == 6) {
2768 ✗ fprintf( pFile, "\n" );
2769 ✗ count=0;
2770 }
2771 }
2772 ✗ fclose(pFile);
2773
2774 // file .vct
2775 ✗ fileName = fileNameDir+"stresses";
2776 ✗ std::stringstream oss;
2777 ✗ oss << indexTime;
2778 ✗ std::string aux = oss.str();
2779 ✗ if ( indexTime < 10 ) {
2780 ✗ aux = "0000" + aux;
2781 }
2782 ✗ else if ( indexTime < 100 ) {
2783 ✗ aux = "000" + aux;
2784 }
2785 ✗ else if ( indexTime < 1000 ) {
2786 ✗ aux = "00" + aux;
2787 }
2788 ✗ else if ( indexTime < 10000 ) {
2789 ✗ aux = "0" + aux;
2790 }
2791
2792 ✗ fileName = fileName + "." + aux + ".vct";
2793
2794 ✗ pFile = fopen (fileName.c_str(),"w");
2795
2796 ✗ if (pFile==NULL) {
2797 ✗ FEL_ERROR("Impossible to write " + fileName + " vector ensight solution");
2798 }
2799
2800 ✗ count=0;
2801 ✗ int nDimensions = dimension();
2802 ✗ fprintf( pFile, "Vector per node\n");
2803 // [u_x_1, u_y_1, u_z_1,u_x_2, u_y_2, u_z_2, ..., u_x_dim, u_y_dim, u_z_dim ]
2804 ✗ for ( int i=0; i < nbp; ++i ) {
2805 ✗ for ( int j=0; j < nDimensions; ++j ) {
2806 ✗ if(indexTime==0){
2807 ✗ fprintf(pFile,"%12.5e", 0.);
2808 ✗ ++count;
2809 } else {
2810 ✗ fprintf(pFile,"%12.5e", (*m_intfForc)[ i*nDimensions+j ] );
2811 ✗ ++count;
2812 }
2813
2814 ✗ if ( count == 6 ) {
2815 ✗ fprintf(pFile, "\n" );
2816 ✗ count=0;
2817 }
2818 }
2819 }
2820 ✗ fclose( pFile );
2821 }
2822
2823 ✗ void LinearProblemNitscheXFEM::computeElementMatrixRHS()
2824 {
2825 // mass coefficient
2826 ✗ const double coef = m_density/m_fstransient->timeStep;
2827
2828 //================================
2829 // matrix
2830 //================================
2831 // Laplacian operator
2832 ✗ if (m_useSymmetricStress) {
2833 // \mu \int \epsilon u^{n+1} : \epsilon v
2834 ✗ m_elementMat[0]->eps_phi_i_eps_phi_j(2*m_viscosity,*m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
2835 } else {
2836 // \mu \int \nabla u^{n+1} : \nabla v
2837 ✗ m_elementMat[0]->grad_phi_i_grad_phi_j(m_viscosity, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
2838 }
2839
2840 // div and grad operator
2841 // \int q \nabla \cdot u^{n+1} and \int p \nabla \cdot v
2842 ✗ m_elementMat[0]->psi_j_div_phi_i_and_psi_i_div_phi_j(*m_feVel, *m_fePres, m_velBlock, m_preBlock, -1, -1);
2843
2844 // mass matrix from the time scheme
2845 // \int u^{n+1} \cdot v
2846 ✗ m_elementMat[0]->phi_i_phi_j(m_bdf->coeffDeriv0()*coef, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
2847 }
2848
2849 ✗ void LinearProblemNitscheXFEM::computeElementMatrixRHSPartDependingOnOldTimeStep(felInt iel, felInt ielOld, felInt idTimeStep)
2850 {
2851 // --------------- //
2852 // convective term //
2853 // --------------- //
2854 ✗ switch( m_useNSEquation ) {
2855 ✗ case 0:
2856 // Stoke's model, no convective term
2857 ✗ break;
2858
2859 ✗ case 1:
2860 // classic Navier Stokes convective term : \rho (u \cdot \grad) u
2861 // Natural condition: \sigma(u,p) \cdot n = 0
2862 // discretization : \rho (u^{n} \cdot \grad) u^{n+1}
2863 // DEFAULT METHOD
2864 ✗ if(FelisceParam::instance().useProjectionForNitscheXFEM)
2865 ✗ m_elemFieldAdv.setValue(m_seqVelExtrapol[idTimeStep], *m_feVel, iel, m_iVelocity, m_ao, dof());
2866 else
2867 ✗ m_elemFieldAdv.setValue(m_seqVelExtrapol[idTimeStep], *m_feVel, ielOld, m_iVelocity, m_aoOld[idTimeStep], m_dofOld[idTimeStep]);
2868
2869 ✗ m_elementMat[0]->u_grad_phi_j_phi_i(m_density, m_elemFieldAdv, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
2870
2871 // Temam's operator to stabilize the convective term
2872 // Convective term in energy balance :
2873 // ... + \rho/2 \int_Gamma (u^n.n)|u^n+1|^2 - \rho/2 \int_Omega (\nabla . u^n) |u^n+1|^2 + ...
2874 // Add \rho/2 \int_Omega (\nabla . u^n) u^n+1 v to delete the second term.
2875 // it's consistent with the exact solution (divergence = 0)
2876 ✗ m_elementMat[0]->div_u_phi_j_phi_i(0.5*m_density, m_elemFieldAdv, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
2877 ✗ break;
2878
2879 ✗ default:
2880 ✗ FEL_ERROR("Error: Wrong or not implemented convection method, check the param NSequationFlag in the data file");
2881 }
2882
2883 //================================
2884 // RHS
2885 //================================
2886 // bdf
2887 // \frac{\rho}{\dt} \sum_{i=1}^k \alpha_i u_{n+1-i}
2888 // 'm_density' and not 'coef' because the time step is integrated into the bdf term
2889 ✗ if(FelisceParam::instance().useProjectionForNitscheXFEM)
2890 ✗ m_elemFieldRHSbdf.setValue(m_seqBdfRHS[idTimeStep], *m_feVel, iel, m_iVelocity, m_ao, dof());
2891 else
2892 ✗ m_elemFieldRHSbdf.setValue(m_seqBdfRHS[idTimeStep], *m_feVel, ielOld, m_iVelocity, m_aoOld[idTimeStep], m_dofOld[idTimeStep]);
2893
2894 ✗ m_elementVector[1]->source(m_density, *m_feVel, m_elemFieldRHSbdf, m_velBlock, m_velocity->numComponent());
2895 }
2896
2897 ✗ void LinearProblemNitscheXFEM::computeNitscheSigma_u_v(double velCoeff, double preCoeff)
2898 {
2899 // \int_K (velCoeff \grad u vec + preCoeff p vec) \cdot v
2900 ✗ if(m_useSymmetricStress)
2901 ✗ m_elementMat[0]->phi_i_eps_phi_j_dot_vec(2.*velCoeff, m_elemFieldNormal, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
2902 else
2903 ✗ m_elementMat[0]->phi_i_grad_phi_j_dot_vec(velCoeff, m_elemFieldNormal, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
2904
2905 // same finite element for the pressure and the velocity
2906 ✗ for(std::size_t idim=0; idim<m_velocity->numComponent(); ++idim)
2907 ✗ m_elementMat[0]->phi_i_phi_j(preCoeff*m_elemFieldNormal.val(idim,0), *m_feVel, m_velBlock+idim, m_preBlock, m_pressure->numComponent());
2908 }
2909
2910 ✗ void LinearProblemNitscheXFEM::computeNitscheSigma_v_u(double velCoeff, double preCoeff)
2911 {
2912 // \int_K (velCoeff \grad v vec + preCoeff q vec) \cdot u
2913 ✗ if(m_useSymmetricStress)
2914 ✗ m_elementMat[0]->eps_phi_i_dot_vec_phi_j(2.*velCoeff, m_elemFieldNormal, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
2915 else
2916 ✗ m_elementMat[0]->grad_phi_i_dot_vec_phi_j(velCoeff, m_elemFieldNormal, *m_feVel, m_velBlock, m_velBlock, m_velocity->numComponent());
2917
2918 // same finite element for the pressure and velocity
2919 ✗ for(std::size_t idim=0; idim<m_velocity->numComponent(); ++idim)
2920 ✗ m_elementMat[0]->phi_i_phi_j(preCoeff*m_elemFieldNormal.val(idim,0), *m_feVel, m_preBlock, m_velBlock+idim, m_pressure->numComponent());
2921 }
2922
2923 ✗ void LinearProblemNitscheXFEM::computeNitscheDpoint_Sigma_v(double velCoeff, double preCoeff)
2924 {
2925 // \int_K ddot \cdot (velCoeff \grad v vec - preCoeff q vec)
2926 ✗ if(m_useSymmetricStress)
2927 ✗ m_elementVector[0]->eps_phi_i_dot_vec1_dot_vec2(2.*velCoeff, m_elemFieldNormal, m_strucVelQuadPts, *m_feVel, m_velBlock, m_velocity->numComponent());
2928 else
2929 ✗ m_elementVector[0]->grad_phi_i_dot_vec1_dot_vec2(velCoeff, m_elemFieldNormal, m_strucVelQuadPts, *m_feVel, m_velBlock, m_velocity->numComponent());
2930
2931 ✗ for(std::size_t idim=0; idim<m_velocity->numComponent(); ++idim) {
2932 ✗ for(felInt ig=0; ig<m_feVel->numQuadraturePoint(); ++ig) {
2933 ✗ m_ddotComp.val(0, ig) = m_strucVelQuadPts.val(idim, ig);
2934 }
2935 ✗ m_elementVector[0]->source(preCoeff*m_elemFieldNormal.val(idim,0), *m_feVel, m_ddotComp, m_preBlock, m_pressure->numComponent());
2936 }
2937 }
2938
2939 ✗ void LinearProblemNitscheXFEM::updateStructureVel(felInt strucId)
2940 {
2941 // update on the structure dof the local structure velocity from the master at the struct element strucId
2942 ✗ FEL_CHECK(m_intfVelo != nullptr, "Error: the pointer of the interface velocity has not been initialized yet.");
2943 // velocity comming from the FSI master ( The order is not the same! )
2944 // the local ordering is done by component while the global is by point
2945 ✗ std::vector<felInt> idDof(m_curvFeStruc->numDof());
2946 ✗ m_interfMesh->getOneElement(strucId,idDof);
2947 ✗ for(int ic=0; ic< m_curvFeStruc->numCoor(); ++ic) {
2948 ✗ for(int il=0; il< m_curvFeStruc->numDof(); ++il) {
2949 ✗ m_strucVelDofPts.val(ic, il) = (*m_intfVelo)[m_curvFeStruc->numCoor()*idDof[il] +ic];
2950 }
2951 }
2952 }
2953
2954 ✗ void LinearProblemNitscheXFEM::updateSubStructureVel(const std::vector<Point*>& ptElem, const std::vector<Point*>& ptSubElem)
2955 {
2956 // update on the SubStructure dof the local structure velocity from m_strucVelDofPts
2957 ✗ FEL_ASSERT(ptElem.size() == 2 || ptElem.size() == 3);
2958
2959 // compute the value of the displacement at the new integration points
2960 // get the coordinate of the integration point in the sub element
2961 ✗ m_curvFeStruc->update(0, ptSubElem);
2962 ✗ m_curvFeStruc->computeCurrentQuadraturePoint();
2963
2964 // compute the coordinate of these integration points in the reference element with the mapping of the
2965 // structure element (not the structure sub element)
2966 ✗ std::vector<Point> intPoints(m_curvFeStruc->numQuadraturePoint()); // Replace with m_computeIntegrationPoints
2967 ✗ if( ptElem.size() == 2){
2968 double dx, dy;
2969 ✗ dx = abs(ptElem[1]->x() - ptElem[0]->x());
2970 ✗ dy = abs(ptElem[1]->y() - ptElem[0]->y());
2971 ✗ if(dx >= dy) {
2972 ✗ for(int ig=0; ig<m_curvFeStruc->numQuadraturePoint(); ++ig) {
2973 ✗ intPoints[ig][0] = 2.*(m_curvFeStruc->currentQuadPoint[ig].x() - ptElem[0]->x())/(ptElem[1]->x() - ptElem[0]->x()) - 1.;
2974 ✗ intPoints[ig][1] = 0.;
2975 }
2976 } else {
2977 ✗ for(int ig=0; ig<m_curvFeStruc->numQuadraturePoint(); ++ig) {
2978 ✗ intPoints[ig][0] = 2.*(m_curvFeStruc->currentQuadPoint[ig].y() - ptElem[0]->y())/(ptElem[1]->y() - ptElem[0]->y()) - 1.;
2979 ✗ intPoints[ig][1] = 0;
2980 }
2981 }
2982 } else {
2983
2984 ✗ double det1 = (ptElem[2]->y() - ptElem[0]->y())*(ptElem[1]->x() - ptElem[0]->x())
2985 ✗ - (ptElem[2]->x() - ptElem[0]->x())*(ptElem[1]->y() - ptElem[0]->y());
2986
2987 ✗ double det2 = (ptElem[2]->x() - ptElem[0]->x())*(ptElem[1]->z() - ptElem[0]->z())
2988 ✗ - (ptElem[2]->z() - ptElem[0]->z())*(ptElem[1]->x() - ptElem[0]->x());
2989
2990 ✗ double det3 = (ptElem[2]->z() - ptElem[0]->z())*(ptElem[1]->y() - ptElem[0]->y())
2991 ✗ - (ptElem[2]->y() - ptElem[0]->y())*(ptElem[1]->z() - ptElem[0]->z());
2992
2993 ✗ if(abs(det1) >= abs(det2) && abs(det1) >= abs(det3)){
2994 ✗ for(int ig=0; ig<m_curvFeStruc->numQuadraturePoint(); ++ig) {
2995 ✗ intPoints[ig][0] = (ptElem[2]->y() - ptElem[0]->y())*(m_curvFeStruc->currentQuadPoint[ig].x() - ptElem[0]->x())/det1
2996 ✗ + (ptElem[0]->x() - ptElem[2]->x())*(m_curvFeStruc->currentQuadPoint[ig].y() - ptElem[0]->y())/det1;
2997
2998 ✗ intPoints[ig][1] = (ptElem[0]->y() - ptElem[1]->y())*(m_curvFeStruc->currentQuadPoint[ig].x() - ptElem[0]->x())/det1
2999 ✗ + (ptElem[1]->x() - ptElem[0]->x())*(m_curvFeStruc->currentQuadPoint[ig].y() - ptElem[0]->y())/det1;
3000
3001 ✗ intPoints[ig][2] = 0.;
3002
3003 }
3004 ✗ } else if(abs(det2) >= abs(det1) && abs(det2) >= abs(det3)){
3005 ✗ for(int ig=0; ig<m_curvFeStruc->numQuadraturePoint(); ++ig) {
3006 ✗ intPoints[ig][0] = (ptElem[2]->x() - ptElem[0]->x())*(m_curvFeStruc->currentQuadPoint[ig].z() - ptElem[0]->z())/det2
3007 ✗ + (ptElem[0]->z() - ptElem[2]->z())*(m_curvFeStruc->currentQuadPoint[ig].x() - ptElem[0]->x())/det2;
3008
3009 ✗ intPoints[ig][1] = (ptElem[0]->x() - ptElem[1]->x())*(m_curvFeStruc->currentQuadPoint[ig].z() - ptElem[0]->z())/det2
3010 ✗ + (ptElem[1]->z() - ptElem[0]->z())*(m_curvFeStruc->currentQuadPoint[ig].x() - ptElem[0]->x())/det2;
3011
3012 ✗ intPoints[ig][2] = 0.;
3013 }
3014 } else {
3015 ✗ for(int ig=0; ig<m_curvFeStruc->numQuadraturePoint(); ++ig) {
3016 ✗ intPoints[ig][0] = (ptElem[2]->z() - ptElem[0]->z())*(m_curvFeStruc->currentQuadPoint[ig].y() - ptElem[0]->y())/det3
3017 ✗ + (ptElem[0]->y() - ptElem[2]->y())*(m_curvFeStruc->currentQuadPoint[ig].z() - ptElem[0]->z())/det3;
3018
3019 ✗ intPoints[ig][1] = (ptElem[0]->z() - ptElem[1]->z())*(m_curvFeStruc->currentQuadPoint[ig].y() - ptElem[0]->y())/det3
3020 ✗ + (ptElem[1]->y() - ptElem[0]->y())*(m_curvFeStruc->currentQuadPoint[ig].z() - ptElem[0]->z())/det3;
3021
3022 ✗ intPoints[ig][2] = 0.;
3023 }
3024 }
3025 }
3026
3027 // update the structure finite element with the element (really needed ? we only want to access the basis
3028 // function directly to evaluate them at the computed integration points).
3029 ✗ m_curvFeStruc->update(0, ptElem);
3030
3031 // Now we can evaluate the basis function of the reference element at the integration points of
3032 // the sub element
3033 ✗ for(int icomp=0; icomp<m_feVel->numCoor(); ++icomp) {
3034 ✗ for(int ig=0; ig<m_curvFeStruc->numQuadraturePoint(); ++ig) {
3035 ✗ m_strucVelQuadPts.val(icomp, ig) = 0;
3036 ✗ for(int idof=0; idof<m_curvFeStruc->numDof(); ++idof) {
3037 ✗ m_strucVelQuadPts.val(icomp, ig) += m_strucVelDofPts.val(icomp, idof) * m_curvFeStruc->refEle().basisFunction().phi(idof, intPoints[ig]);
3038 }
3039 }
3040 }
3041 }
3042
3043 ✗ void LinearProblemNitscheXFEM::initPerElementTypeDarcy(ElementType& eltType, FlagMatrixRHS flagMatrixRHS)
3044 {
3045 IGNORE_UNUSED_ELT_TYPE;
3046 IGNORE_UNUSED_FLAG_MATRIX_RHS;
3047
3048 // Get the curvilinear finite elements
3049 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
3050 ✗ m_curvFePress = m_listCurvilinearFiniteElement[m_iPressure];
3051 ✗ m_curvFeDarcy = m_listCurvilinearFiniteElement[m_iPreDarcy];
3052 }
3053
3054 ✗ void LinearProblemNitscheXFEM::computeElementArrayDarcy(const std::vector<Point*>& elemPoint)
3055 {
3056
3057 // call to the function in linear problem that updates the curvilinear finite element
3058 ✗ if(m_curvFeVel->hasOriginalQuadPoint()) {
3059 ✗ m_curvFeVel->updateMeasNormalContra(0, elemPoint);
3060 ✗ m_curvFeDarcy->updateMeasNormalContra(0, elemPoint);
3061 } else {
3062 ✗ m_curvFeVel->updateBasisAndNormalContra(0, elemPoint);
3063 ✗ m_curvFeDarcy->updateBasisAndNormalContra(0, elemPoint);
3064 }
3065
3066 //if ( FelisceParam::instance().contactDarcy == 2 ) // TODO
3067 // m_curvFeVel->updateMeasNormalTangent(0, elemPoint);
3068 }
3069
3070 ✗ void LinearProblemNitscheXFEM::elementComputeDarcyOnBoundary(const std::vector<Point*>& elemPoint)
3071 {
3072 // Check if we are in the labels where we want to consider the porous media
3073 // No need to reupdate the FEM, already done in computeElementArrayDarcy.
3074
3075 ✗ if ( dimension() == 2 ) {
3076 ✗ const double eps = FelisceParam::instance().epsDarcy;
3077 ✗ const double K = FelisceParam::instance().kDarcy;
3078 ✗ const double alpha = 1.;
3079
3080 // FLUID BOUNDARY TERMS
3081 // + eps /4K uf_j . n vf_i . n
3082 ✗ m_elementMatBD[0]->phi_i_dot_n_phi_j_dot_n(0.25*eps/K, *m_curvFeVel, m_velBlock, m_velBlock, m_velocity->numComponent());
3083
3084 // + pl_j vf_i.n
3085 ✗ m_elementMatBD[0]->psi_j_phi_i_dot_n(1, *m_curvFeVel, *m_curvFeDarcy, m_velBlock, m_darBlock, 0);
3086
3087 // BJS condition
3088 ✗ if( FelisceParam::instance().contactDarcy == 2){
3089 ✗ const double coeff = -alpha/sqrt(eps*K);
3090 ✗ std::vector<double> edgeTan(2);
3091
3092 // tangent vector in 2D
3093 ✗ edgeTan[1] = elemPoint[1]->y() - elemPoint[0]->y();
3094 ✗ edgeTan[0] = elemPoint[1]->x() - elemPoint[0]->x();
3095
3096 ✗ const double tmpTang = 1/sqrt(edgeTan[0]*edgeTan[0] + edgeTan[1]*edgeTan[1]); // TODO
3097 ✗ edgeTan[1] *= tmpTang;
3098 ✗ edgeTan[0] *= tmpTang;
3099
3100 ✗ for(std::size_t idim=0; idim<m_velocity->numComponent(); ++idim)
3101 ✗ m_elementMatBD[0]->phi_i_phi_j(coeff*edgeTan[idim], *m_curvFeVel, m_velBlock+idim, m_velBlock+idim, 1);
3102
3103 // m_elementMatBD[0]->phi_i_dot_t_phi_j_dot_t(coeff, *m_curvFeVel, m_velBlock+idim, m_velBlock+idim, m_velocity->numComponent());
3104 }
3105
3106 // DARCY TERMS
3107 // - uf_j.n ql_i
3108 ✗ m_elementMatBD[0]->phi_i_psi_j_scalar_n(-1, *m_curvFeDarcy, *m_curvFeVel, m_darBlock, m_velBlock, m_velocity->numComponent());
3109
3110 // + eps * Kt grad_t Pl_j . grad_t ql_i
3111 ✗ m_elementMatBD[0]->grad_phi_i_grad_phi_j( eps*K, *m_curvFeDarcy, m_darBlock, m_darBlock, m_presDarcy->numComponent());
3112
3113 } else {
3114 ✗ FEL_ERROR("DarcyForContact model is not yet available for 3D");
3115 }
3116 }
3117
3118 ✗ void LinearProblemNitscheXFEM::m_computeElementaryInterfaceStress(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/, const ElementType eltType, const felInt iel, ElementVector& elemVec)
3119 {
3120
3121 // loop over the sub structure to compute the integrals involving the fluid velocity and pressure
3122 ✗ felInt numSubElement = m_duplicateSupportElements->getNumSubItfEltPerItfElt(iel);
3123 ✗ if ( numSubElement > 0 ) {
3124
3125 // integration points
3126 ✗ m_intPoints.clear();
3127 ✗ m_intPoints.resize(m_curvFeStruc->numQuadraturePoint());
3128
3129 // Assembling of the terms coming from the coupling with the fluid and the Nitsche-XFEM formulation
3130 ✗ std::vector<felInt> vectorIdSupport;
3131
3132 felInt subElemId;
3133 felInt fluidElemId;
3134 ✗ std::vector<double> tmpPt(3);
3135 ✗ std::vector<double> itfNormal(dimension());
3136
3137 ✗ GeometricMeshRegion::ElementType eltTypeFluid = this->mesh()->bagElementTypeDomain()[0];
3138 ✗ const int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltTypeFluid];
3139 ✗ std::vector<Point*> fluidElemPts(numPointsPerElt);
3140 ✗ std::vector<felInt> fluidElemPtsId(numPointsPerElt);
3141
3142 ✗ if( m_skipSmallVolume ){
3143 ✗ m_totalArea = 0;
3144 ✗ m_skippedArea = 0;
3145 }
3146
3147 ✗ updateStructureVel(iel);
3148
3149 // get the normal to this sub element
3150 ✗ m_interfMesh->listElementNormal(eltType, iel).getCoor(itfNormal);
3151
3152 ✗ for(felInt iSubElt = 0; iSubElt < numSubElement; ++iSubElt) {
3153 // get the id of the solid sub element
3154 ✗ subElemId = m_duplicateSupportElements->getSubItfEltIdxItf(iel, iSubElt);
3155
3156 // fill the vector with the points of the sub-element
3157 ✗ for(std::size_t iPt=0; iPt < m_numVerPerSolidElt; ++iPt) {
3158 ✗ m_duplicateSupportElements->getSubItfEltVerCrd(subElemId, iPt, tmpPt);
3159 ✗ *m_subItfPts[iPt] = Point(tmpPt[0], tmpPt[1], tmpPt[2]);
3160 }
3161
3162 ✗ if( m_skipSmallVolume ){
3163 ✗ const double area = compAreaVolume(m_subItfPts);
3164 ✗ m_totalArea += area;
3165 ✗ if ( area <= m_eps*compAreaVolume(elemPoint) ) {
3166 ✗ m_skippedArea += std::abs(area);
3167 ✗ continue;
3168 }
3169 }
3170
3171 // update measure
3172 ✗ m_curvFeStruc->updateSubElementMeasNormal(0, elemPoint, m_subItfPts);
3173
3174 // get the coordinate of the fluid element
3175 ✗ fluidElemId = m_duplicateSupportElements->getMshEltIdxOfSubItfElt(subElemId);
3176
3177 ✗ mesh()->getOneElement(eltTypeFluid, fluidElemId, fluidElemPtsId, fluidElemPts, 0);
3178
3179 // update the element for the fluid finite element
3180 ✗ m_feVel->updateFirstDeriv(0, fluidElemPts);
3181
3182 // compute the integration point for this element
3183 ✗ m_computeIntegrationPoints(fluidElemPts);
3184
3185 // get all the support elements for the fluid
3186 ✗ m_supportDofUnknown[0].getIdElementSupport(fluidElemId, vectorIdSupport);
3187
3188 // Nitsche penalty term (derivative of the displacement)
3189 // This integral is here because the penalty parameter depends of the fluid element
3190 ✗ const double penaltyTerm = m_viscosity * m_nitschePar / m_feVel->diameter();
3191
3192 // stress coming from the nitsche penalty term
3193 // - int_\Sigma gamma.mu/h dotd . w
3194 ✗ elemVec.source(-2.*penaltyTerm, *m_curvFeStruc, m_strucVelDofPts, 0, m_curvFeStruc->numCoor());
3195
3196 ✗ for(std::size_t iSup = 0; iSup < vectorIdSupport.size(); ++iSup) {
3197 // get the velocity at the dof of the structure
3198 ✗ m_seqVelDofPts.setValue(sequentialSolution(), *m_feVel, vectorIdSupport[iSup], m_iVelocity, m_ao, dof());
3199 ✗ m_computeFluidVelOnSubStructure();
3200
3201 // int_\Sgima gamma.mu/h u^n. w
3202 ✗ elemVec.source(penaltyTerm, *m_curvFeStruc, m_elemFieldSolidRHS, 0, m_curvFeStruc->numCoor());
3203 }
3204
3205 // \int sigma(u,p) n . w
3206 ✗ for(std::size_t iSup = 0; iSup < vectorIdSupport.size(); ++iSup) {
3207 // get the velocity and pressure at the dof of the structure
3208 ✗ m_seqVelDofPts.setValue(sequentialSolution(), *m_feVel, vectorIdSupport[iSup], m_iVelocity, m_ao, dof());
3209 ✗ m_seqPreDofPts.setValue(sequentialSolution(), *m_fePres, vectorIdSupport[iSup], m_iPressure, m_ao, dof());
3210
3211 // evaluate stress at quadrature points
3212 ✗ m_computeFluidStress(2.*iSup - 1., itfNormal);
3213 ✗ elemVec.source(-1, *m_curvFeStruc, m_elemFieldSolidRHS, 0, m_curvFeStruc->numCoor());
3214 }
3215 }
3216
3217 ✗ if( m_skipSmallVolume && MpiInfo::rankProc() == 0 ) {
3218 ✗ if ( std::abs(m_skippedArea) > m_eps * std::abs(m_totalArea) )
3219 ✗ std::cout << "In computeElementaryInterfaceStress, Skipped surface is : " << std::fixed << std::setprecision(16) << m_skippedArea / m_totalArea * 100. << " %" << std::defaultfloat << std::endl;
3220 }
3221 }
3222 }
3223
3224 ✗ void LinearProblemNitscheXFEM::m_computeIntegrationPoints(std::vector<Point*>& ptElem)
3225 { // TODO D.C. this function should be moved to finite element class
3226 // compute the integration of the substructure in the reference element of the structure element.
3227 ✗ this->m_curvFeStruc->update(0, m_subItfPts);
3228 ✗ this->m_curvFeStruc->computeCurrentQuadraturePoint();
3229
3230 ✗ TypeShape fluidCellShape = (TypeShape) m_feVel->geoEle().shape().typeShape();
3231
3232 // mapping from the current fluid element to the fluid reference element
3233 ✗ if(fluidCellShape == Quadrilateral) {
3234 ✗ for(int ig=0; ig<this->m_curvFeStruc->numQuadraturePoint(); ++ig) {
3235 ✗ for(int icomp=0; icomp<m_feVel->numCoor(); ++icomp) {
3236 ✗ m_intPoints[ig][icomp] = 2*(this->m_curvFeStruc->currentQuadPoint[ig][icomp] - (*ptElem[0])[icomp])/((*ptElem[2])[icomp] - (*ptElem[0])[icomp]) - 1.;
3237 }
3238 }
3239 ✗ } else if(fluidCellShape == Triangle) {
3240 ✗ double det = (ptElem[2]->y() - ptElem[0]->y())*(ptElem[1]->x() - ptElem[0]->x())
3241 ✗ - (ptElem[2]->x() - ptElem[0]->x())*(ptElem[1]->y() - ptElem[0]->y());
3242
3243 ✗ for(int ig=0; ig<this->m_curvFeStruc->numQuadraturePoint(); ++ig) {
3244 ✗ m_intPoints[ig][0] = (ptElem[2]->y() - ptElem[0]->y())*(this->m_curvFeStruc->currentQuadPoint[ig].x() - ptElem[0]->x())/det + (ptElem[0]->x() - ptElem[2]->x()) * (this->m_curvFeStruc->currentQuadPoint[ig].y() - ptElem[0]->y())/det;
3245 ✗ m_intPoints[ig][1] = (ptElem[0]->y() - ptElem[1]->y())*(this->m_curvFeStruc->currentQuadPoint[ig].x() - ptElem[0]->x())/det + (ptElem[1]->x() - ptElem[0]->x()) * (this->m_curvFeStruc->currentQuadPoint[ig].y() - ptElem[0]->y())/det;
3246 }
3247 ✗ } else if(fluidCellShape == Tetrahedron ) {
3248
3249 ✗ UBlasPermutationMatrix permut(m_feVel->numCoor());
3250 ✗ UBlasMatrix FLU(m_feVel->numCoor(),m_feVel->numCoor());
3251 ✗ UBlasMatrix IdM(m_feVel->numCoor(),m_feVel->numCoor());
3252 ✗ UBlasVector refPoint(m_feVel->numCoor());
3253 ✗ UBlasVector tmpRefPoint(m_feVel->numCoor());
3254 ✗ IdM.assign(boost::numeric::ublas::identity_matrix<double>(FLU.size1()));
3255 //filling the matrix
3256 ✗ for(int indx = 0; indx < m_feVel->numCoor() ; ++indx){
3257 ✗ FLU(0,indx) = ptElem[indx+1]->x() - ptElem[0]->x();
3258 ✗ FLU(1,indx) = ptElem[indx+1]->y() - ptElem[0]->y();
3259 ✗ FLU(2,indx) = ptElem[indx+1]->z() - ptElem[0]->z();
3260 }
3261
3262 ✗ lu_factorize(FLU,permut);
3263 ✗ lu_substitute(FLU,permut,IdM);
3264
3265 ✗ for(int ig=0; ig<this->m_curvFeStruc->numQuadraturePoint(); ++ig){
3266
3267 // filling the rhs
3268 ✗ tmpRefPoint(0) = this->m_curvFeStruc->currentQuadPoint[ig].x() - ptElem[0]->x();
3269 ✗ tmpRefPoint(1) = this->m_curvFeStruc->currentQuadPoint[ig].y() - ptElem[0]->y();
3270 ✗ tmpRefPoint(2) = this->m_curvFeStruc->currentQuadPoint[ig].z() - ptElem[0]->z();
3271
3272 ✗ refPoint = prod(IdM,tmpRefPoint);
3273 ✗ m_intPoints[ig][0] = refPoint[0];
3274 ✗ m_intPoints[ig][1] = refPoint[1];
3275 ✗ m_intPoints[ig][2] = refPoint[2];
3276
3277 ✗ refPoint.clear();
3278 ✗ tmpRefPoint.clear();
3279 }
3280 ✗ } else {
3281 ✗ FEL_ERROR("This shape type is not implemented yet");
3282 }
3283 }
3284
3285 ✗ void LinearProblemNitscheXFEM::m_computeFluidVelOnSubStructure()
3286 {
3287 // Compute the velocity at the integration points of the sub structure
3288 // These points are computed in the function computeIntegrationPoints
3289 ✗ for(int icomp=0; icomp<m_curvFeStruc->numCoor(); ++icomp) {
3290 ✗ for(int ig=0; ig<this->m_curvFeStruc->numQuadraturePoint(); ++ig) {
3291 ✗ m_elemFieldSolidRHS.val(icomp, ig) = 0;
3292 ✗ for(int idof=0; idof<m_feVel->numDof(); ++idof) {
3293 ✗ m_elemFieldSolidRHS.val(icomp, ig) += m_seqVelDofPts.val(icomp, idof) * m_feVel->refEle().basisFunction().phi(idof, m_intPoints[ig]);
3294 }
3295 }
3296 }
3297 }
3298
3299 ✗ void LinearProblemNitscheXFEM::m_computeFluidStress(double sideCoeff, const std::vector<double>& strucNormal)
3300 {
3301 // compute the tensor at integration points of the sub structure
3302 // These points are computed in the function computeIntegrationPoints
3303 ✗ double tmpval1 = 0;
3304
3305 ✗ for(int icomp=0; icomp<m_curvFeStruc->numCoor(); ++icomp) {
3306 ✗ for(int ig=0; ig< m_curvFeStruc->numQuadraturePoint(); ++ig) {
3307 ✗ m_elemFieldSolidRHS.val(icomp, ig) = 0;
3308
3309 // velocity
3310 ✗ for(int jcomp=0; jcomp<m_feVel->numCoor(); ++jcomp) {
3311
3312 ✗ tmpval1 = 0;
3313 ✗ for(int idof=0; idof<m_feVel->numDof(); ++idof)
3314 ✗ tmpval1 += m_seqVelDofPts.val(icomp, idof) * m_feVel->dPhi[ig](jcomp, idof);
3315
3316 ✗ m_elemFieldSolidRHS.val(icomp, ig) += m_viscosity * tmpval1 * sideCoeff * strucNormal[jcomp];
3317 }
3318
3319 // pressure
3320 ✗ for(int idof=0; idof<m_feVel->numDof(); ++idof)
3321 ✗ m_elemFieldSolidRHS.val(icomp, ig) -= sideCoeff * m_seqPreDofPts.val(0, idof) * m_fePres->refEle().basisFunction().phi(idof, m_intPoints[ig]) * strucNormal[icomp];
3322 }
3323 }
3324
3325 // if we use the symmetric formulation for the stress
3326 ✗ if(m_useSymmetricStress) {
3327 ✗ for(int idof=0; idof<m_feVel->numDof(); ++idof) {
3328
3329 ✗ tmpval1 = 0;
3330 ✗ for(int jcomp=0; jcomp<m_feVel->numCoor(); ++jcomp)
3331 ✗ tmpval1 += m_seqVelDofPts.val(jcomp, idof) * strucNormal[jcomp];
3332
3333 ✗ for(int icomp=0; icomp<m_curvFeStruc->numCoor(); ++icomp)
3334 ✗ for(int ig=0; ig<this->m_curvFeStruc->numQuadraturePoint(); ++ig)
3335 ✗ m_elemFieldSolidRHS.val(icomp, ig) += m_viscosity * tmpval1 * sideCoeff * m_feVel->dPhi[ig](icomp, idof);
3336 }
3337 }
3338 }
3339
3340 ✗ double LinearProblemNitscheXFEM::compAreaVolume(const std::vector<Point*>& elemPoint) const
3341 {
3342 double area;
3343 Point normal, edge;
3344 ✗ MathUtilities::CrossProduct(normal.getCoor(), (*elemPoint[1] - *elemPoint[0]).getCoor(), (*elemPoint[2] - *elemPoint[0]).getCoor());
3345
3346 ✗ if(elemPoint.size() == 3){
3347 ✗ area = 0.5 * normal.norm();
3348 }
3349 else {
3350 ✗ edge = (*elemPoint[3]) - (*elemPoint[0]);
3351 ✗ area = std::inner_product(normal.getCoor().begin(), normal.getCoor().end(), edge.getCoor().begin(), 0.) / 6.;
3352 }
3353 ✗ if(area < 0){
3354 ✗ std::cout << " ERROR: Element with negative area" << std::endl;
3355 }
3356 ✗ return area;
3357 }
3358
3359 ✗ void LinearProblemNitscheXFEM::printSkipVolume(bool printSlipVolume)
3360 {
3361
3362 ✗ m_printSkipVolume = printSlipVolume;
3363 }
3364
3365 ✗ void LinearProblemNitscheXFEM::printMeshPartition(int indexTime, double dt) const
3366 {
3367
3368 {
3369 // write case
3370 FILE * pFile;
3371 ✗ std::string fileName, fileNameDir, fileNameCase;
3372
3373 ✗ fileNameDir = FelisceParam::instance().resultDir;
3374 ✗ fileNameCase = fileNameDir+"intersection.case";
3375
3376 ✗ pFile = fopen (fileNameCase.c_str(),"w");
3377 ✗ if (pFile==NULL) {
3378 ✗ std::string command = "mkdir -p " + fileNameDir;
3379 ✗ int ierr = system(command.c_str());
3380 ✗ if( ierr > 0)
3381 ✗ FEL_ERROR("Impossible to write "+fileNameCase+" case file");
3382 ✗ pFile = fopen (fileNameCase.c_str(),"w");
3383 }
3384 ✗ fprintf( pFile,"FORMAT\n");
3385 // fprintf( pFile,"type: ensight\n");
3386 ✗ fprintf( pFile,"type: ensight gold\n");
3387 ✗ fprintf( pFile,"GEOMETRY\n");
3388 // fprintf( pFile,"model: 1 fluid.geo.*****.geo\n");
3389 ✗ fprintf( pFile,"model: 1 fluid.geo\n");
3390 ✗ fprintf( pFile,"VARIABLE\n");
3391 ✗ fprintf( pFile,"scalar per element: 1 intersection intersection.*****.scl\n");
3392 ✗ fprintf( pFile,"TIME\n");
3393 ✗ fprintf( pFile, "time set: %d\n", 1);
3394 ✗ fprintf( pFile,"number of steps: %d \n", indexTime+1);
3395 ✗ fprintf( pFile, "filename start number: %d\n", 0);
3396 ✗ fprintf( pFile, "filename increment: %d\n", 1);
3397 ✗ fprintf( pFile, "time values:\n");
3398
3399 ✗ int count=0;
3400 ✗ for (int i=0; i < indexTime+1; ++i) {
3401 ✗ fprintf( pFile, "%12.5e",i*dt);
3402 ✗ ++count;
3403 ✗ if ( count == 6) {
3404 ✗ fprintf( pFile, "\n" );
3405 ✗ count=0;
3406 }
3407 }
3408 ✗ fclose(pFile);
3409 }
3410
3411 {
3412 felInt currentIelGeo;
3413 ✗ std::vector<felInt> vecSupportCurrent;
3414 GeometricMeshRegion::ElementType eltType;
3415
3416 ✗ std::string fileName = FelisceParam::instance().resultDir + "intersection";
3417 ✗ std::stringstream oss;
3418 ✗ oss << indexTime;
3419 ✗ std::string aux = oss.str();
3420 ✗ if ( indexTime < 10 ) {
3421 ✗ aux = "0000" + aux;
3422 ✗ } else if ( indexTime < 100 ) {
3423 ✗ aux = "000" + aux;
3424 ✗ } else if ( indexTime < 1000 ) {
3425 ✗ aux = "00" + aux;
3426 ✗ } else if ( indexTime < 10000 ) {
3427 ✗ aux = "0" + aux;
3428 }
3429
3430 ✗ fileName = fileName + "." + aux + ".scl";
3431 FILE * pFile;
3432 ✗ pFile = fopen (fileName.c_str(),"w");
3433
3434 ✗ if (pFile==nullptr) {
3435 ✗ FEL_ERROR("Impossible to write " + fileName + " scalar ensight solution");
3436 }
3437
3438 // int count=0;
3439 ✗ fprintf( pFile, "Scalar per element\n");
3440 ✗ fprintf(pFile, "part\n");
3441 ✗ fprintf(pFile, "%10d\n",0);
3442 ✗ if (mesh()->domainDim() == 2) {
3443 ✗ fprintf(pFile, "tria3\n");
3444 }
3445 ✗ else if (mesh()->domainDim() == 3) {
3446 ✗ fprintf(pFile, "tetra4\n");
3447 }
3448
3449 ✗ currentIelGeo = 0;
3450 ✗ for (std::size_t i=0; i<mesh()->bagElementTypeDomain().size(); ++i) {
3451 ✗ eltType = mesh()->bagElementTypeDomain()[i];
3452
3453 // second loop on element
3454 ✗ for (felInt iElm = 0; iElm < mesh()->numElements(eltType); ++iElm) {
3455
3456 ✗ fprintf(pFile,"%d\n", static_cast<int>(m_eltPart[m_currentMesh][currentIelGeo]));
3457
3458 ✗ ++currentIelGeo;
3459 }
3460 }
3461 ✗ fprintf(pFile, "\n");
3462 ✗ fclose(pFile);
3463 }
3464 }
3465
3466 ✗ void LinearProblemNitscheXFEM::writeDuplication(int rank, IO& io, double& time, int iteration) const
3467 {
3468
3469 ✗ std::vector<felInt> vec_id(dimension()+1);
3470 ✗ double *solution = new double[supportDofUnknown(0).listNode().size()];
3471 ✗ for (std::size_t i = 0; i < supportDofUnknown(0).listNode().size(); ++i)
3472 ✗ solution[i] = 0;
3473
3474
3475 ✗ for (std::size_t iElm = 0; iElm < supportDofUnknown(0).vectorIdElementSupport().size(); ++iElm) {
3476 ✗ auto& r_vec = supportDofUnknown(0).vectorIdElementSupport()[iElm];
3477
3478 ✗ std::size_t numDpl = r_vec.size();
3479
3480 ✗ if ( numDpl > 1 ) {
3481 ✗ supportDofUnknown(0).getIdPointElementSupport(r_vec[0], vec_id);
3482
3483 ✗ for (auto node_id = vec_id.begin(); node_id < vec_id.end(); ++node_id)
3484 ✗ solution[*node_id] = 1;
3485
3486 ✗ supportDofUnknown(0).getIdPointElementSupport(r_vec[1], vec_id);
3487
3488 ✗ for (auto node_id = vec_id.begin(); node_id < vec_id.end(); ++node_id)
3489 ✗ solution[*node_id] = -1;
3490 }
3491 }
3492
3493 ✗ io.writeSolution(rank, time, iteration, 0, "duplication", solution, supportDofUnknown(0).listNode().size(), dimension(), true);
3494 }
3495 }
3496