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