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 nabla \cdot 2 eta(gamma) D(v) -------------------------------------------------------------
368 For literature, see for example, [On preconditioning of incompressible non-newtonian flow problems, He et al., 2015]
369 D(v) = 0.5( \nabla v + (\nabla v)^T ) is the rate of deformation tensor, which is a symmetric tensor
370 The shear rate is defined by gamma_dot = 2 sqrt ( - PI_D )
371 with PI_D being the second invariant of the rate of deformation tensor defined for incompressible flows as: PI_D = - tr( D^2 )/2 < 0
372 In the derivation of the Gateueux derivative, the expression d/depsilon ( eta ( D(v+ epsioln delta v ) ) |_epsilon = 0 occurs
373 from which with the chain rule we obtain deta/dgammaDot * dgammaDot/dPI_D * dPI_D/dD(v) dD(v)/d epsilon | epsilon = 0
374 Looking at the single terms and setting epsilon = 0 we obtain:
375 deta/dgammaDot: --> derivative of the viscosity model w.r.t. shear rate - Depending on the model
376 dgammaDot/dTau: --> derivative of shear rate w.r.t. second invariant of rate of deformation tensor
377 = d/dPI_D ( 2 sqrt( - PI_D ) ) = -2.0 / (2 * sqrt(-PI_D) ) = -1.0 / ( sqrt(-PI_D) ) = -2.0 / gammaDot
378 dPI_D/dD(v) : --> derivative of second invariant of rate of deformation tensor w.r.t. rate of deformation tensor D(v)
379 = d/dD(v) ( - tr( D^2 )/2 ) = - D(v)
380 dD(v)/d epsilon | epsilon = 0 : --> derivative of rate of deformation tensor w.r.t. epsilon at epsilon = 0
381 = D( delta v )
382 Putting all together we obtain the Gateueux derivative of the shear stress term as:
383 nabla \cdot ( -2.0 *deta/dgammaDot * dgammaDot/dTau * ( D(v^k) : D( delta v) ) * ( D(v^k) ) )
384 and in the weak formulation we have to test with basis functions phi_j leading to the element matrix entries:
385 ( -2.0 *deta/dgammaDot * dgammaDot/dTau * ( D(v^k) : D(\phi_i) ) * ( D(v^k) : D(\phi_j) ) )
386 -----> This tensor is symmetric
387 Importantly we will summarize deta/dgammaDot * dgammaDot/dTau into a single term, as it might be beneficial to cancel gammaDot from the denominator
388 ------------------------------------------------------------------------------------------------------------------------------------------------------------------
389 */
390 template <class SC, class LO, class GO, class NO>
391 void AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::assemblyStressDev(SmallMatrixPtr_Type &elementMatrix) // Same basic structure/ functions as in assemblyStress
392
393 {
394
395 int dim = this->getDim();
396 int numNodes = this->numNodesVelocity_;
397 std::string FEType = this->FETypeVelocity_;
398 int dofs = this->dofsVelocity_; // for pressure it would be 1
399
400 vec3D_dbl_ptr_Type dPhi;
401 vec_dbl_ptr_Type weights = Teuchos::rcp(new vec_dbl_Type(0));
402
403 // essential ingredients: deta/dgamma(nabla(phi)) * nabla(phi) * nabla(phi) * nabla(phi) (nonmathematical notation)
404 UN extraDegree = 2; // As deta/dgamma is a unknown nonlinear function of the velocity gradient we add some extra degree
405 UN degGradPhi = Helper::determineDegree(dim, FEType, Helper::Deriv1);
406 UN deg = 3*degGradPhi + (degGradPhi + extraDegree);
407 Helper::getDPhi(dPhi, weights, dim, FEType, deg);
408
409
410 SC detB;
411 SC absDetB;
412 SmallMatrix<SC> B(dim);
413 SmallMatrix<SC> Binv(dim);
414
415 this->buildTransformation(B);
416 detB = B.computeInverse(Binv);
417 absDetB = std::fabs(detB);
418
419 vec3D_dbl_Type dPhiTrans(dPhi->size(), vec2D_dbl_Type(dPhi->at(0).size(), vec_dbl_Type(dim, 0.)));
420 Helper::applyBTinv(dPhi, dPhiTrans, Binv);
421
422 TEUCHOS_TEST_FOR_EXCEPTION(dim == 1, std::logic_error, "AssemblyStress Not implemented for dim=1");
423
424 if (dim == 2)
425 {
426 //************************************
427 //************************************
428 // Due to the extra term related to the Gaetaeux-derivative there arise prefactors which depend on the velocity gradients solutions
429 // which therefore also have to be computed here therefore we compute it directly here
430 vec_dbl_Type u11(dPhiTrans.size(), -1.); // should correspond to du/dx at each quadrature point
431 vec_dbl_Type u12(dPhiTrans.size(), -1.); // should correspond to du/dy at each quadrature point
432 vec_dbl_Type u21(dPhiTrans.size(), -1.); // should correspond to dv/dx at each quadrature point
433 vec_dbl_Type u22(dPhiTrans.size(), -1.); // should correspond to dv/dy at each quadrature point
434 vec_dbl_ptr_Type gammaDot(new vec_dbl_Type(weights->size(), 0.0)); // gammaDot->at(j) j=0...weights
435
436 vec_dbl_ptr_Type mixed_term_xy(new vec_dbl_Type(weights->size(), 0.0));
437 for (UN w = 0; w < dPhiTrans.size(); w++)
438 { // quads points
439 // set again to zero
440 u11[w] = 0.0; // du_dx
441 u12[w] = 0.0; // du_dy
442 u21[w] = 0.0; // dv_dx
443 u22[w] = 0.0; // dv_dy
444 for (UN i = 0; i < dPhiTrans[0].size(); i++)
445 { // loop unrolling
446 LO index1 = dim * i + 0; // x
447 LO index2 = dim * i + 1; // y
448 // uLoc[d][w] += (*this->solution_)[index] * phi->at(w).at(i);
449 u11[w] += (*this->solution_)[index1] * dPhiTrans[w][i][0]; // u*dphi_dx
450 u12[w] += (*this->solution_)[index1] * dPhiTrans[w][i][1]; // because we are in 2D , 0 and 1
451 u21[w] += (*this->solution_)[index2] * dPhiTrans[w][i][0];
452 u22[w] += (*this->solution_)[index2] * dPhiTrans[w][i][1];
453 }
454 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]));
455 mixed_term_xy->at(w) = 0.5 * (u12[w] + u21[w]);
456 }
457 //*******************************
458
459 // Initialize some helper vectors/matrices
460 double v11, v12, v21, v22, value1_j, value2_j, value1_i, value2_i, deta_dgamma_dgamma_dtau;
461
462 deta_dgamma_dgamma_dtau = 0.;
463
464 // Construct element matrices
465 for (UN i = 0; i < numNodes; i++)
466 {
467 // Teuchos::Array<SC> value(dPhiTrans[0].size(), 0. ); // dPhiTrans[0].size() is 3
468
469 for (UN j = 0; j < numNodes; j++)
470 {
471 // Reset values
472 v11 = 0.0;
473 v12 = 0.0;
474 v21 = 0.0;
475 v22 = 0.0;
476
477 // 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
478 // 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
479 for (UN w = 0; w < dPhiTrans.size(); w++)
480 {
481
482 value1_j = dPhiTrans[w][j][0]; // so this corresponds to d\phi_j/dx
483 value2_j = dPhiTrans[w][j][1]; // so this corresponds to d\phi_j/dy
484
485 value1_i = dPhiTrans[w][i][0]; // so this corresponds to d\phi_i/dx
486 value2_i = dPhiTrans[w][i][1]; // so this corresponds to d\phi_i/dy
487
488 this->viscosityModel->evaluateDerivative(this->params_, gammaDot->at(w), deta_dgamma_dgamma_dtau);
489 /* EInfacher in unausmultiplizierter Form
490 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
491 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
492 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
493 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
494 */
495 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)));
496 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)));
497
498 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]));
499 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]));
500
501 } // loop end quadrature points
502
503 // multiply determinant from transformation
504 v11 *= absDetB;
505 v12 *= absDetB;
506 v21 *= absDetB;
507 v22 *= absDetB;
508
509 // Put values on the right position in element matrix - d=2 because we are in two dimensional case
510 // [v11 v12 ]
511 // [v21 v22 ]
512 (*elementMatrix)[i * dofs][j * dofs] = v11; // d=0, first dimension
513 (*elementMatrix)[i * dofs][j * dofs + 1] = v12;
514 (*elementMatrix)[i * dofs + 1][j * dofs] = v21;
515 (*elementMatrix)[i * dofs + 1][j * dofs + 1] = v22; // d=1, second dimension
516
517 } // loop end over j node
518
519 } // loop end over i node
520
521 } // end if dim 2
522 //***************************************************************************
523 //***************************************************************************
524 else if (dim == 3)
525 {
526 //************************************
527 // Compute shear rate gammaDot, which is a vector because it is evaluated at a gaussian quadrature point
528 // for that compute velocity gradient
529 vec_dbl_Type u11(dPhiTrans.size(), -1.); // should correspond to du/dx at each quadrature point
530 vec_dbl_Type u12(dPhiTrans.size(), -1.); // should correspond to du/dy at each quadrature point
531 vec_dbl_Type u13(dPhiTrans.size(), -1.); // should correspond to du/dz at each quadrature point
532
533 vec_dbl_Type u21(dPhiTrans.size(), -1.); // should correspond to dv/dx at each quadrature point
534 vec_dbl_Type u22(dPhiTrans.size(), -1.); // should correspond to dv/dy at each quadrature point
535 vec_dbl_Type u23(dPhiTrans.size(), -1.); // should correspond to dv/dz at each quadrature point
536
537 vec_dbl_Type u31(dPhiTrans.size(), -1.); // should correspond to dw/dx at each quadrature point
538 vec_dbl_Type u32(dPhiTrans.size(), -1.); // should correspond to dw/dy at each quadrature point
539 vec_dbl_Type u33(dPhiTrans.size(), -1.); // should correspond to dw/dz at each quadrature point
540
541 vec_dbl_ptr_Type gammaDot(new vec_dbl_Type(weights->size(), 0.0)); // gammaDot->at(j) j=0...weights
542
543 vec_dbl_ptr_Type mixed_term_xy(new vec_dbl_Type(weights->size(), 0.0));
544 vec_dbl_ptr_Type mixed_term_xz(new vec_dbl_Type(weights->size(), 0.0));
545 vec_dbl_ptr_Type mixed_term_yz(new vec_dbl_Type(weights->size(), 0.0));
546
547 for (UN w = 0; w < dPhiTrans.size(); w++)
548 { // quads points
549
550 u11[w] = 0.0;
551 u12[w] = 0.0;
552 u13[w] = 0.0;
553 u21[w] = 0.0;
554 u22[w] = 0.0;
555 u23[w] = 0.0;
556 u31[w] = 0.0;
557 u32[w] = 0.0;
558 u33[w] = 0.0;
559
560 for (UN i = 0; i < dPhiTrans[0].size(); i++)
561 {
562 LO index1 = dim * i + 0; // x
563 LO index2 = dim * i + 1; // y
564 LO index3 = dim * i + 2; // z
565 // uLoc[d][w] += (*this->solution_)[index] * phi->at(w).at(i);
566 u11[w] += (*this->solution_)[index1] * dPhiTrans[w][i][0]; // u*dphi_dx
567 u12[w] += (*this->solution_)[index1] * dPhiTrans[w][i][1]; // because we are in 3D , 0 and 1, 2
568 u13[w] += (*this->solution_)[index1] * dPhiTrans[w][i][2];
569 u21[w] += (*this->solution_)[index2] * dPhiTrans[w][i][0]; // v*dphi_dx
570 u22[w] += (*this->solution_)[index2] * dPhiTrans[w][i][1];
571 u23[w] += (*this->solution_)[index2] * dPhiTrans[w][i][2];
572 u31[w] += (*this->solution_)[index3] * dPhiTrans[w][i][0]; // w*dphi_dx
573 u32[w] += (*this->solution_)[index3] * dPhiTrans[w][i][1];
574 u33[w] += (*this->solution_)[index3] * dPhiTrans[w][i][2];
575 }
576 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]));
577
578 mixed_term_xy->at(w) = 0.5 * (u12[w] + u21[w]);
579 mixed_term_xz->at(w) = 0.5 * (u31[w] + u13[w]);
580 mixed_term_yz->at(w) = 0.5 * (u32[w] + u23[w]);
581 }
582
583 // Initialize some helper vectors/matrices
584 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;
585
586 deta_dgamma_dgamma_dtau = 0.;
587
588 // Construct element matrices
589 for (UN i = 0; i < numNodes; i++)
590 {
591 // Teuchos::Array<SC> value(dPhiTrans[0].size(), 0. ); // dPhiTrans[0].size() is 3
592
593 for (UN j = 0; j < numNodes; j++)
594 {
595 // Reset values
596 v11 = 0.0;
597 v12 = 0.0;
598 v13 = 0.0;
599 v21 = 0.0;
600 v22 = 0.0;
601 v23 = 0.0;
602 v31 = 0.0;
603 v32 = 0.0;
604 v33 = 0.0;
605
606 // We compute for the derivative
607 for (UN w = 0; w < dPhiTrans.size(); w++)
608 {
609
610 value1_j = dPhiTrans.at(w).at(j).at(0); // so this corresponds to d\phi_j/dx
611 value2_j = dPhiTrans.at(w).at(j).at(1); // so this corresponds to d\phi_j/dy
612 value3_j = dPhiTrans.at(w).at(j).at(2); // so this corresponds to d\phi_j/dz
613
614 value1_i = dPhiTrans.at(w).at(i).at(0); // so this corresponds to d\phi_i/dx
615 value2_i = dPhiTrans.at(w).at(i).at(1); // so this corresponds to d\phi_i/dy
616 value3_i = dPhiTrans.at(w).at(i).at(2); // so this corresponds to d\phi_i/dz
617
618 this->viscosityModel->evaluateDerivative(this->params_, gammaDot->at(w), deta_dgamma_dgamma_dtau);
619
620 // ( -2.0 *deta/dgammaDot * dgammaDot/dTau * ( D(v^k) : D(\phi_i) ) * ( D(v^k) : D(\phi_j) ) ) ----> Symmetric tensor
621 // Construct entries - we go over all quadrature points and if j is updated we set v11 etc. again to zero
622 v11 = v11 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_i * u11[w] + value2_i * mixed_term_xy->at(w) + value3_i * mixed_term_xz->at(w)) * (value1_j * u11[w] + value2_j * mixed_term_xy->at(w) + value3_j * mixed_term_xz->at(w))) ;
623 v12 = v12 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_i * u11[w] + value2_i * mixed_term_xy->at(w) + value3_i * mixed_term_xz->at(w)) * (value1_j * mixed_term_xy->at(w) + u22[w] * value2_j + mixed_term_yz->at(w) * value3_j));
624 v13 = v13 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_i * u11[w] + value2_i * mixed_term_xy->at(w) + value3_i * mixed_term_xz->at(w)) * (value1_j * mixed_term_xz->at(w) + value2_j * mixed_term_yz->at(w) + value3_j * u33[w]));
625
626 v21 = v21 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_i * mixed_term_xy->at(w) + value2_i * u22[w] + value3_i * mixed_term_yz->at(w)) * (value1_j * u11[w] + value2_j * mixed_term_xy->at(w) + value3_j * mixed_term_xz->at(w)));
627 v22 = v22 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_i * mixed_term_xy->at(w) + value2_i * u22[w] + value3_i * mixed_term_yz->at(w)) * (value1_j * mixed_term_xy->at(w) + u22[w] * value2_j + mixed_term_yz->at(w) * value3_j));
628 v23 = v23 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_i * mixed_term_xy->at(w) + value2_i * u22[w] + value3_i * mixed_term_yz->at(w)) * (value1_j * mixed_term_xz->at(w) + value2_j * mixed_term_yz->at(w) + value3_j * u33[w]));
629
630 v31 = v31 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_i * mixed_term_xz->at(w) + mixed_term_yz->at(w) * value2_i + u33[w] * value3_i) * (value1_j * u11[w] + value2_j * mixed_term_xy->at(w) + value3_j * mixed_term_xz->at(w)));
631 v32 = v32 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_i * mixed_term_xz->at(w) + mixed_term_yz->at(w) * value2_i + u33[w] * value3_i) * (value1_j * mixed_term_xy->at(w) + u22[w] * value2_j + mixed_term_yz->at(w) * value3_j));
632 v33 = v33 + (-2.0) * deta_dgamma_dgamma_dtau * weights->at(w) * ((value1_i * mixed_term_xz->at(w) + mixed_term_yz->at(w) * value2_i + u33[w] * value3_i) * (value1_j * mixed_term_xz->at(w) + value2_j * mixed_term_yz->at(w) + value3_j * u33[w]));
633
634 } // loop end quadrature points
635
636 // multiply determinant from transformation
637 v11 *= absDetB;
638 v12 *= absDetB;
639 v13 *= absDetB;
640 v21 *= absDetB;
641 v22 *= absDetB;
642 v23 *= absDetB;
643 v31 *= absDetB;
644 v32 *= absDetB;
645 v33 *= absDetB;
646
647 // Put values on the right position in element matrix
648 // [v11 v12 v13]
649 // [v21 v22 v23]
650 // [v31 v32 v33]
651 (*elementMatrix)[i * dofs][j * dofs] = v11; // d=0, first dimension
652 (*elementMatrix)[i * dofs][j * dofs + 1] = v12;
653 (*elementMatrix)[i * dofs][j * dofs + 2] = v13;
654 (*elementMatrix)[i * dofs + 1][j * dofs] = v21;
655 (*elementMatrix)[i * dofs + 1][j * dofs + 1] = v22; // d=1, second dimension
656 (*elementMatrix)[i * dofs + 1][j * dofs + 2] = v23; // d=1, second dimension
657 (*elementMatrix)[i * dofs + 2][j * dofs] = v31;
658 (*elementMatrix)[i * dofs + 2][j * dofs + 1] = v32; // d=2, third dimension
659 (*elementMatrix)[i * dofs + 2][j * dofs + 2] = v33; // d=2, third dimension
660
661 } // loop end over j node
662
663 } // loop end over i node
664
665 } // end if dim = 3
666 }
667
668
669 //ToDo: Using new functionalty integrate again assembly of Neumann boundary term
670
671 // "Fixpunkt"- Matrix without jacobian for calculating Ax
672 // Here update please to unlinearized System Matrix accordingly.
673 template <class SC, class LO, class GO, class NO>
674 void AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::assembleRHS()
675 {
676
677 SmallMatrixPtr_Type elementMatrixN = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
678 SmallMatrixPtr_Type elementMatrixNC = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
679
680 this->ANB_.reset(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_)); // A + B + N
681 this->ANB_->add((*this->constantMatrix_), (*this->ANB_));
682
683 // Nonlinear shear stress tensor *******************************
684 this->assemblyStress(elementMatrixN);
685 this->ANB_->add((*elementMatrixN), (*this->ANB_));
686
687 // Nonlinear convection term *******************************
688 this->assemblyAdvection(elementMatrixNC);
689 elementMatrixNC->scale(this->density_);
690 this->ANB_->add((*elementMatrixNC), (*this->ANB_));
691
692 // TODO: If underlying element is an outflow boundary element - will be added when func nonlinear boundar term *******************************
693 /*if (this->surfaceElement == true)
694 {
695 SmallMatrixPtr_Type elementMatrixNB = Teuchos::rcp(new SmallMatrix_Type(this->dofsElementVelocity_ + this->numNodesPressure_));
696 this->assemblyNeumannBoundaryTerm(elementMatrixNB);
697 this->ANB_->add((*elementMatrixNB), ((*this->ANB_)));
698 }
699 */
700
701 this->rhsVec_.reset(new vec_dbl_Type(this->dofsElement_, 0.));
702 // Multiplying ANB_ * solution // ANB Matrix without nonlinear part.
703 int s = 0, t = 0;
704 for (int i = 0; i < this->ANB_->size(); i++)
705 {
706 if (i >= this->dofsElementVelocity_)
707 s = 1;
708 for (int j = 0; j < this->ANB_->size(); j++)
709 {
710 if (j >= this->dofsElementVelocity_)
711 t = 1;
712 (*this->rhsVec_)[i] += (*this->ANB_)[i][j] * (*this->solution_)[j] * this->coeff_[s][t];
713 }
714 t = 0;
715 }
716 }
717
721 /* In 2D B=[ nodesRefConfig(1).at(0)-nodesRefConfig(0).at(0) nodesRefConfig(2).at(0)-nodesRefConfig(0).at(0) ]
722 [ nodesRefConfig(1).at(1)-nodesRefConfig(0).at(1) nodesRefConfig(2).at(1)-nodesRefConfig(0).at(1) ]
723 */
738
739
740 // Compute Shear Rate on quadrature points depending on gradient of velocity solution at nodes
741 template <class SC, class LO, class GO, class NO>
742 void AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::computeShearRate(vec3D_dbl_Type dPhiTrans,
743 vec_dbl_ptr_Type &gammaDot, int dim)
744 {
745
746 //****************** TWO DIMENSIONAL *********************************
747 if (dim == 2)
748 {
749
750 vec_dbl_Type u11(dPhiTrans.size(), -1.); // should correspond to du/dx at each quadrature point
751 vec_dbl_Type u12(dPhiTrans.size(), -1.); // should correspond to du/dy at each quadrature point
752 vec_dbl_Type u21(dPhiTrans.size(), -1.); // should correspond to dv/dx at each quadrature point
753 vec_dbl_Type u22(dPhiTrans.size(), -1.); // should correspond to dv/dy at each quadrature point
754
755 for (UN w = 0; w < dPhiTrans.size(); w++)
756 { // quads points
757 // set again to zero
758 u11[w] = 0.0;
759 u12[w] = 0.0;
760 u21[w] = 0.0;
761 u22[w] = 0.0;
762 for (UN i = 0; i < dPhiTrans[0].size(); i++)
763 { // loop unrolling
764 LO index1 = dim * i + 0; // x
765 LO index2 = dim * i + 1; // y
766 // uLoc[d][w] += (*this->solution_)[index] * phi->at(w).at(i);
767 u11[w] += (*this->solution_)[index1] * dPhiTrans[w][i][0]; // u*dphi_dx
768 u12[w] += (*this->solution_)[index1] * dPhiTrans[w][i][1]; // because we are in 2D , 0 and 1
769 u21[w] += (*this->solution_)[index2] * dPhiTrans[w][i][0];
770 u22[w] += (*this->solution_)[index2] * dPhiTrans[w][i][1];
771 }
772 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]));
773 }
774 } // end if dim == 2
775 //****************** THREE DIMENSIONAL *********************************
776 else if (dim == 3)
777 {
778
779 vec_dbl_Type u11(dPhiTrans.size(), -1.); // should correspond to du/dx at each quadrature point
780 vec_dbl_Type u12(dPhiTrans.size(), -1.); // should correspond to du/dy at each quadrature point
781 vec_dbl_Type u13(dPhiTrans.size(), -1.); // should correspond to du/dz at each quadrature point
782
783 vec_dbl_Type u21(dPhiTrans.size(), -1.); // should correspond to dv/dx at each quadrature point
784 vec_dbl_Type u22(dPhiTrans.size(), -1.); // should correspond to dv/dy at each quadrature point
785 vec_dbl_Type u23(dPhiTrans.size(), -1.); // should correspond to dv/dz at each quadrature point
786
787 vec_dbl_Type u31(dPhiTrans.size(), -1.); // should correspond to dw/dx at each quadrature point
788 vec_dbl_Type u32(dPhiTrans.size(), -1.); // should correspond to dw/dy at each quadrature point
789 vec_dbl_Type u33(dPhiTrans.size(), -1.); // should correspond to dw/dz at each quadrature point
790
791 for (UN w = 0; w < dPhiTrans.size(); w++)
792 { // quads points
793 // set again to zero
794 u11[w] = 0.0;
795 u12[w] = 0.0;
796 u13[w] = 0.0;
797 u21[w] = 0.0;
798 u22[w] = 0.0;
799 u23[w] = 0.0;
800 u31[w] = 0.0;
801 u32[w] = 0.0;
802 u33[w] = 0.0;
803
804 for (UN i = 0; i < dPhiTrans[0].size(); i++)
805 {
806 LO index1 = dim * i + 0; // x
807 LO index2 = dim * i + 1; // y
808 LO index3 = dim * i + 2; // z
809 // uLoc[d][w] += (*this->solution_)[index] * phi->at(w).at(i);
810 u11[w] += (*this->solution_)[index1] * dPhiTrans[w][i][0]; // u*dphi_dx
811 u12[w] += (*this->solution_)[index1] * dPhiTrans[w][i][1]; // because we are in 3D , 0 and 1, 2
812 u13[w] += (*this->solution_)[index1] * dPhiTrans[w][i][2];
813 u21[w] += (*this->solution_)[index2] * dPhiTrans[w][i][0]; // v*dphi_dx
814 u22[w] += (*this->solution_)[index2] * dPhiTrans[w][i][1];
815 u23[w] += (*this->solution_)[index2] * dPhiTrans[w][i][2];
816 u31[w] += (*this->solution_)[index3] * dPhiTrans[w][i][0]; // w*dphi_dx
817 u32[w] += (*this->solution_)[index3] * dPhiTrans[w][i][1];
818 u33[w] += (*this->solution_)[index3] * dPhiTrans[w][i][2];
819 }
820 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]));
821 }
822 } // end if dim == 3
823 }
824
825 /* Based on the current solution (velocity, pressure etc.) we want to be able to compute postprocessing fields
826 like here the viscosity inside an element.
827 */
828 template <class SC, class LO, class GO, class NO>
829 void AssembleFEGeneralizedNewtonian<SC, LO, GO, NO>::computeLocalconstOutputField()
830 {
831 int dim = this->getDim();
832 std::string FEType = this->FETypeVelocity_;
833
834 SC detB;
835 SmallMatrix<SC> B(dim);
836 SmallMatrix<SC> Binv(dim);
837
838 this->buildTransformation(B);
839 detB = B.computeInverse(Binv);
840
841 vec3D_dbl_ptr_Type dPhiAtCM;
842
843 // Compute viscosity at center of mass using nodal values and shape function **********************************************************************************
844 TEUCHOS_TEST_FOR_EXCEPTION(dim == 1, std::logic_error, "computeLocalconstOutputField Not implemented for dim=1");
845
846 Helper::getDPhiAtCM(dPhiAtCM, dim, FEType); // These are the original coordinates of the reference element
847 vec3D_dbl_Type dPhiTransAtCM(dPhiAtCM->size(), vec2D_dbl_Type(dPhiAtCM->at(0).size(), vec_dbl_Type(dim, 0.)));
848 Helper::applyBTinv(dPhiAtCM, dPhiTransAtCM, Binv); // We need transformation because of velocity gradient in shear rate equation
849
850 vec_dbl_ptr_Type gammaDoti(new vec_dbl_Type(dPhiAtCM->size(), 0.0)); // Only one value because size is one
851 computeShearRate(dPhiTransAtCM, gammaDoti, dim); // updates gammaDot using velocity solution
852 this->viscosityModel->evaluateMapping(this->params_, gammaDoti->at(0), this->constOutputField_.at(0));
853 }
854
855}
856#endif
FEDD::AssembleFEGeneralizedNewtonian::assembleRHS
void assembleRHS() override
Assemble the element right hand side vector.
Definition AssembleFEGeneralizedNewtonian_def.hpp:674
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:829
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:742
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