Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
AssembleFEGeneralizedNewtonian_def.hpp
1#ifndef AssembleFEGeneralizedNewtonian_DEF_hpp
2#define AssembleFEGeneralizedNewtonian_DEF_hpp
3
4#include "AssembleFENavierStokes_decl.hpp"
5
6namespace FEDD
7{
8 // All important things are so far defined in AssembleFENavierStokes. Please check there.
9 /* Interesting paper to generalized Newtonian fluids
10 @article{Poole2023,
11 author = {Poole, Robert J.},
12 doi = {10.1016/j.jnnfm.2023.105106},
13 issn = {03770257},
14 journal = {Journal of Non-Newtonian Fluid Mechanics},
15 keywords = {Constitutive modelling,Flow-type,Generalised Newtonian fluids,Inelastic},
16 number = {August},
17 pages = {105106},
18 publisher = {Elsevier B.V.},
19 title = {{Inelastic and flow-type parameter models for non-Newtonian fluids}},
20 url = {https://doi.org/10.1016/j.jnnfm.2023.105106},
21 volume = {320},
22 year = {2023}
23 }
24
25 */
26
27 /* *******************************************************************************
28 This class is for Generalized-Newtonian fluids, where we consider that the viscosity is non-constant. Because the viscosity is no longer constant, the conventional formulation with the Laplacian term cannot be considered
29 (although there is a generalized Laplacian version of the equation, see "On outflow boundary conditions in finite element simulations of non-Newtonian internal flow" 2021).
30 Instead, we use the stress-divergence formulation of the momentum equation and derive from that the element-wise entrie
31 ******************************************************************************* */
32 template <class SC, class LO, class GO, class NO>
33 AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::AssembleFEGeneralizedNewtonian(int flag, vec2D_dbl_Type nodesRefConfig, ParameterListPtr_Type params, tuple_disk_vec_ptr_Type tuple) : AssembleFENavierStokes<SC, LO, GO, NO>(flag, nodesRefConfig, params, tuple)
34 {
35
37 dofsElementViscosity_ = this->dofsPressure_ * this->numNodesVelocity_; // So it is a scalar quantity but as it depend on the velocity it is defined at the nodes of the velocity
38 this->constOutputField_ = vec_dbl_Type(dofsElementViscosity_);
39
40 // Reading through parameterlist
41 shearThinningModel = params->sublist("Material").get("ShearThinningModel", "");
42 // New: We have to check which material model we use
43 if (shearThinningModel == "Carreau-Yasuda")
44 {
45 Teuchos::RCP<CarreauYasuda<SC, LO, GO, NO>> viscosityModelSpecific(new CarreauYasuda<SC, LO, GO, NO>(params));
46 viscosityModel = viscosityModelSpecific;
47 }
48 else if (shearThinningModel == "Power-Law")
49 {
50 Teuchos::RCP<PowerLaw<SC, LO, GO, NO>> viscosityModelSpecific(new PowerLaw<SC, LO, GO, NO>(params));
51 viscosityModel = viscosityModelSpecific;
52 }
53 else if (shearThinningModel == "Dimless-Carreau")
54 {
55 Teuchos::RCP<Dimless_Carreau<SC, LO, GO, NO>> viscosityModelSpecific(new Dimless_Carreau<SC, LO, GO, NO>(params));
56 viscosityModel = viscosityModelSpecific;
57 }
58 else
59 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "No specific implementation for your material model request. Valid are:Carreau-Yasuda, Power-Law, Dimless-Carreau");
60
61 }
62
63 template <class SC, class LO, class GO, class NO>
64 void AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::assembleJacobian()
65 {
66
67 // For nonlinear generalized newtonian stress tensor part
68 SmallMatrixPtr_Type elementMatrixN = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
69 SmallMatrixPtr_Type elementMatrixW = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
70
71 // For nonlinear convection
72 SmallMatrixPtr_Type elementMatrixNC = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
73 SmallMatrixPtr_Type elementMatrixWC = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
74
75 // In the first iteration step we initialize the constant matrices
76 // So in the case of a newtonian fluid we would have the matrix A with the contributions of the Laplacian term
77 // and the matrix B with the mixed-pressure terms. Latter one exists also in the generlized-newtonian case
78 if (this->newtonStep_ == 0)
79 {
80 SmallMatrixPtr_Type elementMatrixA = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
81 SmallMatrixPtr_Type elementMatrixB = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
82
83 this->constantMatrix_.reset(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
84 // Construct the matrix B from FE formulation - as it is equivalent to Newtonian case we call the same function
85 this->assemblyDivAndDivT(elementMatrixB); // For Matrix B
86 elementMatrixB->scale(-1.);
87 this->constantMatrix_->add((*elementMatrixB), (*this->constantMatrix_));
88 }
89
90 // The other element matrices are not constant so we have to update them in each step
91 // As the stress tensor term, considering a generalized-newtonian constitutive equation, is nonlinear we add its contribution here
92
93 // ANB is the FixedPoint Formulation which was named for newtonian fluids.
94 // Matrix A (Laplacian term (here not occuring)), Matrix B for div-Pressure Part, Matrix N for nonlinear parts -
95
96 this->ANB_.reset(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_)); // A + B + N
97 this->ANB_->add(((*this->constantMatrix_)), ((*this->ANB_)));
98
99 // Nonlinear advection term \rho (u \cdot \nabla) u
100 // As this class is derived from NavierStokes class we can call already implemented function
101 //*************** ADVECTION************************
102 this->assemblyAdvection(elementMatrixNC);
103 elementMatrixNC->scale(this->density_);
104 this->ANB_->add((*elementMatrixNC), ((*this->ANB_)));
105
106 // For a generalized-newtonian fluid we add additional element matrix and fill it with specific contribution
107 // Remember that this term is based on the stress-divergence formulation of the momentum equation
108 // \nabla \dot \tau with \tau=\eta(\gammaDot)(\nabla u + (\nabla u)^T)
109 //*************** STRESS TENSOR************************
110 this->assemblyStress(elementMatrixN);
111 this->ANB_->add((*elementMatrixN), ((*this->ANB_)));
112
113 this->jacobian_.reset(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
114 this->jacobian_->add((*this->ANB_), (*this->jacobian_));
115
116 // If linearization is not FixdPoint (so NOX or Newton) we add the derivative to the Jacobian matrix. Otherwise the FixedPoint formulation becomes the jacobian.
117 if (this->linearization_ != "FixedPoint")
118 {
119
120 this->assemblyStressDev(elementMatrixW); // shear stress tensor
121 this->assemblyAdvectionInU(elementMatrixWC); // convection
122 elementMatrixWC->scale(this->density_);
123
124 this->jacobian_->add((*elementMatrixW), (*this->jacobian_));
125 this->jacobian_->add((*elementMatrixWC), (*this->jacobian_)); // int add(SmallMatrix<T> &bMat, SmallMatrix<T> &cMat); //this+B=C elementMatrix + constantMatrix_;
126 }
127
128 //*************** BOUNDARY TERM *******************************
129 /* Because we have stress-divergence form of Navier-Stokes equations in the non-newtonian case
130 we have to add a extra boundary term to get the same outflow boundary condition as in the conventional formulation with
131 the laplacian operator in the equations due to the fact that in the stress-divergence formulation the
132 natural boundary condition is different
133 We have to check whether it is an element which has edges (2D) / surfaces (3D) corresponding to an Outflow Neumann boundary
134 Then we have to compute contribution
135 @ToDo Add if element normal computation is integrated
136 */
137 }
138
139
140
141 template <class SC, class LO, class GO, class NO>
142 void AssembleFEGeneralizedNewtonian<SC,LO,GO,NO>::assembleFixedPoint() {
143
144 SmallMatrixPtr_Type elementMatrixN = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
145 SmallMatrixPtr_Type elementMatrixNC = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
146
147 if(this->newtonStep_ ==0){
148 SmallMatrixPtr_Type elementMatrixB = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
149
150 this->constantMatrix_.reset(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
151 // Construct the matrix B from FE formulation - as it is equivalent to Newtonian case we call the same function
152 this->assemblyDivAndDivT(elementMatrixB); // For Matrix B
153 elementMatrixB->scale(-1.);
154 this->constantMatrix_->add((*elementMatrixB), (*this->constantMatrix_));
155 }
156
157 this->ANB_.reset(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_)); // A + B + N
158 this->ANB_->add(((*this->constantMatrix_)), ((*this->ANB_)));
159
160 // Nonlinear advection term \rho (u \cdot \nabla) u
161 // As this class is derived from NavierStokes class we can call already implemented function
162 //*************** ADVECTION************************
163 this->assemblyAdvection(elementMatrixNC);
164 elementMatrixNC->scale(this->density_);
165 this->ANB_->add((*elementMatrixNC), ((*this->ANB_)));
166
167 // For a generalized-newtonian fluid we add additional element matrix and fill it with specific contribution
168 // Remember that this term is based on the stress-divergence formulation of the momentum equation
169 // \nabla \dot \tau with \tau=\eta(\gammaDot)(\nabla u + (\nabla u)^T)
170 //*************** STRESS TENSOR************************
171 this->assemblyStress(elementMatrixN);
172 this->ANB_->add((*elementMatrixN), ((*this->ANB_)));
173
174}
175
176 // Extra stress term resulting from chosen non-newtonian constitutive model - Compute element matrix entries
177 template <class SC, class LO, class GO, class NO>
178 void AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::assemblyStress(SmallMatrixPtr_Type &elementMatrix)
179 {
180
181 int dim = this->getDim();
182 int numNodes = this->numNodesVelocity_;
183 std::string FEType = this->FETypeVelocity_;
184 int dofs = this->dofsVelocity_; // For pressure it would be 1
185
186 vec3D_dbl_ptr_Type dPhi;
187 vec_dbl_ptr_Type weights = Teuchos::rcp(new vec_dbl_Type(0));
188
189 // essential ingredients: eta(nabla(phi) ) * nabla(phi) * nabla(phi) (nonmathematical notation)
190 UN degGradPhi = Helper::determineDegree(dim, FEType, Helper::Deriv1);
191 UN extraDeg = 2; // As eta is a unknown nonlinear function of the velocity gradient we add some extra degree
192 UN deg = (degGradPhi + extraDeg) + degGradPhi + degGradPhi;
193 Helper::getDPhi(dPhi, weights, dim, FEType, deg);
194 // Example Values: dPhi->size() = 7 if number of quadrature points 7, dPhi->at(0).size() = 3 number of local element points, dPhi->at(0).at(0).size() = 2 as we have dim 2 therefore we have 2 derivatives (xi/eta in natural coordinates)
195 // Phi is defined on reference element
196
197 SC detB;
198 SC absDetB;
199 SmallMatrix<SC> B(dim);
200 SmallMatrix<SC> Binv(dim);
201
202 this->buildTransformation(B);
203 detB = B.computeInverse(Binv); // The function computeInverse returns a double value corrsponding to determinant of B
204 absDetB = std::fabs(detB); // Absolute value of B
205
206 // dPhiTrans are the transorfmed basisfunctions, so B^(-T) * \grad_phi bzw. \grad_phi^T * B^(-1)
207 // Corresponds to \hat{grad_phi}.
208 vec3D_dbl_Type dPhiTrans(dPhi->size(), vec2D_dbl_Type(dPhi->at(0).size(), vec_dbl_Type(dim, 0.)));
209 Helper::applyBTinv(dPhi, dPhiTrans, Binv); // dPhiTrans corresponds now to our basisfunction in natural coordinates
210
211 TEUCHOS_TEST_FOR_EXCEPTION(dim == 1, std::logic_error, "AssemblyStress Not implemented for dim=1");
212 //***************************************************************************
213 if (dim == 2)
214 {
215 //************************************
216 // Compute shear rate gammaDot, which is a vector because it is evaluated at each gaussian quadrature point
217 // for that first compute velocity gradient
218 vec_dbl_ptr_Type gammaDot(new vec_dbl_Type(weights->size(), 0.0)); // gammaDot->at(j) j=0...weights
219 computeShearRate(dPhiTrans, gammaDot, dim); // updates gammaDot using velocity solution
220 //************************************
221 // Compute entries
222 // Initialize some helper vectors/matrices
223 double v11, v12, v21, v22, value1_j, value2_j, value1_i, value2_i, viscosity_atw;
224 viscosity_atw = 0.;
225
226 // Construct element matrices
227 for (UN i = 0; i < numNodes; i++)
228 {
229 // Teuchos::Array<SC> value(dPhiTrans[0].size(), 0. ); // dPhiTrans[0].size() is 3
230 for (UN j = 0; j < numNodes; j++)
231 {
232 // Reset values
233 v11 = 0.0;
234 v12 = 0.0;
235 v21 = 0.0;
236 v22 = 0.0;
237
238 // So in general compute the components of eta*[ dPhiTrans_i : ( dPhiTrans_j + (dPhiTrans_j)^T )]
239 for (UN w = 0; w < dPhiTrans.size(); w++)
240 {
241
242 value1_j = dPhiTrans[w][j][0]; // so this corresponds to d\phi_j/dx
243 value2_j = dPhiTrans[w][j][1]; // so this corresponds to d\phi_j/dy
244
245 value1_i = dPhiTrans[w][i][0]; // so this corresponds to d\phi_i/dx
246 value2_i = dPhiTrans[w][i][1]; // so this corresponds to d\phi_i/dy
247
248 // viscosity function evaluated where we consider the dynamic viscosity!!
249 this->viscosityModel->evaluateMapping(this->params_, gammaDot->at(w), viscosity_atw);
250
251 v11 = v11 + viscosity_atw * weights->at(w) * (2.0 * value1_i * value1_j + value2_i * value2_j);
252 v12 = v12 + viscosity_atw * weights->at(w) * (value2_i * value1_j);
253 v21 = v21 + viscosity_atw * weights->at(w) * (value1_i * value2_j);
254 v22 = v22 + viscosity_atw * weights->at(w) * (2.0 * value2_i * value2_j + value1_i * value1_j);
255
256 } // loop end quadrature points
257
258 // multiply determinant from transformation
259 v11 *= absDetB;
260 v12 *= absDetB;
261 v21 *= absDetB;
262 v22 *= absDetB;
263
264 // Put values on the right position in element matrix - d=2 because we are in two dimensional case
265 // [v11 v12 ]
266 // [v21 v22 ]
267 (*elementMatrix)[i * dofs][j * dofs] = v11; // d=0, first dimension
268 (*elementMatrix)[i * dofs][j * dofs + 1] = v12;
269 (*elementMatrix)[i * dofs + 1][j * dofs] = v21;
270 (*elementMatrix)[i * dofs + 1][j * dofs + 1] = v22; // d=1, second dimension
271
272 } // loop end over j node
273
274 } // loop end over i node
275
276 } // end if dim 2
277 //***************************************************************************
278 else if (dim == 3)
279 {
280 //************************************#
281
282 // Compute shear rate gammaDot, which is a vector because it is evaluated at a gaussian quadrature point
283 // for that compute velocity gradient
284 vec_dbl_ptr_Type gammaDot(new vec_dbl_Type(weights->size(), 0.0)); // gammaDot->at(j) j=0...weights
285 computeShearRate(dPhiTrans, gammaDot, dim); // updates gammaDot using velcoity solution
286
287 // Initialize some helper vectors/matrices
288 double v11, v12, v13, v21, v22, v23, v31, v32, v33, value1_j, value2_j, value3_j, value1_i, value2_i, value3_i, viscosity_atw;
289
290 viscosity_atw = 0.;
291
292 // Construct element matrices
293 for (UN i = 0; i < numNodes; i++)
294 {
295 // Teuchos::Array<SC> value(dPhiTrans[0].size(), 0. ); // dPhiTrans[0].size() is 3
296
297 for (UN j = 0; j < numNodes; j++)
298 {
299 // Reset values
300 v11 = 0.0;
301 v12 = 0.0;
302 v13 = 0.0;
303 v21 = 0.0;
304 v22 = 0.0;
305 v23 = 0.0;
306 v31 = 0.0;
307 v32 = 0.0;
308 v33 = 0.0;
309
310 // So in general compute the components of eta*[ dPhiTrans_i : ( dPhiTrans_j + (dPhiTrans_j)^T )]
311 for (UN w = 0; w < dPhiTrans.size(); w++)
312 {
313
314 value1_j = dPhiTrans.at(w).at(j).at(0); // so this corresponds to d\phi_j/dx
315 value2_j = dPhiTrans.at(w).at(j).at(1); // so this corresponds to d\phi_j/dy
316 value3_j = dPhiTrans.at(w).at(j).at(2); // so this corresponds to d\phi_j/dz
317
318 value1_i = dPhiTrans.at(w).at(i).at(0); // so this corresponds to d\phi_i/dx
319 value2_i = dPhiTrans.at(w).at(i).at(1); // so this corresponds to d\phi_i/dy
320 value3_i = dPhiTrans.at(w).at(i).at(2); // so this corresponds to d\phi_i/dz
321
322 this->viscosityModel->evaluateMapping(this->params_, gammaDot->at(w), viscosity_atw);
323
324 // Construct entries - we go over all quadrature points and if j is updated we set v11 etc. again to zero
325 v11 = v11 + viscosity_atw * weights->at(w) * (2.0 * value1_j * value1_i + value2_j * value2_i + value3_j * value3_i);
326 v12 = v12 + viscosity_atw * weights->at(w) * (value2_i * value1_j);
327 v13 = v13 + viscosity_atw * weights->at(w) * (value3_i * value1_j);
328
329 v21 = v21 + viscosity_atw * weights->at(w) * (value1_i * value2_j);
330 v22 = v22 + viscosity_atw * weights->at(w) * (value1_i * value1_j + 2.0 * value2_j * value2_i + value3_j * value3_i);
331 v23 = v23 + viscosity_atw * weights->at(w) * (value3_i * value2_j);
332
333 v31 = v31 + viscosity_atw * weights->at(w) * (value1_i * value3_j);
334 v32 = v32 + viscosity_atw * weights->at(w) * (value2_i * value3_j);
335 v33 = v33 + viscosity_atw * weights->at(w) * (value1_i * value1_j + value2_i * value2_j + 2.0 * value3_i * value3_j);
336
337 } // loop end quadrature points
338
339 // multiply determinant from transformation
340 v11 *= absDetB;
341 v12 *= absDetB;
342 v13 *= absDetB;
343 v21 *= absDetB;
344 v22 *= absDetB;
345 v23 *= absDetB;
346 v31 *= absDetB;
347 v32 *= absDetB;
348 v33 *= absDetB;
349
350 // Put values on the right position in element matrix
351 // [v11 v12 v13]
352 // [v21 v22 v23]
353 // [v31 v32 v33]
354 (*elementMatrix)[i * dofs][j * dofs] = v11; // d=0, first dimension
355 (*elementMatrix)[i * dofs][j * dofs + 1] = v12;
356 (*elementMatrix)[i * dofs][j * dofs + 2] = v13;
357 (*elementMatrix)[i * dofs + 1][j * dofs] = v21;
358 (*elementMatrix)[i * dofs + 1][j * dofs + 1] = v22; // d=1, second dimension
359 (*elementMatrix)[i * dofs + 1][j * dofs + 2] = v23; // d=1, second dimension
360 (*elementMatrix)[i * dofs + 2][j * dofs] = v31;
361 (*elementMatrix)[i * dofs + 2][j * dofs + 1] = v32; // d=2, third dimension
362 (*elementMatrix)[i * dofs + 2][j * dofs + 2] = v33; // d=2, third dimension
363
364 } // loop end over j node
365 } // loop end over i node
366 } // end if dim==3
367 }
368
369 // Directional Derivative of shear stress term resulting from chosen nonlinear non-newtonian model -----
370 // Same structure and functions as in assemblyStress
371 // ( -2.0*deta/dgammaDot * dgammaDot/dTau * (0.5(dv^k + (dvh^k)^T): 0.5( dPhiTrans_j + (dPhiTrans_j)^T))0.5(dv^k + (dvh^k)^T): 0.5( dPhiTrans_i + (dPhiTrans_i)^T) )
372 template <class SC, class LO, class GO, class NO>
373 void AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::assemblyStressDev(SmallMatrixPtr_Type &elementMatrix)
374 {
375
376 int dim = this->getDim();
377 int numNodes = this->numNodesVelocity_;
378 std::string FEType = this->FETypeVelocity_;
379 int dofs = this->dofsVelocity_; // for pressure it would be 1
380
381 vec3D_dbl_ptr_Type dPhi;
382 vec_dbl_ptr_Type weights = Teuchos::rcp(new vec_dbl_Type(0));
383
384 // essential ingredients: deta/dgamma(nabla(phi)) * nabla(phi) * nabla(phi) * nabla(phi) (nonmathematical notation)
385 UN extraDegree = 2; // As deta/dgamma is a unknown nonlinear function of the velocity gradient we add some extra degree
386 UN degGradPhi = Helper::determineDegree(dim, FEType, Helper::Deriv1);
387 UN deg = 3*degGradPhi + (degGradPhi + extraDegree);
388 Helper::getDPhi(dPhi, weights, dim, FEType, deg);
389
390
391 SC detB;
392 SC absDetB;
393 SmallMatrix<SC> B(dim);
394 SmallMatrix<SC> Binv(dim);
395
396 this->buildTransformation(B);
397 detB = B.computeInverse(Binv);
398 absDetB = std::fabs(detB);
399
400 vec3D_dbl_Type dPhiTrans(dPhi->size(), vec2D_dbl_Type(dPhi->at(0).size(), vec_dbl_Type(dim, 0.)));
401 Helper::applyBTinv(dPhi, dPhiTrans, Binv);
402
403 TEUCHOS_TEST_FOR_EXCEPTION(dim == 1, std::logic_error, "AssemblyStress Not implemented for dim=1");
404
405 if (dim == 2)
406 {
407 //************************************
408 //************************************
409 // Due to the extra term related to the Gaetaeux-derivative there arise prefactors which depend on the velocity gradients solutions
410 // which therefore also have to be computed here therefore we compute it directly here
411 vec_dbl_Type u11(dPhiTrans.size(), -1.); // should correspond to du/dx at each quadrature point
412 vec_dbl_Type u12(dPhiTrans.size(), -1.); // should correspond to du/dy at each quadrature point
413 vec_dbl_Type u21(dPhiTrans.size(), -1.); // should correspond to dv/dx at each quadrature point
414 vec_dbl_Type u22(dPhiTrans.size(), -1.); // should correspond to dv/dy at each quadrature point
415 vec_dbl_ptr_Type gammaDot(new vec_dbl_Type(weights->size(), 0.0)); // gammaDot->at(j) j=0...weights
416
417 vec_dbl_ptr_Type mixed_term_xy(new vec_dbl_Type(weights->size(), 0.0));
418 for (UN w = 0; w < dPhiTrans.size(); w++)
419 { // quads points
420 // set again to zero
421 u11[w] = 0.0; // du_dx
422 u12[w] = 0.0; // du_dy
423 u21[w] = 0.0; // dv_dx
424 u22[w] = 0.0; // dv_dy
425 for (UN i = 0; i < dPhiTrans[0].size(); i++)
426 { // loop unrolling
427 LO index1 = dim * i + 0; // x
428 LO index2 = dim * i + 1; // y
429 // uLoc[d][w] += (*this->solution_)[index] * phi->at(w).at(i);
430 u11[w] += (*this->solution_)[index1] * dPhiTrans[w][i][0]; // u*dphi_dx
431 u12[w] += (*this->solution_)[index1] * dPhiTrans[w][i][1]; // because we are in 2D , 0 and 1
432 u21[w] += (*this->solution_)[index2] * dPhiTrans[w][i][0];
433 u22[w] += (*this->solution_)[index2] * dPhiTrans[w][i][1];
434 }
435 gammaDot->at(w) = std::sqrt(2.0 * u11[w] * u11[w] + 2.0 * u22[w] * u22[w] + (u12[w] + u21[w]) * (u12[w] + u21[w]));
436 mixed_term_xy->at(w) = 0.5 * (u12[w] + u21[w]);
437 }
438 //*******************************
439
440 // Initialize some helper vectors/matrices
441 double v11, v12, v21, v22, value1_j, value2_j, value1_i, value2_i, deta_dgamma_dgamma_dtau;
442
443 deta_dgamma_dgamma_dtau = 0.;
444
445 // Construct element matrices
446 for (UN i = 0; i < numNodes; i++)
447 {
448 // Teuchos::Array<SC> value(dPhiTrans[0].size(), 0. ); // dPhiTrans[0].size() is 3
449
450 for (UN j = 0; j < numNodes; j++)
451 {
452 // Reset values
453 v11 = 0.0;
454 v12 = 0.0;
455 v21 = 0.0;
456 v22 = 0.0;
457
458 // Only the part deta/dgammaDot is different for all shear thinning models (because we make the assumption of incompressibility)? NO ALSO DIFFERENT FOR EXAMPLE FOR CASSON SO CONSIDERING YIELD STRESS
459 // but we put the two terms together because then we can multiply them together and get e.g. for carreau yasuda : gammaDot^{a-2.0} which is for a=2.0 equals 0 and we do not have to worry about the problem what if gammaDot = 0.0
460 for (UN w = 0; w < dPhiTrans.size(); w++)
461 {
462
463 value1_j = dPhiTrans[w][j][0]; // so this corresponds to d\phi_j/dx
464 value2_j = dPhiTrans[w][j][1]; // so this corresponds to d\phi_j/dy
465
466 value1_i = dPhiTrans[w][i][0]; // so this corresponds to d\phi_i/dx
467 value2_i = dPhiTrans[w][i][1]; // so this corresponds to d\phi_i/dy
468
469 this->viscosityModel->evaluateDerivative(this->params_, gammaDot->at(w), deta_dgamma_dgamma_dtau);
470 /* EInfacher in unausmultiplizierter Form
471 v11 = v11 + (-2.0)*deta_dgamma_dgamma_dtau * weights->at(w) *(u11[w]*u11[w]*value1_i*value1_j+u11[w]*mixed_terms->at(w)*(value1_i*value2_j+value2_i*value1_j)+ mixed_terms->at(w)*mixed_terms->at(w)*(value2_i*value2_j)); // xx contribution: (dv_x/dx)^2*dphi_i/dx*dphi_j/dx+dv_x/dx*f*(dphi_i/dx*dphi_j/dy+dphi_i/dy*dphi_j/dx)+f^2*dphi_i/dy*dphi_j/dy
472 v12 = v12 + (-2.0)*deta_dgamma_dgamma_dtau * weights->at(w) *(u11[w]*mixed_terms->at(w)*value1_i*value1_j+u11[w]*u22[w]*(value1_i*value2_j) +mixed_terms->at(w)*mixed_terms->at(w)*value2_i*value1_j+mixed_terms->at(w)*u22[w]*value2_i*value2_j ); // xy contribution: dv_x/dx*f*dphi_i/dx*dphi_j/dx+dv_x/dx*dv_y/dy*dphi_i/dx*dphi_j/dy+f^2*dphi_i_dy*dphi_j/dx+f*dv_y/dy*dphi_i/dy*dphi_j/dy
473 v21 = v21 + (-2.0)*deta_dgamma_dgamma_dtau * weights->at(w) *(u11[w]*mixed_terms->at(w)*value1_i*value1_j+mixed_terms->at(w)*mixed_terms->at(w)*value1_i*value2_j+u11[w]*u22[w]*value2_i*value1_j +mixed_terms->at(w)*u22[w]*value2_i*value2_j ); // yx contribution: dv_x/dx*f*dphi_i/dx*dphi_j/dx+dv_x/dx*dv_y/dy*dphi_i/dy*dphi_j/dx+f^2*dphi_i_dx*dphi_j/dy+f*dv_y/dy*dphi_i/dy*dphi_j/dy
474 v22 = v22 + (-2.0)*deta_dgamma_dgamma_dtau * weights->at(w) *(u22[w]*u22[w]*value2_i*value2_j+u22[w]*mixed_terms->at(w)*(value1_i*value2_j+value2_i*value1_j)+ mixed_terms->at(w)*mixed_terms->at(w)*(value1_i*value1_j) ); // yy contribution: (dv_y/dy)^2*dphi_i/dy*dphi_j/dy+dv_y/dy*f*(dphi_i/dx*dphi_j/dy+dphi_i/dy*dphi_j/dx)+f^2*dphi_i/dx*dphi_j/dx
475 */
476 v11 = v11 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_j * u11[w] + value2_j * mixed_term_xy->at(w)) * (value1_i * u11[w] + value2_i * mixed_term_xy->at(w)));
477 v12 = v12 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_j * mixed_term_xy->at(w) + u22[w] * value2_j) * (value1_i * u11[w] + value2_i * mixed_term_xy->at(w)));
478
479 v21 = v21 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_j * u11[w] + value2_j * mixed_term_xy->at(w)) * (value1_i * mixed_term_xy->at(w) + value2_i * u22[w]));
480 v22 = v22 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_j * mixed_term_xy->at(w) + u22[w] * value2_j) * (value1_i * mixed_term_xy->at(w) + value2_i * u22[w]));
481
482 } // loop end quadrature points
483
484 // multiply determinant from transformation
485 v11 *= absDetB;
486 v12 *= absDetB;
487 v21 *= absDetB;
488 v22 *= absDetB;
489
490 // Put values on the right position in element matrix - d=2 because we are in two dimensional case
491 // [v11 v12 ]
492 // [v21 v22 ]
493 (*elementMatrix)[i * dofs][j * dofs] = v11; // d=0, first dimension
494 (*elementMatrix)[i * dofs][j * dofs + 1] = v12;
495 (*elementMatrix)[i * dofs + 1][j * dofs] = v21;
496 (*elementMatrix)[i * dofs + 1][j * dofs + 1] = v22; // d=1, second dimension
497
498 } // loop end over j node
499
500 } // loop end over i node
501
502 } // end if dim 2
503 //***************************************************************************
504 //***************************************************************************
505 else if (dim == 3)
506 {
507 //************************************
508 // Compute shear rate gammaDot, which is a vector because it is evaluated at a gaussian quadrature point
509 // for that compute velocity gradient
510 vec_dbl_Type u11(dPhiTrans.size(), -1.); // should correspond to du/dx at each quadrature point
511 vec_dbl_Type u12(dPhiTrans.size(), -1.); // should correspond to du/dy at each quadrature point
512 vec_dbl_Type u13(dPhiTrans.size(), -1.); // should correspond to du/dz at each quadrature point
513
514 vec_dbl_Type u21(dPhiTrans.size(), -1.); // should correspond to dv/dx at each quadrature point
515 vec_dbl_Type u22(dPhiTrans.size(), -1.); // should correspond to dv/dy at each quadrature point
516 vec_dbl_Type u23(dPhiTrans.size(), -1.); // should correspond to dv/dz at each quadrature point
517
518 vec_dbl_Type u31(dPhiTrans.size(), -1.); // should correspond to dw/dx at each quadrature point
519 vec_dbl_Type u32(dPhiTrans.size(), -1.); // should correspond to dw/dy at each quadrature point
520 vec_dbl_Type u33(dPhiTrans.size(), -1.); // should correspond to dw/dz at each quadrature point
521
522 vec_dbl_ptr_Type gammaDot(new vec_dbl_Type(weights->size(), 0.0)); // gammaDot->at(j) j=0...weights
523
524 vec_dbl_ptr_Type mixed_term_xy(new vec_dbl_Type(weights->size(), 0.0));
525 vec_dbl_ptr_Type mixed_term_xz(new vec_dbl_Type(weights->size(), 0.0));
526 vec_dbl_ptr_Type mixed_term_yz(new vec_dbl_Type(weights->size(), 0.0));
527
528 for (UN w = 0; w < dPhiTrans.size(); w++)
529 { // quads points
530
531 u11[w] = 0.0;
532 u12[w] = 0.0;
533 u13[w] = 0.0;
534 u21[w] = 0.0;
535 u22[w] = 0.0;
536 u23[w] = 0.0;
537 u31[w] = 0.0;
538 u32[w] = 0.0;
539 u33[w] = 0.0;
540
541 for (UN i = 0; i < dPhiTrans[0].size(); i++)
542 {
543 LO index1 = dim * i + 0; // x
544 LO index2 = dim * i + 1; // y
545 LO index3 = dim * i + 2; // z
546 // uLoc[d][w] += (*this->solution_)[index] * phi->at(w).at(i);
547 u11[w] += (*this->solution_)[index1] * dPhiTrans[w][i][0]; // u*dphi_dx
548 u12[w] += (*this->solution_)[index1] * dPhiTrans[w][i][1]; // because we are in 3D , 0 and 1, 2
549 u13[w] += (*this->solution_)[index1] * dPhiTrans[w][i][2];
550 u21[w] += (*this->solution_)[index2] * dPhiTrans[w][i][0]; // v*dphi_dx
551 u22[w] += (*this->solution_)[index2] * dPhiTrans[w][i][1];
552 u23[w] += (*this->solution_)[index2] * dPhiTrans[w][i][2];
553 u31[w] += (*this->solution_)[index3] * dPhiTrans[w][i][0]; // w*dphi_dx
554 u32[w] += (*this->solution_)[index3] * dPhiTrans[w][i][1];
555 u33[w] += (*this->solution_)[index3] * dPhiTrans[w][i][2];
556 }
557 gammaDot->at(w) = std::sqrt(2.0 * u11[w] * u11[w] + 2.0 * u22[w] * u22[w] + 2.0 * u33[w] * u33[w] + (u12[w] + u21[w]) * (u12[w] + u21[w]) + (u13[w] + u31[w]) * (u13[w] + u31[w]) + (u23[w] + u32[w]) * (u23[w] + u32[w]));
558
559 mixed_term_xy->at(w) = 0.5 * (u12[w] + u21[w]);
560 mixed_term_xz->at(w) = 0.5 * (u31[w] + u13[w]);
561 mixed_term_yz->at(w) = 0.5 * (u32[w] + u23[w]);
562 }
563
564 // Initialize some helper vectors/matrices
565 double v11, v12, v13, v21, v22, v23, v31, v32, v33, value1_j, value2_j, value3_j, value1_i, value2_i, value3_i, deta_dgamma_dgamma_dtau;
566
567 deta_dgamma_dgamma_dtau = 0.;
568
569 // Construct element matrices
570 for (UN i = 0; i < numNodes; i++)
571 {
572 // Teuchos::Array<SC> value(dPhiTrans[0].size(), 0. ); // dPhiTrans[0].size() is 3
573
574 for (UN j = 0; j < numNodes; j++)
575 {
576 // Reset values
577 v11 = 0.0;
578 v12 = 0.0;
579 v13 = 0.0;
580 v21 = 0.0;
581 v22 = 0.0;
582 v23 = 0.0;
583 v31 = 0.0;
584 v32 = 0.0;
585 v33 = 0.0;
586
587 // So in general compute the components of eta*[ dPhiTrans_i : ( dPhiTrans_j + (dPhiTrans_j)^T )]
588 for (UN w = 0; w < dPhiTrans.size(); w++)
589 {
590
591 value1_j = dPhiTrans.at(w).at(j).at(0); // so this corresponds to d\phi_j/dx
592 value2_j = dPhiTrans.at(w).at(j).at(1); // so this corresponds to d\phi_j/dy
593 value3_j = dPhiTrans.at(w).at(j).at(2); // so this corresponds to d\phi_j/dz
594
595 value1_i = dPhiTrans.at(w).at(i).at(0); // so this corresponds to d\phi_i/dx
596 value2_i = dPhiTrans.at(w).at(i).at(1); // so this corresponds to d\phi_i/dy
597 value3_i = dPhiTrans.at(w).at(i).at(2); // so this corresponds to d\phi_i/dz
598
599 this->viscosityModel->evaluateDerivative(this->params_, gammaDot->at(w), deta_dgamma_dgamma_dtau);
600
601 // Construct entries - we go over all quadrature points and if j is updated we set v11 etc. again to zero
602 v11 = v11 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_j * u11[w] + value2_j * mixed_term_xy->at(w) + value3_j * mixed_term_xz->at(w)) * (value1_i * u11[w] + value2_i * mixed_term_xy->at(w) + value3_i * mixed_term_xz->at(w)));
603 v12 = v12 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_j * mixed_term_xy->at(w) + u22[w] * value2_j + mixed_term_yz->at(w) * value3_j) * (value1_i * u11[w] + value2_i * mixed_term_xy->at(w) + value3_i * mixed_term_xz->at(w)));
604 v13 = v13 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value3_j * mixed_term_xz->at(w) + value2_j * mixed_term_yz->at(w) + value3_j * u33[w]) * (value1_i * u11[w] + value2_i * mixed_term_xy->at(w) + value3_i * mixed_term_xz->at(w)));
605
606 v21 = v21 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_j * u11[w] + value2_j * mixed_term_xy->at(w) + value3_j * mixed_term_xz->at(w)) * (value1_i * mixed_term_xy->at(w) + value2_i * u22[w] + value3_i * mixed_term_yz->at(w)));
607 v22 = v22 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_j * mixed_term_xy->at(w) + u22[w] * value2_j + mixed_term_yz->at(w) * value3_j) * (value1_i * mixed_term_xy->at(w) + value2_i * u22[w] + value3_i * mixed_term_yz->at(w)));
608 v23 = v23 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value3_j * mixed_term_xz->at(w) + value2_j * mixed_term_yz->at(w) + value3_j * u33[w]) * (value1_i * mixed_term_xy->at(w) + value2_i * u22[w] + value3_i * mixed_term_yz->at(w)));
609
610 v31 = v31 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_j * u11[w] + value2_j * mixed_term_xy->at(w) + value3_j * mixed_term_xz->at(w)) * (value1_i * mixed_term_xz->at(w) + mixed_term_yz->at(w) * value2_i + u33[w] * value3_i));
611 v32 = v32 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_j * mixed_term_xy->at(w) + u22[w] * value2_j + mixed_term_yz->at(w) * value3_j) * (value1_i * mixed_term_xz->at(w) + mixed_term_yz->at(w) * value2_i + u33[w] * value3_i));
612 v33 = v33 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value3_j * mixed_term_xz->at(w) + value2_j * mixed_term_yz->at(w) + value3_j * u33[w]) * (value1_i * mixed_term_xz->at(w) + mixed_term_yz->at(w) * value2_i + u33[w] * value3_i));
613
614 } // loop end quadrature points
615
616 // multiply determinant from transformation
617 v11 *= absDetB;
618 v12 *= absDetB;
619 v13 *= absDetB;
620 v21 *= absDetB;
621 v22 *= absDetB;
622 v23 *= absDetB;
623 v31 *= absDetB;
624 v32 *= absDetB;
625 v33 *= absDetB;
626
627 // Put values on the right position in element matrix
628 // [v11 v12 v13]
629 // [v21 v22 v23]
630 // [v31 v32 v33]
631 (*elementMatrix)[i * dofs][j * dofs] = v11; // d=0, first dimension
632 (*elementMatrix)[i * dofs][j * dofs + 1] = v12;
633 (*elementMatrix)[i * dofs][j * dofs + 2] = v13;
634 (*elementMatrix)[i * dofs + 1][j * dofs] = v21;
635 (*elementMatrix)[i * dofs + 1][j * dofs + 1] = v22; // d=1, second dimension
636 (*elementMatrix)[i * dofs + 1][j * dofs + 2] = v23; // d=1, second dimension
637 (*elementMatrix)[i * dofs + 2][j * dofs] = v31;
638 (*elementMatrix)[i * dofs + 2][j * dofs + 1] = v32; // d=2, third dimension
639 (*elementMatrix)[i * dofs + 2][j * dofs + 2] = v33; // d=2, third dimension
640
641 } // loop end over j node
642
643 } // loop end over i node
644
645 } // end if dim = 3
646 }
647
648
649 //ToDo: Using new functionalty integrate again assembly of Neumann boundary term
650
651 // "Fixpunkt"- Matrix without jacobian for calculating Ax
652 // Here update please to unlinearized System Matrix accordingly.
653 template <class SC, class LO, class GO, class NO>
654 void AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::assembleRHS()
655 {
656
657 SmallMatrixPtr_Type elementMatrixN = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
658 SmallMatrixPtr_Type elementMatrixNC = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
659
660 this->ANB_.reset(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_)); // A + B + N
661 this->ANB_->add((*this->constantMatrix_), (*this->ANB_));
662
663 // Nonlinear shear stress tensor *******************************
664 this->assemblyStress(elementMatrixN);
665 this->ANB_->add((*elementMatrixN), (*this->ANB_));
666
667 // Nonlinear convection term *******************************
668 this->assemblyAdvection(elementMatrixNC);
669 elementMatrixNC->scale(this->density_);
670 this->ANB_->add((*elementMatrixNC), (*this->ANB_));
671
672 // @ToDo If underlying element is an outflow boundary element - will be added when func nonlinear boundar term *******************************
673 /*if (this->surfaceElement == true)
674 {
675 SmallMatrixPtr_Type elementMatrixNB = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
676 this->assemblyNeumannBoundaryTerm(elementMatrixNB);
677 this->ANB_->add((*elementMatrixNB), ((*this->ANB_)));
678 }
679 */
680
681 this->rhsVec_.reset(new vec_dbl_Type(this->dofsElement_, 0.));
682 // Multiplying ANB_ * solution // ANB Matrix without nonlinear part.
683 int s = 0, t = 0;
684 for (int i = 0; i < this->ANB_->size(); i++)
685 {
686 if (i >= this->dofsElementVelocity_)
687 s = 1;
688 for (int j = 0; j < this->ANB_->size(); j++)
689 {
690 if (j >= this->dofsElementVelocity_)
691 t = 1;
692 (*this->rhsVec_)[i] += (*this->ANB_)[i][j] * (*this->solution_)[j] * this->coeff_[s][t];
693 }
694 t = 0;
695 }
696 }
697
701 /* In 2D B=[ nodesRefConfig(1).at(0)-nodesRefConfig(0).at(0) nodesRefConfig(2).at(0)-nodesRefConfig(0).at(0) ]
702 [ nodesRefConfig(1).at(1)-nodesRefConfig(0).at(1) nodesRefConfig(2).at(1)-nodesRefConfig(0).at(1) ]
717
718
719 // Compute Shear Rate on quadrature points depending on gradient of velocity solution at nodes
720 template <class SC, class LO, class GO, class NO>
721 void AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::computeShearRate(vec3D_dbl_Type dPhiTrans,
722 vec_dbl_ptr_Type &gammaDot, int dim)
723 {
724
725 //****************** TWO DIMENSIONAL *********************************
726 if (dim == 2)
727 {
728
729 vec_dbl_Type u11(dPhiTrans.size(), -1.); // should correspond to du/dx at each quadrature point
730 vec_dbl_Type u12(dPhiTrans.size(), -1.); // should correspond to du/dy at each quadrature point
731 vec_dbl_Type u21(dPhiTrans.size(), -1.); // should correspond to dv/dx at each quadrature point
732 vec_dbl_Type u22(dPhiTrans.size(), -1.); // should correspond to dv/dy at each quadrature point
733
734 for (UN w = 0; w < dPhiTrans.size(); w++)
735 { // quads points
736 // set again to zero
737 u11[w] = 0.0;
738 u12[w] = 0.0;
739 u21[w] = 0.0;
740 u22[w] = 0.0;
741 for (UN i = 0; i < dPhiTrans[0].size(); i++)
742 { // loop unrolling
743 LO index1 = dim * i + 0; // x
744 LO index2 = dim * i + 1; // y
745 // uLoc[d][w] += (*this->solution_)[index] * phi->at(w).at(i);
746 u11[w] += (*this->solution_)[index1] * dPhiTrans[w][i][0]; // u*dphi_dx
747 u12[w] += (*this->solution_)[index1] * dPhiTrans[w][i][1]; // because we are in 2D , 0 and 1
748 u21[w] += (*this->solution_)[index2] * dPhiTrans[w][i][0];
749 u22[w] += (*this->solution_)[index2] * dPhiTrans[w][i][1];
750 }
751 gammaDot->at(w) = std::sqrt(2.0 * u11[w] * u11[w] + 2.0 * u22[w] * u22[w] + (u12[w] + u21[w]) * (u12[w] + u21[w]));
752 }
753 } // end if dim == 2
754 //****************** THREE DIMENSIONAL *********************************
755 else if (dim == 3)
756 {
757
758 vec_dbl_Type u11(dPhiTrans.size(), -1.); // should correspond to du/dx at each quadrature point
759 vec_dbl_Type u12(dPhiTrans.size(), -1.); // should correspond to du/dy at each quadrature point
760 vec_dbl_Type u13(dPhiTrans.size(), -1.); // should correspond to du/dz at each quadrature point
761
762 vec_dbl_Type u21(dPhiTrans.size(), -1.); // should correspond to dv/dx at each quadrature point
763 vec_dbl_Type u22(dPhiTrans.size(), -1.); // should correspond to dv/dy at each quadrature point
764 vec_dbl_Type u23(dPhiTrans.size(), -1.); // should correspond to dv/dz at each quadrature point
765
766 vec_dbl_Type u31(dPhiTrans.size(), -1.); // should correspond to dw/dx at each quadrature point
767 vec_dbl_Type u32(dPhiTrans.size(), -1.); // should correspond to dw/dy at each quadrature point
768 vec_dbl_Type u33(dPhiTrans.size(), -1.); // should correspond to dw/dz at each quadrature point
769
770 for (UN w = 0; w < dPhiTrans.size(); w++)
771 { // quads points
772 // set again to zero
773 u11[w] = 0.0;
774 u12[w] = 0.0;
775 u13[w] = 0.0;
776 u21[w] = 0.0;
777 u22[w] = 0.0;
778 u23[w] = 0.0;
779 u31[w] = 0.0;
780 u32[w] = 0.0;
781 u33[w] = 0.0;
782
783 for (UN i = 0; i < dPhiTrans[0].size(); i++)
784 {
785 LO index1 = dim * i + 0; // x
786 LO index2 = dim * i + 1; // y
787 LO index3 = dim * i + 2; // z
788 // uLoc[d][w] += (*this->solution_)[index] * phi->at(w).at(i);
789 u11[w] += (*this->solution_)[index1] * dPhiTrans[w][i][0]; // u*dphi_dx
790 u12[w] += (*this->solution_)[index1] * dPhiTrans[w][i][1]; // because we are in 3D , 0 and 1, 2
791 u13[w] += (*this->solution_)[index1] * dPhiTrans[w][i][2];
792 u21[w] += (*this->solution_)[index2] * dPhiTrans[w][i][0]; // v*dphi_dx
793 u22[w] += (*this->solution_)[index2] * dPhiTrans[w][i][1];
794 u23[w] += (*this->solution_)[index2] * dPhiTrans[w][i][2];
795 u31[w] += (*this->solution_)[index3] * dPhiTrans[w][i][0]; // w*dphi_dx
796 u32[w] += (*this->solution_)[index3] * dPhiTrans[w][i][1];
797 u33[w] += (*this->solution_)[index3] * dPhiTrans[w][i][2];
798 }
799 gammaDot->at(w) = std::sqrt(2.0 * u11[w] * u11[w] + 2.0 * u22[w] * u22[w] + 2.0 * u33[w] * u33[w] + (u12[w] + u21[w]) * (u12[w] + u21[w]) + (u13[w] + u31[w]) * (u13[w] + u31[w]) + (u23[w] + u32[w]) * (u23[w] + u32[w]));
800 }
801 } // end if dim == 3
802 }
803
804 /* Based on the current solution (velocity, pressure etc.) we want to be able to compute postprocessing fields
805 like here the viscosity inside an element.
806 */
807 template <class SC, class LO, class GO, class NO>
808 void AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::computeLocalconstOutputField()
809 {
810 int dim = this->getDim();
811 std::string FEType = this->FETypeVelocity_;
812
813 SC detB;
814 SmallMatrix<SC> B(dim);
815 SmallMatrix<SC> Binv(dim);
816
817 this->buildTransformation(B);
818 detB = B.computeInverse(Binv);
819
820 vec3D_dbl_ptr_Type dPhiAtCM;
821
822 // Compute viscosity at center of mass using nodal values and shape function **********************************************************************************
823 TEUCHOS_TEST_FOR_EXCEPTION(dim == 1, std::logic_error, "computeLocalconstOutputField Not implemented for dim=1");
824
825 Helper::getDPhiAtCM(dPhiAtCM, dim, FEType); // These are the original coordinates of the reference element
826 vec3D_dbl_Type dPhiTransAtCM(dPhiAtCM->size(), vec2D_dbl_Type(dPhiAtCM->at(0).size(), vec_dbl_Type(dim, 0.)));
827 Helper::applyBTinv(dPhiAtCM, dPhiTransAtCM, Binv); // We need transformation because of velocity gradient in shear rate equation
828
829 vec_dbl_ptr_Type gammaDoti(new vec_dbl_Type(dPhiAtCM->size(), 0.0)); // Only one value because size is one
830 computeShearRate(dPhiTransAtCM, gammaDoti, dim); // updates gammaDot using velocity solution
831 this->viscosityModel->evaluateMapping(this->params_, gammaDoti->at(0), this->constOutputField_.at(0));
832 }
833
834}
835#endif
FEDD::AssembleFEGeneralizedNewtonian::assembleRHS
void assembleRHS() override
Assemble the element right hand side vector.
Definition AssembleFEGeneralizedNewtonian_def.hpp:654
FEDD::AssembleFEGeneralizedNewtonian::assembleJacobian
void assembleJacobian() override
Assemble the element Jacobian matrix.
Definition AssembleFEGeneralizedNewtonian_def.hpp:64
FEDD::AssembleFEGeneralizedNewtonian::AssembleFEGeneralizedNewtonian
AssembleFEGeneralizedNewtonian(int flag, vec2D_dbl_Type nodesRefConfig, ParameterListPtr_Type parameters, tuple_disk_vec_ptr_Type tuple)
Constructor for AssembleFEAceNavierStokes.
Definition AssembleFEGeneralizedNewtonian_def.hpp:33
FEDD::AssembleFEGeneralizedNewtonian< SC, LO, GO, NO >::assemblyStress
void assemblyStress(SmallMatrixPtr_Type &elementMatrix)
FEDD::AssembleFEGeneralizedNewtonian::computeLocalconstOutputField
void computeLocalconstOutputField() override
Compute the viscosity for an element depending on the knwon velocity solution.
Definition AssembleFEGeneralizedNewtonian_def.hpp:808
FEDD::AssembleFEGeneralizedNewtonian::computeShearRate
void computeShearRate(vec3D_dbl_Type dPhiTrans, vec_dbl_ptr_Type &gammaDot, int dim)
Computation of shear rate using the current velocity solution at the nodes and the derivative of the ...
Definition AssembleFEGeneralizedNewtonian_def.hpp:721
FEDD::AssembleFEGeneralizedNewtonian< SC, LO, GO, NO >::assemblyStressDev
void assemblyStressDev(SmallMatrixPtr_Type &elementMatrix)
FEDD::AssembleFENavierStokes< default_sc, default_lo, default_go, default_no >::dofsVelocity_
int dofsVelocity_
Definition AssembleFENavierStokes_decl.hpp:110
FEDD::AssembleFENavierStokes< default_sc, default_lo, default_go, default_no >::buildTransformation
void buildTransformation(SmallMatrix< default_sc > &B)
Definition AssembleFENavierStokes_def.hpp:478
FEDD::AssembleFENavierStokes< SC, LO, GO, NO >::assemblyDivAndDivT
void assemblyDivAndDivT(SmallMatrixPtr_Type &elementMatrix)
FEDD::AssembleFENavierStokes< SC, LO, GO, NO >::assemblyAdvectionInU
void assemblyAdvectionInU(SmallMatrixPtr_Type &elementMatrix)
FEDD::AssembleFENavierStokes< SC, LO, GO, NO >::assemblyAdvection
void assemblyAdvection(SmallMatrixPtr_Type &elementMatrix)
FEDD::AssembleFENavierStokes< default_sc, default_lo, default_go, default_no >::AssembleFENavierStokes
AssembleFENavierStokes(int flag, vec2D_dbl_Type nodesRefConfig, ParameterListPtr_Type parameters, tuple_disk_vec_ptr_Type tuple)
Definition AssembleFENavierStokes_def.hpp:9
FEDD::CarreauYasuda
Definition CarreauYasuda_decl.hpp:67
FEDD::Dimless_Carreau
This class is derived from the abstract class DifferentiableFuncClass and should provide functionalit...
Definition Dimless_Carreau_decl.hpp:68
FEDD::Helper::determineDegree
static UN determineDegree(UN dim, std::string FEType, VarType orderOfDerivative)
Determine polynomial degree of a finite element basis function or its gradient that is required to se...
Definition Helper.cpp:68
FEDD::Helper::applyBTinv
static void applyBTinv(vec3D_dbl_ptr_Type &dPhiIn, vec3D_dbl_Type &dPhiOut, const SmallMatrix< SC > &Binv)
Applying the transformation matriX B to the gradient of phi, as is done in when transforming the grad...
Definition Helper.cpp:200
FEDD::Helper::Deriv1
@ Deriv1
order 1, gradient(f(x))
Definition Helper.hpp:28
FEDD::Helper::getDPhiAtCM
static void getDPhiAtCM(vec3D_dbl_ptr_Type &DPhi, int dim, std::string FEType)
Definition Helper.cpp:1470
FEDD::Helper::getDPhi
static int getDPhi(vec3D_dbl_ptr_Type &DPhi, vec_dbl_ptr_Type &weightsDPhi, int Dimension, std::string FEType, int Degree)
Full matrix representation of gradient of a basis function for each quadrature point.
Definition Helper.cpp:215
FEDD::PowerLaw
This class is derived from the abstract class DifferentiableFuncClass and should provide functionalit...
Definition PowerLaw_decl.hpp:69
FEDD::SmallMatrix
This class represents a templated small Matrix of type T. Primarily created for 2x2 and 3x3 matrices....
Definition SmallMatrix.hpp:22
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement.cpp:5