GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemLinearElasticity.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 91 509 17.9%
Branches: 135 1495 9.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:
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/linearProblemLinearElasticity.hpp"
21 #include "DegreeOfFreedom/boundaryCondition.hpp"
22 #include "Core/felisceTransient.hpp"
23
24 namespace felisce {
25 4 LinearProblemLinearElasticity::LinearProblemLinearElasticity():
26
11/22
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 11 taken 4 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 4 times.
✗ Branch 15 not taken.
✓ Branch 21 taken 4 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 4 times.
✗ Branch 25 not taken.
✓ Branch 30 taken 4 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 4 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 4 times.
✗ Branch 37 not taken.
✓ Branch 39 taken 4 times.
✗ Branch 40 not taken.
✓ Branch 42 taken 4 times.
✗ Branch 43 not taken.
 4 LinearProblem("Linear Elasticity",/*numberOfMatrices*/2),m_coupled(false) {
27
28
29 4 std::vector<std::string> listPara; //list
30
31 // The petsc vectors normalField and measStar are necessary for the computation
32 // of the P1-normal field.
33
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listPara.emplace_back("normalField");
34
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listPara.emplace_back("measStar");
35 // Old quantities
36
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listPara.emplace_back("displacementOld");
37
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 if (FelisceParam::instance().timeScheme == 0 ) {
38
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listPara.emplace_back("displacementOlder");
39 } else {
40 ✗ listPara.emplace_back("accelerationCurrent");
41 ✗ listPara.emplace_back("accelerationOld");
42 ✗ listPara.emplace_back("velocityCurrent");
43 ✗ listPara.emplace_back("velocityOld");
44 }
45
46 // Stress coming from external problem
47
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listPara.emplace_back("externalStress");
48
49
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 std::vector<std::string> listSeq(listPara);
50
51
2/2
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 4 times.
 24 for( std::size_t k(0); k<listPara.size();++k) {
52
2/4
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
 20 m_vecs.Init(listPara[k]);
53 }
54
2/2
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 4 times.
 24 for( std::size_t k(0); k<listSeq.size();++k) {
55
2/4
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
 20 m_seqVecs.Init(listSeq[k]);
56 }
57
58
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
 4 if ( FelisceParam::instance().withCVG ) {
59 #ifdef FELISCE_WITH_CVGRAPH
60 ✗ m_cvgDirichletVariable = "cvgraphDisplacement";
61 #endif
62 ✗ m_seqVecs.Init("cvgraphDisplacementOld");
63 ✗ m_seqVecs.Init("cvgraphDisplacement");
64 ✗ m_seqVecs.Init("externalVelocity");
65 ✗ m_seqVecs.Init("cvgraphSTRESS");
66 ✗ m_seqVecs.Init("currentVelocity");
67 }
68
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 m_interfaceLabels=FelisceParam::instance().fsiInterfaceLabel;
69 4 }
70
71 ✗ void LinearProblemLinearElasticity::initPetscVecsForCoupling() {
72 ✗ m_seqVecs.Init("previousFixedPointVelocity");
73 ✗ m_vecs.Init("velocityCurrent");
74 ✗ m_seqVecs.Init("velocityCurrent");
75 ✗ this->initFixedPointAcceleration();
76 }
77
78 void
79 ✗ LinearProblemLinearElasticity::initializeDofBoundaryAndBD2VolMaps() {
80 /// this function builds all the mappings
81 /// to deal with petsc objects defined at the interface
82 ✗ m_dofBD[0/*iBD*/].initialize(this);
83 ✗ std::vector<int> components(3);
84 ✗ for ( std::size_t icomp(0); icomp<(std::size_t)this->dimension(); icomp++) {
85 ✗ this->dofBD(/*iBD*/0).buildListOfBoundaryPetscDofs(this, m_interfaceLabels, /*iunknown*/ 0 , icomp);
86 ✗ components[icomp]=icomp;
87 }
88 ✗ m_dofBD[0/*iBD*/].buildBoundaryVolumeMapping(/*iunknown*/0, components);
89 }
90 4 void LinearProblemLinearElasticity::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) {
91
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 LinearProblem::initialize(mesh,comm,doUseSNES);
92
93 4 m_fstransient = fstransient;
94 // you can either specify the lame coefficients or young and poisson, by default I assume to use ( E, \nu)
95
5/10
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 4 times.
 4 if ( FelisceParam::instance().young < 0 || FelisceParam::instance().poisson < 0 ) {
96 ✗ m_mu = FelisceParam::instance().mu_lame;
97 ✗ m_lambda = FelisceParam::instance().lambda_lame;
98 } else {
99
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 m_mu = FelisceParam::instance().young/2/(1+FelisceParam::instance().poisson);
100
4/8
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
 4 m_lambda = FelisceParam::instance().young*FelisceParam::instance().poisson/(1.+FelisceParam::instance().poisson)/(1.-2.*FelisceParam::instance().poisson);
101 }
102
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_deltaT = FelisceParam::instance().timeStep;
103
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_density = FelisceParam::instance().densitySolid;
104
105
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_betaNM = FelisceParam::instance().gammaNM;
106
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_gammaNM = FelisceParam::instance().betaNM;
107
108
4/8
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
 4 std::cout<<"E: " << FelisceParam::instance().young << std::endl;
109
4/8
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
 4 std::cout<<"nu: "<< FelisceParam::instance().poisson << std::endl;
110
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 4 std::cout<<"mu: "<< m_mu << std::endl;
111
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 4 std::cout<<"lambda: "<< m_lambda << std::endl;
112
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 4 std::cout<<"deltaT: "<< m_deltaT << std::endl;
113
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 4 std::cout<<"density: "<< m_density << std::endl;
114
115 4 std::vector<PhysicalVariable> listVariable;
116 4 std::vector<std::size_t> listNumComp;
117
118
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listVariable.push_back(displacement);
119
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 listNumComp.push_back(this->dimension());
120
121
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 userAddOtherVariables(listVariable,listNumComp);
122
123 //define unknowns of the linear system.
124
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_listUnknown.push_back(displacement);
125
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 this->definePhysicalVariable(listVariable,listNumComp);
126
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 m_iDisplacement = this->listVariable().getVariableIdList(displacement);
127 4 }
128
129 124 void LinearProblemLinearElasticity::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS) {
130 (void) eltType;
131 (void) flagMatrixRHS;
132 124 m_fe = m_listCurrentFiniteElement[m_iDisplacement];
133
134
3/6
✓ Branch 2 taken 124 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 124 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 124 times.
✗ Branch 10 not taken.
 124 m_elemFields["displacementOld"].initialize(DOF_FIELD,*m_fe,this->dimension());
135
1/2
✓ Branch 1 taken 124 times.
✗ Branch 2 not taken.
 124 if (FelisceParam::instance().timeScheme == 0) {
136
3/6
✓ Branch 2 taken 124 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 124 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 124 times.
✗ Branch 10 not taken.
 124 m_elemFields["displacementOlder"].initialize(DOF_FIELD,*m_fe,this->dimension());
137 } else {
138 ✗ m_elemFields["velocityOld"].initialize(DOF_FIELD,*m_fe,this->dimension());
139 ✗ m_elemFields["accelerationOld"].initialize(DOF_FIELD,*m_fe,this->dimension());
140 }
141
142
143 124 }
144 120 void LinearProblemLinearElasticity::initPerElementTypeBoundaryCondition(ElementType& eltType, FlagMatrixRHS flagMatrixRHS) {
145 IGNORE_UNUSED_FLAG_MATRIX_RHS;
146 IGNORE_UNUSED_ELT_TYPE;
147
148 120 m_curvFe = m_listCurvilinearFiniteElement[m_iDisplacement];
149 #ifdef FELISCE_WITH_CVGRAPH
150
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 120 times.
 120 if (FelisceParam::instance().withCVG ) {
151 ✗ if ( slave()->thereIsAtLeastOneRobinCondition() ){
152 ✗ m_robinAux.initialize( DOF_FIELD, *m_curvFe, this->dimension() );
153 }
154 }
155 #endif
156 120 }
157
158 35640 void LinearProblemLinearElasticity::computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/, felInt& iel, FlagMatrixRHS /*flagMatrixRHS*/) {
159
160 (void)iel;
161
162 35640 m_curvFe->updateMeas(0,elemPoint);
163 #ifdef FELISCE_WITH_CVGRAPH
164
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 35640 times.
 35640 if (FelisceParam::instance().withCVG ) {
165 ✗ this->cvgraphNaturalBC(iel);
166 }
167 #endif
168 35640 }
169
170 155403 void LinearProblemLinearElasticity::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
171 (void) elemIdPoint;
172 (void) elemPoint;
173
174 155403 m_fe->updateFirstDeriv(0, elemPoint);
175
176
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 155403 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 155403 if ( flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_matrix ) {
177
2/2
✓ Branch 1 taken 5013 times.
✓ Branch 2 taken 150390 times.
 155403 if ( m_fstransient->iteration == 0 ) {
178 //! Divergence of the stress-tensor
179 5013 m_elementMat[1]->eps_phi_i_eps_phi_j( 2*m_mu, *m_fe, 0, 0, this->dimension());
180 5013 m_elementMat[1]->div_phi_j_div_phi_i( m_lambda, *m_fe, 0, 0, this->dimension());
181
182 //! Mass term for the time scheme
183
1/2
✓ Branch 1 taken 5013 times.
✗ Branch 2 not taken.
 5013 if (FelisceParam::instance().timeScheme == 0 ) {
184 5013 m_elementMat[1]->phi_i_phi_j(m_density/(m_deltaT*m_deltaT),*m_fe,0,0,this->dimension());
185 } else {
186 ✗ m_elementMat[1]->phi_i_phi_j(m_density/(m_deltaT*m_deltaT*m_betaNM),*m_fe,0,0,this->dimension());
187 }
188 }
189 }
190
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 155403 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 155403 if ( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_rhs)) {
191 155403 this->computeRHS(iel);
192 }
193 155403 }
194
195 155403 void LinearProblemLinearElasticity::computeRHS(felInt iel) {
196
1/2
✓ Branch 1 taken 155403 times.
✗ Branch 2 not taken.
 155403 if (FelisceParam::instance().timeScheme == 0 ) {
197 //! We read the values we need
198
5/10
✓ Branch 2 taken 155403 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 155403 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 155403 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 155403 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 155403 times.
✗ Branch 17 not taken.
 155403 m_elemFields.Get("displacementOld") .setValue(m_seqVecs.Get("displacementOld") ,*m_fe,iel,m_iDisplacement,m_ao,dof());
199
5/10
✓ Branch 2 taken 155403 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 155403 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 155403 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 155403 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 155403 times.
✗ Branch 17 not taken.
 155403 m_elemFields.Get("displacementOlder") .setValue(m_seqVecs.Get("displacementOlder") ,*m_fe,iel,m_iDisplacement,m_ao,dof());
200
201 //! We scale the values with the correct coefficients and we sum them into the first elementField
202
3/6
✓ Branch 2 taken 155403 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 155403 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 155403 times.
✗ Branch 9 not taken.
 155403 m_elemFields.Get("displacementOld") .val *= 2./m_deltaT/m_deltaT;
203
3/6
✓ Branch 2 taken 155403 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 155403 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 155403 times.
✗ Branch 9 not taken.
 155403 m_elemFields.Get("displacementOlder") .val *= -1./m_deltaT/m_deltaT;
204
5/10
✓ Branch 2 taken 155403 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 155403 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 155403 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 155403 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 155403 times.
✗ Branch 16 not taken.
 155403 m_elemFields.Get("displacementOld") .val += m_elemFields.Get("displacementOlder").val;
205
206 //! Everything is considered as a source term
207
3/6
✓ Branch 5 taken 155403 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 155403 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 155403 times.
✗ Branch 12 not taken.
 155403 m_elementVector[0]->source(m_density, *m_fe, m_elemFields.Get("displacementOld"), /*iblock*/0,/*numComp*/this->dimension());
208
209 } else {
210 //! We read the values we need
211 ✗ m_elemFields.Get("displacementOld") .setValue(m_seqVecs.Get("displacementOld"), *m_fe,iel,m_iDisplacement,m_ao,dof());
212 ✗ m_elemFields.Get("velocityOld") .setValue(m_seqVecs.Get("velocityOld"), *m_fe,iel,m_iDisplacement,m_ao,dof());
213 ✗ m_elemFields.Get("accelerationOld") .setValue(m_seqVecs.Get("accelerationOld"), *m_fe,iel,m_iDisplacement,m_ao,dof());
214
215 //! We scale the values with the correct coefficients and we sum them into the first elementField
216 ✗ m_elemFields.Get("displacementOld").val += m_deltaT*m_elemFields.Get("velocityOld").val;
217 ✗ m_elemFields.Get("displacementOld").val += (m_deltaT*m_deltaT*(1-2*m_betaNM)/2)*m_elemFields.Get("accelerationOld").val;
218
219 //! Everything is considered as a source term
220 ✗ m_elementVector[0]->source(m_density/(m_betaNM*m_deltaT*m_deltaT), *m_fe, m_elemFields.Get("displacementOld"), /*iblock*/0,/*numComp*/this->dimension());
221 }
222 155403 }
223
224 120 void LinearProblemLinearElasticity::updateCurrentVelocity() {
225
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 120 times.
 120 if ( FelisceParam::instance().timeScheme == 1) {
226 ✗ if ( m_coupled ) {
227 ✗ m_seqVecs.Get("previousFixedPointVelocity").copyFrom(m_seqVecs.Get("velocityCurrent"));
228 }
229 // The current acceleration:
230 ✗ m_vecs.Get("accelerationCurrent").copyFrom(m_vecs.Get("accelerationOld"));
231 ✗ m_vecs.Get("accelerationCurrent").axpbypcz(/*a*/1, /*b*/m_deltaT,/*c*/m_deltaT*m_deltaT/2*(1-2*m_betaNM), /*x*/m_vecs.Get("displacementOld"), /*y*/m_vecs.Get("velocityOld"));
232 ✗ m_vecs.Get("accelerationCurrent").axpby(/*a*/1./m_betaNM/m_deltaT/m_deltaT, /*b*/ -1./m_betaNM/m_deltaT/m_deltaT, /*x*/this->solution());
233 ✗ this->gatherVector(m_vecs.Get("accelerationCurrent"),m_seqVecs.Get("accelerationCurrent"));
234
235 ✗ m_vecs.Get("velocityCurrent").copyFrom(m_vecs.Get("velocityOld"));
236 ✗ m_vecs.Get("velocityCurrent").axpbypcz(m_deltaT*(1-m_gammaNM), m_deltaT*m_gammaNM,/*c*/1, m_vecs.Get("accelerationOld"),m_vecs.Get("accelerationCurrent"));
237 ✗ this->gatherVector(m_vecs.Get("velocityCurrent"),m_seqVecs.Get("velocityCurrent"));
238 } else {
239 // The current displacement is the solution...
240
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 120 times.
 120 if ( m_coupled ) {
241 ✗ m_seqVecs.Get("previousFixedPointVelocity").copyFrom(m_seqVecs.Get("velocityCurrent"));
242 // The current velocity:
243 ✗ m_vecs.Get("velocityCurrent").copyFrom(this->solution());// currentVelocity = currentDisplacement
244 ✗ m_vecs.Get("velocityCurrent").axpby(-1./m_deltaT, 1./m_deltaT, m_vecs.Get("displacementOld")); // currentVelocity = (currentDisplacent - oldDisplacement)/delta T
245 ✗ this->gatherVector(m_vecs.Get("velocityCurrent"),m_seqVecs.Get("velocityCurrent"));
246 }
247 }
248
249
250
251 120 }
252
253 120 void LinearProblemLinearElasticity::updateOldDisplacement() {
254
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 if ( FelisceParam::instance().timeScheme == 0 ) {
255
256
5/10
✓ Branch 2 taken 120 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 120 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 120 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 120 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 120 times.
✗ Branch 16 not taken.
 120 m_vecs.Get("displacementOlder").copyFrom(m_vecs.Get("displacementOld"));
257
5/10
✓ Branch 2 taken 120 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 120 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 120 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 120 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 120 times.
✗ Branch 16 not taken.
 120 this->gatherVector(m_vecs.Get("displacementOlder"),m_seqVecs.Get("displacementOlder"));
258
259 } else {
260
261 ✗ m_vecs.Get("velocityOld").copyFrom(m_vecs.Get("velocityCurrent"));
262 ✗ m_vecs.Get("accelerationOld").copyFrom(m_vecs.Get("accelerationCurrent"));
263
264 ✗ this->gatherVector(m_vecs.Get("velocityOld"),m_seqVecs.Get("velocityOld"));
265 ✗ this->gatherVector(m_vecs.Get("accelerationOld"),m_seqVecs.Get("accelerationOld"));
266
267 }
268
3/6
✓ Branch 2 taken 120 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 120 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 120 times.
✗ Branch 10 not taken.
 120 m_vecs.Get("displacementOld").copyFrom(this->solution());
269
5/10
✓ Branch 2 taken 120 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 120 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 120 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 120 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 120 times.
✗ Branch 16 not taken.
 120 this->gatherVector(m_vecs.Get("displacementOld"),m_seqVecs.Get("displacementOld"));
270 #ifdef FELISCE_WITH_CVGRAPH
271
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 120 times.
 120 if( FelisceParam::instance().withCVG ) { //&& m_readVariable == "VELOCITY"
272 ✗ this->gatherSolution();
273 ✗ m_seqVecs.Get("cvgraphDisplacementOld").copyFrom(this->sequentialSolution());
274 }
275 #endif
276 120 }
277
278 void
279 ✗ LinearProblemLinearElasticity::readStressInterface(std::map<int, std::map<int,std::map<int,int> > > externalMap, const PetscVector& externalSeqStress, felInt iUC0) {
280 ✗ if ( FelisceParam::verbose() > 0 ) {
281 ✗ std::stringstream out;
282 ✗ out<<"reading data at the interface...";
283 ✗ PetscPrintf(MpiInfo::petscComm(), "%s",out.str().c_str());
284 }
285 ✗ felInt iDofGlobScalarPressure(0);
286 ✗ felInt iDofGlobScalarIOP(0);
287 ✗ felInt iDofGlobScalarPressureAus(0);
288 ✗ double valueAus(0.0);
289 ✗ std::vector<felInt> labels = FelisceParam::instance().fsiInterfaceLabel;
290 ✗ for ( std::size_t iComp(0); iComp< (std::size_t) this->dimension(); ++iComp, ++iUC0) {
291 ✗ for ( std::size_t iLab(0); iLab < m_interfaceLabels.size(); ++iLab) {
292 //Contrary to the application PerfectFluid we assume the labels are the same, therefore we do not need a mapBetweenLabels
293 ✗ for(auto it_suppDofIop = externalMap.at(iUC0)[m_interfaceLabels[iLab]].begin(); it_suppDofIop != externalMap.at(iUC0)[m_interfaceLabels[iLab]].end(); it_suppDofIop++) {
294 ✗ iDofGlobScalarIOP = it_suppDofIop->second;
295 ✗ iDofGlobScalarPressure = m_dofBD[0/*iBD*/].getPetscDofs( this->dof().getNumGlobComp(m_iDisplacement,iComp), m_interfaceLabels[iLab], it_suppDofIop->first );
296 ✗ iDofGlobScalarPressureAus = iDofGlobScalarPressure;
297 ✗ AOPetscToApplication(m_ao,/*counter*/1,&iDofGlobScalarPressureAus);
298 ✗ if ( m_dofPart[iDofGlobScalarPressureAus] == MpiInfo::rankProc() ) {
299 ✗ externalSeqStress.getValues(1,&iDofGlobScalarIOP, &valueAus);
300 ✗ m_vecs.Get("externalStress").setValue(iDofGlobScalarPressure, valueAus, INSERT_VALUES);
301 }
302 }
303 }
304 }
305 ✗ m_vecs.Get("externalStress").assembly();
306 ✗ gatherVector(m_vecs.Get("externalStress"), m_seqVecs.Get("externalStress"));
307
308 ✗ if ( FelisceParam::verbose() > 0 ) {
309 ✗ std::stringstream out;
310 ✗ out<<"done"<<std::endl;
311 ✗ PetscPrintf(MpiInfo::petscComm(), "%s",out.str().c_str());
312 }
313 }
314 ✗ double LinearProblemLinearElasticity::computeL2ScalarProduct(std::string l, std::string r) {
315 ✗ m_auxiliaryStrings.clear();
316 ✗ m_auxiliaryStrings.push_back(l);
317 ✗ m_auxiliaryStrings.push_back(r);
318
319 ✗ m_auxiliaryDoubles.clear();
320 ✗ m_auxiliaryDoubles.push_back(0.0);
321 ✗ assemblyLoopBoundaryGeneral(&LinearProblemLinearElasticity::scalarProductComputer,
322 ✗ m_interfaceLabels,
323 &LinearProblemLinearElasticity::initPerET,
324 &LinearProblemLinearElasticity::updateFe);
325 ✗ double result(0);
326 ✗ MPI_Allreduce(&m_auxiliaryDoubles[0],&result,1,MPI_DOUBLE,MPI_SUM,MpiInfo::petscComm());
327 ✗ return result;
328 }
329 ✗ void LinearProblemLinearElasticity::scalarProductComputer(felInt ielSupportDof) {
330 double l,r;
331 ✗ for ( std::size_t iComp(0); iComp < (std::size_t) this->dimension(); ++iComp ) {
332 ✗ for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_curvFe->numDof(); iSurfDof++) {
333 felInt iDofGlobScalar;
334 ✗ dof().loc2glob( ielSupportDof, iSurfDof, m_iDisplacement, iComp, iDofGlobScalar);
335 ✗ AOApplicationToPetsc(m_ao,1,&iDofGlobScalar);
336 ✗ m_seqVecs.Get( m_auxiliaryStrings[0] ).getValues(1,&iDofGlobScalar, &r);
337 ✗ for ( std::size_t jSurfDof(0); jSurfDof < ( std::size_t ) m_curvFe->numDof(); jSurfDof++) {
338 felInt jDofGlobScalar;
339 ✗ dof().loc2glob( ielSupportDof, jSurfDof, m_iDisplacement, iComp, jDofGlobScalar);
340 ✗ AOApplicationToPetsc(m_ao,1,&jDofGlobScalar);
341 ✗ m_seqVecs.Get( m_auxiliaryStrings[1] ).getValues(1,&iDofGlobScalar, &l);
342 ✗ for(int ig=0; ig < m_curvFe->numQuadraturePoint(); ig++) {
343 ✗ m_auxiliaryDoubles[0] += l*r * m_curvFe->phi[ig](iSurfDof) * m_curvFe->phi[ig](jSurfDof) * m_curvFe->weightMeas(ig);
344 }
345 }
346 }
347 }
348 }
349
350 std::pair<double,double>
351 ✗ LinearProblemLinearElasticity::computeTestQuantities() {
352 ✗ std::pair<double,double> pairOfTestQuantities;
353 ✗ m_auxiliaryDoubles.clear();
354 ✗ m_auxiliaryDoubles.resize(4,0.0);
355 ✗ assemblyLoopBoundaryGeneral( &LinearProblemLinearElasticity::testQuantitiesComputer,
356 ✗ m_interfaceLabels,
357 &LinearProblemLinearElasticity::initPerET,
358 &LinearProblemLinearElasticity::updateFe);
359 ✗ std::vector<double> result(4,0.0);
360 ✗ MPI_Allreduce(m_auxiliaryDoubles.data(),result.data(),4,MPI_DOUBLE,MPI_SUM,MpiInfo::petscComm());
361 using std::sqrt;
362 ✗ pairOfTestQuantities.first = std::sqrt(result[0])/( std::sqrt(result[1]) + 1e-12);
363 ✗ pairOfTestQuantities.second = std::sqrt(result[2])/( std::sqrt(result[3]) + 1e-12);
364 ✗ return pairOfTestQuantities;
365 }
366 ✗ void LinearProblemLinearElasticity::initPerET() {
367 ✗ m_curvFe = m_listCurvilinearFiniteElement[m_iDisplacement];
368 }
369 ✗ void LinearProblemLinearElasticity::updateFe(const std::vector<Point*>& elemPoint,const std::vector<int>& /*elemIdPoint*/) {
370 ✗ m_curvFe ->updateMeasNormal(0, elemPoint);
371 }
372
373 void
374 ✗ LinearProblemLinearElasticity::testQuantitiesComputer(felInt ielSupportDof) {
375
376 double uni,unoldi,ciopi,x1i;
377 double value0i,value1i,value2i,value3i;
378
379 double unj,unoldj,ciopj,x1j;
380 double value0j,value1j,value2j,value3j;
381
382 felInt iDofGlobScalar;
383
384 ✗ for ( std::size_t iComp(0); iComp < ( std::size_t )this->dimension(); iComp++) {
385 ✗ for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_curvFe->numDof(); iSurfDof++) {
386
387 ✗ dof().loc2glob( ielSupportDof, iSurfDof, m_iDisplacement, iComp, iDofGlobScalar);
388 ✗ AOApplicationToPetsc(m_ao,1,&iDofGlobScalar);
389
390 ✗ m_seqVecs.Get( "velocityCurrent" ).getValues(1,&iDofGlobScalar, &uni);
391 ✗ m_seqVecs.Get( "previousFixedPointVelocity" ).getValues(1,&iDofGlobScalar, &unoldi);
392 ✗ m_seqVecs.Get( "externalStress" ).getValues(1,&iDofGlobScalar, &ciopi);
393 ✗ m_seqVecs.Get( "x1" ).getValues(1,&iDofGlobScalar, &x1i);
394
395 ✗ value0i = uni-unoldi;
396 ✗ value1i = uni;
397 ✗ value2i = ciopi-x1i;
398 ✗ value3i = ciopi;
399
400
401 ✗ for ( std::size_t jSurfDof(0); jSurfDof < ( std::size_t ) m_curvFe->numDof(); jSurfDof++) {
402 felInt jDofGlobScalar;
403 ✗ dof().loc2glob( ielSupportDof, jSurfDof, m_iDisplacement, iComp, jDofGlobScalar);
404 ✗ AOApplicationToPetsc(m_ao,1,&jDofGlobScalar);
405
406 ✗ m_seqVecs.Get( "velocityCurrent" ).getValues(1,&jDofGlobScalar, &unj);
407 ✗ m_seqVecs.Get( "previousFixedPointVelocity" ).getValues(1,&jDofGlobScalar, &unoldj);
408 ✗ m_seqVecs.Get( "externalStress" ).getValues(1,&jDofGlobScalar, &ciopj);
409 ✗ m_seqVecs.Get( "x1" ).getValues(1,&jDofGlobScalar, &x1j);
410
411 ✗ value0j = unj-unoldj;
412 ✗ value1j = unj;
413 ✗ value2j = ciopj-x1j;
414 ✗ value3j = ciopj;
415
416 ✗ for(int ig=0; ig < m_curvFe->numQuadraturePoint(); ig++) {
417 ✗ double w = m_curvFe->phi[ig](iSurfDof) * m_curvFe->phi[ig](jSurfDof) * m_curvFe->weightMeas(ig);
418 ✗ m_auxiliaryDoubles[0] += value0i*value0j * w;
419 ✗ m_auxiliaryDoubles[1] += value1i*value1j * w;
420 ✗ m_auxiliaryDoubles[2] += value2i*value2j * w;
421 ✗ m_auxiliaryDoubles[3] += value3i*value3j * w;
422 }
423 }
424 }
425 }
426 }
427
428 120 void LinearProblemLinearElasticity::applyEssentialBoundaryConditionDerivedProblem(int rank, FlagMatrixRHS flagMatrixRHS) {
429
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 120 times.
 120 if( FelisceParam::instance().useEssDerivedBoundaryCondition ) {
430 ✗ if ( FelisceParam::verbose() > 2 )
431 ✗ std::cout<<"["<<rank<<"]--Applying zero-Dirichlet BCs on the velocity at the interface between inlet outlet and surface."<<std::endl;
432 ✗ getRings();
433 ✗ double bigValue=1.e20;
434
435 ✗ for(auto itDofBC = m_idDofRings.begin(); itDofBC != m_idDofRings.end(); ++itDofBC) {
436 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix))
437 ✗ matrix(0).setValue(*itDofBC, *itDofBC, bigValue, INSERT_VALUES);
438 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_rhs))
439 ✗ vector().setValue(*itDofBC,0.0, INSERT_VALUES);
440 }
441 }
442 120 }
443
444 #ifdef FELISCE_WITH_CVGRAPH
445 ✗ void LinearProblemLinearElasticity::readData() {
446 ✗ bool dispUpdateNeeded(false);
447 ✗ bool dispDirectly(false);
448 ✗ for ( std::size_t iConn(0); iConn<this->slave()->numConnections(); ++iConn ) {
449 ✗ std::vector<PetscVector> vecs;
450 ✗ for ( std::size_t cVar(0); cVar<this->slave()->numVarToRead(iConn); ++cVar) {
451 ✗ std::string varToRead = this->slave()->readVariable(iConn,cVar);
452 ✗ if ( varToRead == "STRESS" ) {
453 ✗ vecs.push_back(m_seqVecs.Get("cvgraphSTRESS"));
454 ✗ } else if ( varToRead == "VELOCITY" ) {
455 ✗ dispUpdateNeeded = true;
456 ✗ vecs.push_back(m_seqVecs.Get("externalVelocity"));
457 ✗ } else if ( varToRead == "DISPLACEMENT" ) {
458 ✗ dispDirectly = true;
459 ✗ vecs.push_back(m_seqVecs.Get("cvgraphDisplacement"));
460 } else {
461 ✗ std::cout<<"Connection: "<<iConn<<"--variable to read: "<<varToRead<<std::endl;
462 ✗ FEL_ERROR("Error in variable to read");
463 }
464 }
465 ✗ this->slave()->receiveData(vecs,iConn);
466 }
467 ✗ if ( dispUpdateNeeded && dispDirectly ) {
468 ✗ FEL_ERROR("TODO: implement a better readData function");
469 }
470 ✗ if ( dispUpdateNeeded ) {
471 ✗ waxpy(m_seqVecs.Get("cvgraphDisplacement"), m_deltaT, m_seqVecs.Get("externalVelocity"),m_seqVecs.Get("cvgraphDisplacementOld"));
472 }
473 }
474
475 ✗ void LinearProblemLinearElasticity::sendData() {
476 ✗ this->gatherSolution();
477
478 ✗ m_seqVecs.Get("currentVelocity").copyFrom(this->sequentialSolution());// currentVelocity = currentDisplacement
479 ✗ m_seqVecs.Get("currentVelocity").axpby(-1./m_deltaT, 1./m_deltaT, m_seqVecs.Get("displacementOld")); // currentVelocity = (currentDisplacent - oldDisplacement)/delta T
480
481 ✗ for ( std::size_t iConn(0); iConn<this->slave()->numConnections(); ++iConn ) {
482 ✗ std::vector<PetscVector> vecs;
483 ✗ for ( std::size_t cVar(0); cVar<this->slave()->numVarToSend(iConn); ++cVar) {
484 ✗ std::string varToSend = this->slave()->sendVariable(iConn, cVar);
485 ✗ if ( varToSend == "STRESS" ) {
486 ✗ this->prepareResidual(iConn);
487 ✗ vecs.push_back(m_seqStressOut);
488 ✗ } else if ( varToSend == "DISPLACEMENT" ) {
489 ✗ vecs.push_back(this->sequentialSolution());
490 ✗ } else if ( varToSend == "VELOCITY" ) {
491 ✗ vecs.push_back(m_seqVecs.Get("currentVelocity"));
492 } else {
493 ✗ std::cout<<"Connection: "<<iConn<<"--variable to send: "<<varToSend<<std::endl;
494 ✗ FEL_ERROR("Error in variable to send");
495 }
496 }
497 ✗ this->slave()->sendData(vecs,iConn);
498 }
499 }
500 ✗ void LinearProblemLinearElasticity::assembleMassBoundaryAndInitKSP( std::size_t iConn ) {
501 ✗ LinearProblem::assembleMassBoundaryAndInitKSP(&LinearProblemLinearElasticity::massMatrixComputer,
502 ✗ this->slave()->interfaceLabels(iConn),
503 &LinearProblemLinearElasticity::initPerETMass,
504 &LinearProblemLinearElasticity::updateFE,
505 iConn);
506 }
507
508 /*! \brief Function to assemble the mass matrix
509 * Function to be called in the assemblyLoopBoundaryGeneral
510 */
511 ✗ void LinearProblemLinearElasticity::massMatrixComputer(felInt ielSupportDof) {
512 ✗ this->m_elementMatBD[0]->zero();
513 ✗ this->m_elementMatBD[0]->phi_i_phi_j(/*coef*/1.,*m_curvFe,/*iblock*/0,/*iblock*/0,this->dimension());
514 ✗ this->setValueMatrixBD(ielSupportDof);
515 }
516
517 void
518 ✗ LinearProblemLinearElasticity::initPerETMass() {
519 ✗ felInt numDofTotal = 0;
520 //pay attention this numDofTotal is bigger than expected since it counts for all the VARIABLES
521 ✗ for (std::size_t iFe = 0; iFe < this->m_listCurvilinearFiniteElement.size(); iFe++) {//this loop is on variables while it should be on unknown
522 ✗ numDofTotal += this->m_listCurvilinearFiniteElement[iFe]->numDof()*this->m_listVariable[iFe].numComponent();
523 }
524 ✗ m_globPosColumn.resize(numDofTotal); m_globPosRow.resize(numDofTotal); m_matrixValues.resize(numDofTotal*numDofTotal);
525
526 ✗ m_curvFe = this->m_listCurvilinearFiniteElement[ m_iDisplacement ];
527 }
528
529 void
530 ✗ LinearProblemLinearElasticity::updateFE(const std::vector<Point*>& elemPoint, const std::vector<int>&) {
531 ✗ m_curvFe->updateMeasNormal(0, elemPoint);
532 }
533 #endif
534
535 ✗ void LinearProblemLinearElasticity::computeBoundaryStress(const std::vector<int>& labels) {
536 ✗ if ( m_computedStress.isNull() )
537 ✗ m_computedStress.duplicateFrom(this->solution());
538 ✗ if ( m_seqComputedStress.isNull() )
539 ✗ m_seqComputedStress.duplicateFrom(this->sequentialSolution());
540 ✗ m_computedStress.zeroEntries();
541 ✗ this->gatherSolution();
542 ✗ m_grad.resize(this->dimension(),this->dimension());
543 ✗ m_stress.resize(this->dimension());
544 ✗ m_previousVolEltType=GeometricMeshRegion::Nod;
545 ✗ if(m_mesh[m_currentMesh]->statusFaces() == false) {
546 ✗ m_mesh[m_currentMesh]->buildFaces();
547 }
548 ✗ assemblyLoopBoundaryGeneral(&LinearProblemLinearElasticity::boundaryStressComputer,
549 labels,
550 &LinearProblemLinearElasticity::initPerETBoundaryStress,
551 &LinearProblemLinearElasticity::updateFEBoundaryStress);
552 ✗ this->desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_rhs);
553 ✗ m_computedStress.assembly();
554 ✗ this->gatherVector(m_computedStress,m_seqComputedStress);
555 }
556
557 /*! \brief Function to compute boundary stress
558 * Function to be called in the assemblyLoopBoundaryGeneral
559 */
560 ✗ void LinearProblemLinearElasticity::boundaryStressComputer(felInt ielSupportDof) {
561 ✗ std::vector<felInt> bd2vol;
562 ✗ m_fe->identifyLocBDDof(*m_curvFe,bd2vol);
563 ✗ m_elemFieldStress.val*=0.0;
564 ✗ for (int hbdDof(0); hbdDof<m_curvFe->numDof(); hbdDof++) {
565 ✗ m_grad = prod(m_fe->dPhiOnDofs[bd2vol.at(hbdDof)],trans(m_elemFieldDisplacementVolumeDofs.val));
566 ✗ m_stress = m_mu*( prod(m_grad,m_curvFe->normal[0]) + prod(trans(m_grad),m_curvFe->normal[0]) ) + m_lambda*Tools::trace(m_grad)*m_curvFe->normal[0];
567 ✗ for (int jcomp(0); jcomp<this->dimension(); jcomp++) {
568 ✗ m_elemFieldStress.val(jcomp,hbdDof) = m_stress(jcomp);
569 }
570 }
571 ✗ m_elementVectorBD[0]->zero();
572 ✗ m_elementVectorBD[0]->source(-1.0, *m_curvFe, m_elemFieldStress, 0, this->dimension());
573 ✗ this->setValueCustomVectorBD(ielSupportDof, ADD_VALUES, m_computedStress );
574 }
575 void
576 ✗ LinearProblemLinearElasticity::initPerETBoundaryStress() {
577 ✗ this->desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS::only_rhs);
578 ✗ this->allocateArrayForAssembleMatrixRHSBD(FlagMatrixRHS::only_rhs);
579 ✗ m_curvFe = m_listCurvilinearFiniteElement[m_iDisplacement];
580 ✗ m_elemFieldStress.initialize(DOF_FIELD, *m_curvFe,this->dimension());
581 }
582 void
583 ✗ LinearProblemLinearElasticity::updateFEBoundaryStress(const std::vector<Point*>& elemPoint, const std::vector<int>& elemIdPoint) {
584 ✗ std::vector<Point*> elemPointVol;
585 ✗ std::vector<felInt> elemIdPointVol;
586 ✗ int idLocFace = -1;
587 ✗ felInt idElemVol = -1;
588 ✗ m_mesh[m_currentMesh]->getElementFromBdElem(elemIdPoint,elemIdPointVol,elemPointVol,idLocFace,idElemVol);
589 int dummy;
590 ✗ m_mesh[m_currentMesh]->getTypeElemFromIdElem(idElemVol, m_volEltType, dummy);
591 ✗ if (m_volEltType != m_previousVolEltType ) {
592 ✗ this->defineFiniteElement(m_volEltType);
593 ✗ this->initElementArray();
594 ✗ m_previousVolEltType=m_volEltType;
595 ✗ m_fe = m_listCurrentFiniteElement[m_iDisplacement];
596 ✗ m_elemFieldDisplacementVolumeDofs.initialize(DOF_FIELD, *m_fe,this->dimension());
597 ✗ m_fe->allocateOnDofsDataStructures();
598 }
599 ✗ m_fe->updateFirstDerivOnDofs(0, elemPointVol);
600 ✗ m_curvFe->updateMeasNormal(0,elemPoint);
601 ✗ m_elemFieldDisplacementVolumeDofs.setValue(this->sequentialSolution(),*m_fe,idElemVol,m_iDisplacement,m_ao,dof());
602 }
603
604 /***********************************************************************************/
605 /***********************************************************************************/
606
607 ✗ void LinearProblemLinearElasticity::computeNormalField(const std::vector<int> & labels, felInt idVecVariable)
608 {
609 ✗ m_auxiliaryInts.resize(1);
610 ✗ m_auxiliaryInts[0]=idVecVariable;//for instance m_iVelocity;
611 /*!
612 For each dof d of the surface the sum: \f$\sum_{K} m(K)\vec n(K)\f$ is computed.
613 - K are all the boundary elements that share the dof d.
614 - m(K) is the measure of the element
615 - n(K) is the P0 normal defined on the element
616 */
617 ✗ this->assemblyLoopBoundaryGeneral(&LinearProblemLinearElasticity::normalFieldComputer,labels,&LinearProblemLinearElasticity::initPerETNormVel,&LinearProblemLinearElasticity::updateFeNormVel);
618
619 ✗ m_vecs.Get("normalField").assembly();
620 ✗ m_vecs.Get("measStar").assembly();
621
622 /// Than the vector is divided by \f$\sum_{K} m(K)\f$ to obtain the P1 field.
623 ✗ pointwiseDivide(m_vecs.Get("normalField"), m_vecs.Get("normalField"), m_vecs.Get("measStar"));
624 ✗ gatherVector( m_vecs.Get("normalField"), m_seqVecs.Get("normalField") );
625
626 /// using normalizeComputer we normalize the result.
627 ✗ this->assemblyLoopBoundaryGeneral(&LinearProblemLinearElasticity::normalizeComputer,labels,&LinearProblemLinearElasticity::initPerETNormVel,&LinearProblemLinearElasticity::updateFeNormVel);
628
629 ✗ m_vecs.Get("normalField").assembly();
630 ✗ gatherVector( m_vecs.Get("normalField"), m_seqVecs.Get("normalField") );
631 }
632
633 /***********************************************************************************/
634 /***********************************************************************************/
635
636 ✗ void LinearProblemLinearElasticity::normalFieldComputer(felInt ielSupportDof)
637 {
638 double value;
639 ✗ for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numDof(); iSurfDof++) {
640 ✗ double feMeas(m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->measure());
641 ✗ for ( std::size_t iCoor(0); iCoor < ( std::size_t ) m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numCoor(); iCoor++) {
642 felInt iDofGlob;
643 ✗ dof().loc2glob( ielSupportDof, iSurfDof, m_auxiliaryInts[0], iCoor, iDofGlob);
644 ✗ AOApplicationToPetsc(m_ao,1,&iDofGlob);
645 ✗ value=0.0;
646 ✗ for ( std::size_t iQuadBd(0.0) ; iQuadBd < ( std::size_t )m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numQuadraturePoint(); iQuadBd++) {
647 ✗ value += m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->weightMeas( iQuadBd ) * m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->normal[iQuadBd](iCoor);
648 }
649 ✗ m_vecs.Get("normalField").setValue(iDofGlob, value, ADD_VALUES);
650 ✗ m_vecs.Get("measStar").setValue(iDofGlob, feMeas, ADD_VALUES);
651 }
652 }
653 }
654
655 /***********************************************************************************/
656 /***********************************************************************************/
657
658 ✗ void LinearProblemLinearElasticity::updateFeNormVel(const std::vector<Point*>& elemPoint,const std::vector<int>& /*elemIdPoint*/)
659 {
660 ✗ m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->updateMeasNormal(0,elemPoint);
661 }
662
663 /***********************************************************************************/
664 /***********************************************************************************/
665
666 ✗ void LinearProblemLinearElasticity::normalizeComputer(felInt ielSupportDof)
667 {
668 ✗ const int numComp = m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numCoor( );
669 ✗ const int numDof = m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numDof( );
670
671 ✗ std::vector< double > norm(numDof,0.0);
672 ✗ std::vector< std::vector< double > > NN(numDof, std::vector< double >(numComp, 0.0));
673 ✗ for ( int idof = 0; idof < numDof; idof++ ) {
674 ✗ for ( int icomp = 0; icomp < numComp; icomp++ ) {
675 int iDofGlob;
676 double value;
677 ✗ dof().loc2glob( ielSupportDof, idof, m_auxiliaryInts[0], icomp, iDofGlob);
678 ✗ AOApplicationToPetsc(m_ao, 1, &iDofGlob);
679 ✗ m_seqVecs.Get("normalField").getValues(1, &iDofGlob, &value);
680 ✗ norm[idof] += value*value;
681 ✗ NN[idof][icomp] = value;
682 }
683 ✗ norm[idof] = std::sqrt( norm[idof] );
684 ✗ for ( int icomp = 0; icomp < numComp; icomp++ )
685 ✗ NN[idof][icomp] = NN[idof][icomp]/norm[idof];
686 }
687 ✗ for ( int idof = 0; idof < numDof; idof++ ) {
688 ✗ for ( int icomp = 0; icomp < numComp; icomp++ ) {
689 int iDofGlob;
690 ✗ dof().loc2glob( ielSupportDof, idof, m_auxiliaryInts[0], icomp, iDofGlob);
691 ✗ AOApplicationToPetsc(m_ao, 1, &iDofGlob);
692 ✗ m_vecs.Get("normalField").setValue(iDofGlob, NN[idof][icomp], INSERT_VALUES);
693 }
694 }
695 }
696
697 /***********************************************************************************/
698 /***********************************************************************************/
699
700 ✗ void LinearProblemLinearElasticity::initFixedPointAcceleration()
701 {
702 ✗ m_accMethod = accelerationType(FelisceParam::instance(this->instanceIndex()).accelerationMethod);
703 ✗ m_omegaAcceleration = FelisceParam::instance(this->instanceIndex()).omegaAcceleration;
704 ✗ m_seqVecs.Init("x1");
705 ✗ switch ( m_accMethod ) {
706 ✗ case FixedPoint:
707 case Relaxation:
708 ✗ break;
709 ✗ case Aitken:
710 case IronsTuck:
711 ✗ m_seqVecs.Init("x0");
712 ✗ m_seqVecs.Init("fx0");
713 ✗ m_seqVecs.Init("fx1");
714 ✗ break;
715 }
716 }
717
718 /***********************************************************************************/
719 /***********************************************************************************/
720
721 ✗ void LinearProblemLinearElasticity::accelerationPreStep(PetscVector& seqCurrentInput)
722 {
723 // pre and post with respect to reading the data at the interface
724
725 // method explanation in accelerationPostStep
726 ✗ switch(m_accMethod) {
727 ✗ case FixedPoint:
728 case Relaxation:
729 ✗ break;
730 ✗ case Aitken:
731 case IronsTuck:
732 // in these two cases we have to move forward x0 and x1
733 ✗ m_seqVecs.Get("x0" ).copyFrom(m_seqVecs.Get("x1" ));
734 ✗ m_seqVecs.Get("fx0").copyFrom(m_seqVecs.Get("fx1"));
735 ✗ break;
736 }
737 // in any case we need to the save the current value of iop which was the last point
738 // where we have have evaluated the functional F.
739 ✗ m_seqVecs.Get("x1").copyFrom(seqCurrentInput);
740 // The following evaluation of the system will
741 // be saved in fx1 in the function acceleration
742 }
743
744 /***********************************************************************************/
745 /***********************************************************************************/
746
747 ✗ void LinearProblemLinearElasticity::accelerationPostStep( PetscVector& seqCurrentInput, int nbOfCurrentIteration )
748 {
749 // We are solving a fixed point problem.
750
751 // here an example for IOPcouplModel:
752 // the input to the coupled system is currentiop
753 // then you solve NS you compute the displacement
754 // than you solve PF given the displacement and the current iop
755 // the solution of PF at the boarder is the outcome of function.
756 // P_interface = F ( currentIop )
757 // an evaluation of F implies the solution of a NS pb + a PF pb.
758 // P_interface will be modified by this function and reused for a new iteration.
759
760 // Four (4) different methods:
761 // 0. fixed point iterations
762 // default method: we simly iterate F over the same value.
763 // new x1=fx1;
764 // 1. relaxation:
765 // new x1= w * fx1 + (1 - w) x1
766 // w=1 ==> back to fixed point iterations
767 // 2. Aitken acceleration
768 // new x1= w * fx1 + (1 - w) x1, but w is a funcion of ( x0, fx0, x1, fx1 )
769 // 3. Irons-Tuck: the same as Aitken, should be more "stable"
770 ✗ switch(m_accMethod) {
771 ✗ case FixedPoint:
772 case Relaxation:
773 // there is no need to save explicitly fx1 as done in aitken and irons tuck
774 // since it is stored in currentIop and it will be overridden with the correct value.
775 ✗ break;
776 ✗ case Aitken:
777 {
778 // at this point we just read the results from the PF problem
779 // this has been stored in the petscvec currentIop, since in a
780 // fixed point approach it is already the new input (new x1 = currentIop = fx1).
781 // in all the other approaches this is just the lates outcome
782 // of the F.
783 ✗ m_seqVecs.Get("fx1").copyFrom( seqCurrentInput );
784
785 ✗ if ( nbOfCurrentIteration == 1 /*first iteration: w=1*/ ) {
786 ✗ m_omegaAcceleration = 1;
787 } else {
788 // we compute the value for the omega parameter
789 ✗ m_seqVecs.Init("dx"); // map element initialization
790 ✗ m_seqVecs.Get("dx").duplicateFrom(m_seqVecs.Get("x1")); // structure of dx copied from x1
791 ✗ m_seqVecs.Get("dx").copyFrom(m_seqVecs.Get("x1")); // dx = x1
792 ✗ m_seqVecs.Get("dx").axpy(-1.,m_seqVecs.Get("x0")); // dx = x1 - x0
793
794 ✗ m_seqVecs.Init("tmp"); // map element initialization
795 ✗ m_seqVecs.Get("tmp").duplicateFrom(m_seqVecs.Get("fx0")); // structure of tmp copied from fx0
796 ✗ m_seqVecs.Get("tmp").copyFrom(m_seqVecs.Get("fx0")); // tmp = fx0
797 ✗ m_seqVecs.Get("tmp").axpy(-1.,m_seqVecs.Get("fx1")); // tmp = fx0 - fx1
798 ✗ m_seqVecs.Get("tmp").axpy(1.,m_seqVecs.Get("dx")); // tmp = fx0 - fx1 + x1 - x0
799
800 ✗ const double omega = computeL2ScalarProduct("dx" ,"tmp");
801 ✗ const double xxnorm = computeL2ScalarProduct("tmp","tmp");
802 ✗ if ( Tools::equal( xxnorm, 0.0 ) )
803 ✗ m_omegaAcceleration = 1;
804 else
805 ✗ m_omegaAcceleration = omega / xxnorm;
806
807 // we destroy those auxiliary vectors
808 ✗ m_seqVecs.Get("tmp").destroy();
809 ✗ m_seqVecs.erase("tmp");
810 ✗ m_seqVecs.Get("dx").destroy();
811 ✗ m_seqVecs.erase("dx");
812
813 }
814 }
815 ✗ break;
816 ✗ case IronsTuck:
817 {
818 ✗ m_seqVecs.Get("fx1").copyFrom(seqCurrentInput);
819
820 ✗ if ( nbOfCurrentIteration == 1 /*first iteration: w=1*/ ) {
821 ✗ m_omegaAcceleration = 1;
822 } else {
823
824 // we compute the value for the omega parameter
825 ✗ m_seqVecs.Init("d1"); // map element initialization
826 ✗ m_seqVecs.Get("d1").duplicateFrom(m_seqVecs.Get("x1")); // structure of d1 copied from x1
827 ✗ m_seqVecs.Get("d1").copyFrom(m_seqVecs.Get("x1")); // d1 = x1
828 ✗ m_seqVecs.Get("d1").axpy(-1.,m_seqVecs.Get("fx1")); // d1 = x1 - fx1
829
830 ✗ m_seqVecs.Init("tmp"); // map element initialization
831 ✗ m_seqVecs.Get("tmp").duplicateFrom(m_seqVecs.Get("x0")); // structure of tmp copied from x0
832 ✗ m_seqVecs.Get("tmp").copyFrom(m_seqVecs.Get("x0")); // tmp = x0
833 ✗ m_seqVecs.Get("tmp").axpy(-1.,m_seqVecs.Get("fx0")); // tmp = x0 - fx0
834 ✗ m_seqVecs.Get("tmp").axpy(-1.,m_seqVecs.Get("d1")); // tmp = x0 - fx0 - (x1 - fx1 )
835
836 ✗ const double xxnorm ( computeL2ScalarProduct("tmp","tmp") );
837 ✗ const double omega ( computeL2ScalarProduct("d1","tmp") );
838 ✗ if ( Tools::equal( xxnorm, 0.0 ) ) {
839 ✗ m_omegaAcceleration = 1;
840 } else {
841 ✗ m_omegaAcceleration = ( 1 + omega / xxnorm ) * m_omegaAcceleration;
842 }
843
844 // we destroy those auxiliary vectors
845 ✗ m_seqVecs.Get("tmp").destroy();
846 ✗ m_seqVecs.Get("d1").destroy();
847 ✗ m_seqVecs.erase("d1");
848 ✗ m_seqVecs.erase("tmp");
849 }
850 }
851 ✗ break;
852 }
853
854 ✗ if ( FelisceParam::verbose() > 1 ) {
855 ✗ std::stringstream msg;
856 ✗ msg<<"----> using omega = "<<m_omegaAcceleration<<std::endl;
857 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
858 }
859
860 //for the fixed point this step is not needed
861 ✗ if ( m_accMethod != FixedPoint ) {
862 ✗ std::stringstream msg;
863 ✗ msg<<" updating fixed point input "<<std::endl;
864 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
865 // update of currentIop, i.e. the new x1, the input for the next evaluation of F
866 ✗ seqCurrentInput.axpby( 1. - m_omegaAcceleration, m_omegaAcceleration, m_seqVecs.Get("x1") );
867 }
868 }
869
870 /***********************************************************************************/
871 /***********************************************************************************/
872
873 ✗ void LinearProblemLinearElasticity::sumOnBoundary(PetscVector& v, PetscVector& b, const DofBoundary& dofBD)
874 {
875 ✗ felInt numLocalDofInterface = dofBD.numLocalDofInterface();
876 ✗ std::vector<double> tmp(numLocalDofInterface);
877 ✗ b.getValues(numLocalDofInterface, dofBD.loc2PetscVolPtr(), tmp.data() );
878 ✗ v.setValues(numLocalDofInterface, dofBD.loc2PetscVolPtr(), tmp.data(),ADD_VALUES);
879 ✗ v.assembly();
880 }
881
882 /***********************************************************************************/
883 /***********************************************************************************/
884
885 #ifdef FELISCE_WITH_CVGRAPH
886 void
887 ✗ LinearProblemLinearElasticity::cvgraphNaturalBC(felInt iel) {
888 BoundaryCondition* BC;
889 // Loop over the cvgraph connections
890 ✗ for ( std::size_t iConn(0); iConn<this->slave()->numConnections(); ++iConn ) {
891 // Extract the labels of the current connection
892 ✗ std::vector<int> labels = this->slave()->interfaceLabels(iConn);
893 // Verify if the current label is part of the interface of this connection
894 ✗ if ( std::find( labels.begin(), labels.end(), m_currentLabel ) != labels.end() ) {
895 ✗ switch ( this->slave()->numVarToRead(iConn) ) {
896 ✗ case 1: // only one var to read: neumann interface BC
897 // Verify that we have to apply the stress in its strong form.
898 ✗ if ( this->slave()->readVariable(iConn,/*cVar*/0) == "STRESS" && this->slave()->residualInterpolationType(iConn) == 0 ) {
899 // Look through the different neumann conditions to find the correct one.
900 ✗ for (std::size_t iNeumann=0; iNeumann<m_boundaryConditionList.numNeumannBoundaryCondition(); iNeumann++ ) {
901 ✗ BC = m_boundaryConditionList.Neumann(iNeumann);//get the the BC
902 ✗ if ( std::find( BC->listLabel().begin(), BC->listLabel().end(),m_currentLabel) != BC->listLabel().end() ) {
903 //compute the force term.
904 ✗ m_elemFieldNeumann[iNeumann].setValue( m_seqVecs.Get("cvgraphSTRESS"), *m_curvFe, iel, m_iDisplacement, m_ao, dof());
905 }
906 }
907 }
908 ✗ break;
909 ✗ case 2: //two variables to read: robin interface BC
910 // Look through the different robin conditions to find the correct one.
911 ✗ for (std::size_t iRobin=0; iRobin<m_boundaryConditionList.numRobinBoundaryCondition(); iRobin++ ) {
912 ✗ BC = m_boundaryConditionList.Robin(iRobin);
913 ✗ if ( std::find( BC->listLabel().begin(), BC->listLabel().end(),m_currentLabel ) != BC->listLabel().end() ) {
914 ✗ if ( this->slave()->residualInterpolationType(iConn) == 0 ) {
915 ✗ m_elemFieldRobin[iRobin].setValue( m_seqVecs.Get("cvgraphSTRESS"), *m_curvFe , iel, m_iDisplacement, m_ao, dof());
916 } else {
917 ✗ m_elemFieldRobin[iRobin].val *=0 ;
918 }
919 ✗ m_robinAux.setValue(m_seqVecs.Get("cvgraphDisplacement"),*m_curvFe, iel, m_iDisplacement, m_ao, dof());
920 ✗ m_elemFieldRobin[iRobin].val += FelisceParam::instance().alphaRobin[iRobin]*m_robinAux.val;
921 }
922 }
923 ✗ break;
924 ✗ default:
925 ✗ FEL_ERROR("Three variables to read not yet implemented for this linear problem.")
926 ✗ break;
927 }
928 }
929 }
930 }
931 #endif
932 }
933