dune-composites (2.5.1)

getStress.hh
1#ifndef GET_STRESS_HH
2#define GET_STRESS_HH
3
4#include <dune/geometry/type.hh>
5#include <dune/geometry/referenceelements.hh>
6#include <dune/geometry/quadraturerules.hh>
7
8namespace Dune {
9 namespace PDELab {
10
11 template<typename GV, typename PROBLEM, typename DGF, int dofel>
12 class getStressLinear :
21 public NumericalJacobianApplyVolume<getStressLinear<GV,PROBLEM,DGF,dofel> >,
22 public FullVolumePattern,
23 public LocalOperatorDefaultFlags,
24 public InstationaryLocalOperatorDefaultMethods<double>,
25 public NumericalJacobianVolume<getStressLinear<GV,PROBLEM,DGF, dofel> >
26 {
27 public:
28 // pattern assembly flags
29 enum { doPatternVolume = false };
30
31 // residual assembly flags
32 enum { doAlphaVolume = true };
33
34 //Constructor
35 getStressLinear (const GV& gv_, const PROBLEM& prob_, const DGF& dgf_, double& maxDisp_) : gv(gv_), problem(prob_) , maxDisp(maxDisp_), dgf(dgf_){}
36 //
37 // alpha_volume for getStress
38 //
39 // Volume Integral Depending on Test and Ansatz Functions
40 template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
41 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
42 {
43 const int dim = 3;
44
45 //Unwrap function spaces
46 typedef typename LFSU::template Child<0>::Type LFS1;
47 typedef typename LFS1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::JacobianType JacobianType;
48 //Unwrap displacements
49 const auto& lfs1 = lfsu.template child<0>();
50 const auto& lfs2 = lfsu.template child<1>();
51 const auto& lfs3 = lfsu.template child<2>();
52 //Unwrap stress
53 const auto& lfsv1 = lfsv.template child<0>();
54 const auto& lfsv2 = lfsv.template child<1>();
55 const auto& lfsv3 = lfsv.template child<2>();
56 const auto& lfsv4 = lfsv.template child<3>();
57 const auto& lfsv5 = lfsv.template child<4>();
58 const auto& lfsv6 = lfsv.template child<5>();
59
60 //Evaluate elasticity tensor
61 Dune::FieldMatrix<double,6,6> C(0.0),Cthermal(0.0);
62 problem.evaluate(eg.entity(),C);
63 problem.evaluate_strain(eg.entity(),Cthermal);
64
65 //Residual heat stresses
66 Dune::FieldVector<double,6> delta(0.0);
67 if(problem.getMaterialType(eg.entity()) == 0)//resin
68 problem.evaluateHeat(delta,true);
69 else
70 problem.evaluateHeat(delta,false); //ply
71
72 //Calculate input heat stress
73 Dune::FieldVector<double,6> appliedStress(0.0);
74 Cthermal.mv(delta, appliedStress);
75
76 int nodes_per_element = lfsv1.size();
77 int dof = lfs1.size();
78
79 // select quadrature rule
80 for (int iter = 0; iter < nodes_per_element; iter++){
81 Dune::FieldVector<double,1> num_elem(1.0);
82 Dune::FieldVector<double,dim> pt = eg.geometry().local(eg.geometry().corner(iter));
83 Dune::FieldVector<double,dim> pt_cen = eg.geometry().local(eg.geometry().center());
84
85 auto theta = problem.getThetaFct(eg.entity(),pt);
86 auto Cijkl = getMatrix(C,theta);
87
88 // Evaluate Jacobian
89 std::vector<JacobianType> js(dof);
90 lfs1.finiteElement().localBasis().evaluateJacobian(pt,js);
91
92 // Transform gradient to real element and evaluate at integration points
93 const auto jac = eg.geometry().jacobianInverseTransposed(pt);
94 std::vector<Dune::FieldVector<double,dim> > gradphi(dof);
95 for (int i=0; i < dof; i++){
96 gradphi[i] = 0.0;
97 jac.umv(js[i][0],gradphi[i]);
98 }
99 Dune::FieldMatrix<double,6,dofel> B(0.0);
100 for (int i = 0; i < dof; i++)
101 {
102 B[0][i] = gradphi[i][0]; // E11
103 B[1][i + dof] = gradphi[i][1]; // E22
104 B[2][i + 2 * dof] = gradphi[i][2]; // E33
105 B[3][i + dof] = gradphi[i][2]; B[3][i + 2 * dof] = gradphi[i][1]; // E23
106 B[4][i] = gradphi[i][2]; B[4][i + 2 * dof] = gradphi[i][0]; // E13
107 B[5][i] = gradphi[i][1]; B[5][i + dof] = gradphi[i][0]; // E12
108 }
109
110 Dune::FieldVector<double,6> e(0.0),s(0.0);
111 Dune::FieldVector<double,dofel> xlc(0.0);
112 for(int i = 0; i < dof; i++){
113 xlc[i] = x(lfs1,i);
114 xlc[i+dof] = x(lfs2,i);
115 xlc[i+2*dof] = x(lfs3,i);
116 }
117
118 B.mv(xlc,e); // e = B * x
119 Cijkl.mv(e,s); // s = Cijkl * e
120
121 s -= appliedStress;
122
123 // Evaluate Jacobian
124 std::vector<JacobianType> js_cen(dof);
125 lfs1.finiteElement().localBasis().evaluateJacobian(pt_cen,js_cen);
126
127 // Transform gradient to real element and evaluate at integration points
128 const auto jac_cen = eg.geometry().jacobianInverseTransposed(pt_cen);
129 std::vector<Dune::FieldVector<double,dim> > gradphi_cen(dof);
130 for (int i=0; i < dof; i++){
131 gradphi_cen[i] = 0.0;
132 jac_cen.umv(js_cen[i][0],gradphi_cen[i]);
133 }
134 Dune::FieldMatrix<double,6,dofel> B_cen(0.0);
135 for (int i = 0; i < dof; i++)
136 {
137 B_cen[0][i] = gradphi_cen[i][0]; // E11
138 B_cen[1][i + dof] = gradphi_cen[i][1]; // E22
139 B_cen[2][i + 2 * dof] = gradphi_cen[i][2]; // E33
140 B_cen[3][i + dof] = gradphi_cen[i][2]; B_cen[3][i + 2 * dof] = gradphi_cen[i][1]; // E23
141 B_cen[4][i] = gradphi_cen[i][2]; B_cen[4][i + 2 * dof] = gradphi_cen[i][0]; // E13
142 B_cen[5][i] = gradphi_cen[i][1]; B_cen[5][i + dof] = gradphi_cen[i][0]; // E12
143 }
144
145 Dune::FieldVector<double,6> e_cen(0.0),s_cen(0.0);
146 B_cen.mv(xlc,e_cen); // e = B * x
147 Cijkl.mv(e_cen,s_cen); // s = Cijkl * e
148
149 s_cen -= appliedStress;
150
151 for(int i = 0; i < dof; i++){
152 double maxUpdate = sqrt(xlc[i]*xlc[i]+xlc[i+dof]*xlc[i+dof]+xlc[i+2*dof]*xlc[i+2*dof]);
153 if(maxUpdate > maxDisp){ maxDisp = maxUpdate; }
154 }
155
156 //Stresses
157 dgf.evaluate(eg.entity(),pt,num_elem);
158 r.accumulate(lfsv1, iter, s[0]);
159 r.accumulate(lfsv2, iter, s[1] /num_elem[0]);
160 r.accumulate(lfsv3, iter, s_cen[2]/num_elem[0]);
161 r.accumulate(lfsv4, iter, s[3]);
162 r.accumulate(lfsv5, iter, s[4] /num_elem[0]);
163 r.accumulate(lfsv6, iter, s[5]);
164 }
165
166
167 }
168 private:
169
170 Dune::FieldMatrix<double,6,6> getR(Dune::FieldMatrix<double,3,3> A) const{
171 Dune::FieldMatrix<double,6,6> R;
172 R[0][0] = A[0][0]*A[0][0]; R[0][1] = A[0][1]*A[0][1]; R[0][2] = A[0][2]*A[0][2];
173 R[1][0] = A[1][0]*A[1][0]; R[1][1] = A[1][1]*A[1][1]; R[1][2] = A[1][2]*A[1][2];
174 R[2][0] = A[2][0]*A[2][0]; R[2][1] = A[2][1]*A[2][1]; R[2][2] = A[2][2]*A[2][2];
175
176 R[0][3] = 2.0*A[0][1]*A[0][2]; R[0][4] = 2.0*A[0][0]*A[0][2]; R[0][5] = 2.0*A[0][0]*A[0][1];
177 R[1][3] = 2.0*A[1][1]*A[1][2]; R[1][4] = 2.0*A[1][0]*A[1][2]; R[1][5] = 2.0*A[1][0]*A[1][1];
178 R[2][3] = 2.0*A[2][1]*A[2][2]; R[2][4] = 2.0*A[2][0]*A[2][2]; R[2][5] = 2.0*A[2][0]*A[2][1];
179
180 R[3][0] = A[1][0]*A[2][0]; R[3][1] = A[1][1]*A[2][1]; R[3][2] = A[1][2]*A[2][2];
181 R[4][0] = A[0][0]*A[2][0]; R[4][1] = A[0][1]*A[2][1]; R[4][2] = A[0][2]*A[2][2];
182 R[5][0] = A[0][0]*A[1][0]; R[5][1] = A[0][1]*A[1][1]; R[5][2] = A[0][2]*A[1][2];
183
184 R[3][3] = A[1][1]*A[2][2]+A[1][2]*A[2][1]; R[3][4] = A[1][0]*A[2][2]+A[1][2]*A[2][0]; R[3][5] = A[1][0]*A[2][1]+A[1][1]*A[2][0];
185 R[4][3] = A[0][1]*A[2][2]+A[2][1]*A[0][2]; R[4][4] = A[0][0]*A[2][2]+A[2][0]*A[0][2]; R[4][5] = A[2][0]*A[0][1]+A[2][1]*A[0][0];
186 R[5][3] = A[0][1]*A[1][2]+A[0][2]*A[1][1]; R[5][4] = A[0][0]*A[1][2]+A[0][2]*A[1][0]; R[5][5] = A[0][0]*A[1][1]+A[0][1]*A[1][0];
187 return R;
188 }
189
190 Dune::FieldMatrix<double,6,6> getMatrix(Dune::FieldMatrix<double,6,6> C, std::vector<double> theta) const{
191 double c = cos(theta[1]*M_PI/180.0);
192 double s = sin(theta[1]*M_PI/180.0);
193 Dune::FieldMatrix<double,3,3> A = {{c,0,-s},{0,1,0},{s,0,c}};
194
195 Dune::FieldMatrix<double,6,6> R = getR(A);
196
197 auto Ch = C;
198 Dune::FieldMatrix<double,6,6> RT(0.0);
199 Composites::transpose<Dune::FieldMatrix<double,6,6>,6>(R,RT);
200 Ch.rightmultiply(RT);
201 return Ch;
202 }
203
204 const GV& gv;
205 const PROBLEM& problem;
206 const DGF& dgf;
207 double& maxDisp;
208 };
209
210
211 template<typename GV, typename MODEL, int dofel>
212 class getStress :
213 public NumericalJacobianApplyVolume<getStress<GV,MODEL,dofel> >,
214 public FullVolumePattern,
215 public LocalOperatorDefaultFlags,
216 public InstationaryLocalOperatorDefaultMethods<double>,
217 public NumericalJacobianVolume<getStress<GV,MODEL, dofel> >
218 {
219 public:
220 // pattern assembly flags
221 enum { doPatternVolume = false };
222
223 // residual assembly flags
224 enum { doAlphaVolume = true };
225
226 //Constructor
227 getStress (const GV& gv_, const MODEL& model_) : gv(gv_), model(model_) {}
228 //
229 // alpha_volume for getStress
230 //
231 // Volume Integral Depending on Test and Ansatz Functions
232 template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
233 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
234 {
235 const int dim = 3;
236
237 //Unwrap function spaces
238 typedef typename LFSU::template Child<0>::Type LFS1;
239 typedef typename LFS1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::DomainFieldType DF;
240 typedef typename LFS1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::JacobianType JacobianType;
241 //Unwrap displacements
242 const auto& lfs1 = lfsu.template child<0>();
243 const auto& lfs2 = lfsu.template child<1>();
244 const auto& lfs3 = lfsu.template child<2>();
245 //Unwrap stress
246 const auto& lfsv1 = lfsv.template child<0>();
247 const auto& lfsv2 = lfsv.template child<1>();
248 const auto& lfsv3 = lfsv.template child<2>();
249 const auto& lfsv4 = lfsv.template child<3>();
250 const auto& lfsv5 = lfsv.template child<4>();
251 const auto& lfsv6 = lfsv.template child<5>();
252
253 // Get element id
254 int id = gv.indexSet().index(eg.entity());
255
256 //Evaluate elasticity tensor
257 Dune::FieldMatrix<double,6,6> C = model.getElasticTensor(id);
258 auto theta = model.returnTheta(eg.geometry().center());
259 auto Cthermal = model.TensorRotationRight(C,-theta[2],2);
260
261 int nodes_per_element = lfsv1.size();
262 int dof = lfs1.size();
263
264 //Residual heat stresses
265 Dune::FieldVector<double,6> delta(0.0);
266 model.evaluateHeat(delta, id);
267
268 //Calculate input heat stress
269 Dune::FieldVector<double,6> appliedStress(0.0);
270 Cthermal.mv(delta, appliedStress);
271
272 // select quadrature rule
273 GeometryType gt = eg.geometry().type();
274 const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,1);
275 for (typename Dune::QuadratureRule<DF,dim>::const_iterator it = rule.begin(),endit = rule.end(); it != endit; ++it){
276 Dune::FieldVector<double,dim> pt = it->position();
277 Dune::FieldVector<DF,dim> local = eg.geometry().global(pt);
278
279 auto theta = model.returnTheta(local);
280 auto Cijkl = model.TensorRotationRight(C,theta[1],1);
281
282 // Evaluate Jacobian
283 std::vector<JacobianType> js(dof);
284 lfs1.finiteElement().localBasis().evaluateJacobian(pt,js);
285
286 // Transform gradient to real element and evaluate at element center
287 const auto jac = eg.geometry().jacobianInverseTransposed(pt);
288 std::vector<Dune::FieldVector<double,dim> > gradphi(dof);
289 for (int i=0; i < dof; i++){
290 gradphi[i] = 0.0;
291 jac.umv(js[i][0],gradphi[i]);
292 }
293 Dune::FieldMatrix<double,6,dofel> B(0.0);
294 for (int i = 0; i < dof; i++)
295 {
296 B[0][i] = gradphi[i][0]; // E11
297 B[1][i + dof] = gradphi[i][1]; // E22
298 B[2][i + 2 * dof] = gradphi[i][2]; // E33
299 B[3][i + dof] = gradphi[i][2]; B[3][i + 2 * dof] = gradphi[i][1]; // E23
300 B[4][i] = gradphi[i][2]; B[4][i + 2 * dof] = gradphi[i][0]; // E13
301 B[5][i] = gradphi[i][1]; B[5][i + dof] = gradphi[i][0]; // E12
302 }
303
304 Dune::FieldVector<double,6> e(0.0),s(0.0);
305 Dune::FieldVector<double,dofel> xlc(0.0);
306 for(int i = 0; i < dof; i++){
307 xlc[i] = x(lfs1,i);
308 xlc[i+dof] = x(lfs2,i);
309 xlc[i+2*dof] = x(lfs3,i);
310 }
311
312 B.mv(xlc,e); // e = B * x
313 Cijkl.mv(e,s); // s = Cijkl * e
314
315 s-=appliedStress;
316
317 for (int i = 0; i < nodes_per_element; i++){
318 r.accumulate(lfsv1, i, s[0] * it->weight());
319 r.accumulate(lfsv2, i, s[1] * it->weight());
320 r.accumulate(lfsv3, i, s[2] * it->weight());
321 r.accumulate(lfsv4, i, s[3] * it->weight());
322 r.accumulate(lfsv5, i, s[4] * it->weight());
323 r.accumulate(lfsv6, i, s[5] * it->weight());
324 }
325 }
326 }
327 const GV& gv;
328 const MODEL& model;
329
330 };
331
332 class countElem :
333 public NumericalJacobianApplyVolume<countElem >,
334 public FullVolumePattern,
335 public LocalOperatorDefaultFlags,
336 public InstationaryLocalOperatorDefaultMethods<double>,
337 public NumericalJacobianVolume<countElem >
338 {
339 public:
340 // pattern assembly flags
341 enum { doPatternVolume = false };
342
343 // residual assembly flags
344 enum { doLambdaVolume = true };
345
346 //Constructor
347 countElem (){}
348
349 // alpha_volume for countElem
350 template<typename EG, typename LFSV, typename R>
351 void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
352 {
353 for (int i = 0; i < lfsv.size(); i++){
354 r.accumulate(lfsv, i, 1.);
355 }
356 }
357 };
358
359
360 template<typename GV, typename PROBLEM, typename DGF, int dofel>
361 class getStressUG :
362 public NumericalJacobianApplyVolume<getStressUG<GV,PROBLEM,DGF,dofel> >,
363 public FullVolumePattern,
364 public LocalOperatorDefaultFlags,
365 public InstationaryLocalOperatorDefaultMethods<double>,
366 public NumericalJacobianVolume<getStressUG<GV,PROBLEM,DGF, dofel> >
367 {
368 public:
369 // pattern assembly flags
370 enum { doPatternVolume = false };
371
372 // residual assembly flags
373 enum { doAlphaVolume = true };
374
375 //Constructor
376 getStressUG (const GV& gv_, const PROBLEM& prob_, const DGF& dgf_, double& maxDisp_) : gv(gv_), problem(prob_) , maxDisp(maxDisp_), dgf(dgf_){}
377 //
378 // alpha_volume for getStress
379 //
380 // Volume Integral Depending on Test and Ansatz Functions
381 template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
382 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
383 {
384 const int dim = 3;
385
386 //Unwrap function spaces
387 typedef typename LFSU::template Child<0>::Type LFS1;
388 typedef typename LFS1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::JacobianType JacobianType;
389 //Unwrap displacements
390 const auto& lfs1 = lfsu.template child<0>();
391 const auto& lfs2 = lfsu.template child<1>();
392 const auto& lfs3 = lfsu.template child<2>();
393 //Unwrap stress
394 const auto& lfsv1 = lfsv.template child<0>();
395 const auto& lfsv2 = lfsv.template child<1>();
396 const auto& lfsv3 = lfsv.template child<2>();
397 const auto& lfsv4 = lfsv.template child<3>();
398 const auto& lfsv5 = lfsv.template child<4>();
399 const auto& lfsv6 = lfsv.template child<5>();
400
401 //Evaluate elasticity tensor
402 Dune::FieldMatrix<double,6,6> C(0.0),Cthermal(0.0);
403 problem.evaluate(eg.entity(),C);
404 problem.evaluate_strain(eg.entity(),Cthermal);
405
406 //Residual heat stresses
407 Dune::FieldVector<double,6> delta(0.0);
408 if(problem.getMaterialType(eg.entity()) == 0)//resin
409 problem.evaluateHeat(delta,true);
410 else
411 problem.evaluateHeat(delta,false); //ply
412
413 //Calculate input heat stress
414 Dune::FieldVector<double,6> appliedStress(0.0);
415 Cthermal.mv(delta, appliedStress);
416
417 int nodes_per_element = lfsv1.size();
418 int dof = lfs1.size();
419
420 // select quadrature rule
421 for (int iter = 0; iter < nodes_per_element; iter++){
422 Dune::FieldVector<double,1> num_elem(1.0);
423 Dune::FieldVector<double,dim> pt = eg.geometry().local(eg.geometry().corner(iter));
424
425 auto theta = problem.getThetaFct(eg.entity(),pt);
426 auto Cijkl = getMatrix(C,theta);
427
428 // Evaluate Jacobian
429 std::vector<JacobianType> js(dof);
430 lfs1.finiteElement().localBasis().evaluateJacobian(pt,js);
431
432 // Transform gradient to real element and evaluate at integration points
433 const auto jac = eg.geometry().jacobianInverseTransposed(pt);
434 std::vector<Dune::FieldVector<double,dim> > gradphi(dof);
435 for (int i=0; i < dof; i++){
436 gradphi[i] = 0.0;
437 jac.umv(js[i][0],gradphi[i]);
438 }
439 Dune::FieldMatrix<double,6,dofel> B(0.0);
440 for (int i = 0; i < dof; i++)
441 {
442 B[0][i] = gradphi[i][0]; // E11
443 B[1][i + dof] = gradphi[i][1]; // E22
444 B[2][i + 2 * dof] = gradphi[i][2]; // E33
445 B[3][i + dof] = gradphi[i][2]; B[3][i + 2 * dof] = gradphi[i][1]; // E23
446 B[4][i] = gradphi[i][2]; B[4][i + 2 * dof] = gradphi[i][0]; // E13
447 B[5][i] = gradphi[i][1]; B[5][i + dof] = gradphi[i][0]; // E12
448 }
449
450 Dune::FieldVector<double,6> e(0.0),s(0.0);
451 Dune::FieldVector<double,dofel> xlc(0.0);
452 for(int i = 0; i < dof; i++){
453 xlc[i] = x(lfs1,i);
454 xlc[i+dof] = x(lfs2,i);
455 xlc[i+2*dof] = x(lfs3,i);
456 }
457
458 B.mv(xlc,e); // e = B * x
459 Cijkl.mv(e,s); // s = Cijkl * e
460
461 s -= appliedStress;
462
463 for(int i = 0; i < dof; i++){
464 double maxUpdate = sqrt(xlc[i]*xlc[i]+xlc[i+dof]*xlc[i+dof]+xlc[i+2*dof]*xlc[i+2*dof]);
465 if(maxUpdate > maxDisp){ maxDisp = maxUpdate; }
466 }
467
468 //Stresses
469 dgf.evaluate(eg.entity(),pt,num_elem);
470 r.accumulate(lfsv1, iter, s[0]/num_elem[0]);
471 r.accumulate(lfsv2, iter, s[1]/num_elem[0]);
472 r.accumulate(lfsv3, iter, s[2]/num_elem[0]);
473 r.accumulate(lfsv4, iter, s[3]/num_elem[0]);
474 r.accumulate(lfsv5, iter, s[4]/num_elem[0]);
475 r.accumulate(lfsv6, iter, s[5]/num_elem[0]);
476 }
477
478
479 }
480 private:
481
482 Dune::FieldMatrix<double,6,6> getR(Dune::FieldMatrix<double,3,3> A) const{
483 Dune::FieldMatrix<double,6,6> R;
484 R[0][0] = A[0][0]*A[0][0]; R[0][1] = A[0][1]*A[0][1]; R[0][2] = A[0][2]*A[0][2];
485 R[1][0] = A[1][0]*A[1][0]; R[1][1] = A[1][1]*A[1][1]; R[1][2] = A[1][2]*A[1][2];
486 R[2][0] = A[2][0]*A[2][0]; R[2][1] = A[2][1]*A[2][1]; R[2][2] = A[2][2]*A[2][2];
487
488 R[0][3] = 2.0*A[0][1]*A[0][2]; R[0][4] = 2.0*A[0][0]*A[0][2]; R[0][5] = 2.0*A[0][0]*A[0][1];
489 R[1][3] = 2.0*A[1][1]*A[1][2]; R[1][4] = 2.0*A[1][0]*A[1][2]; R[1][5] = 2.0*A[1][0]*A[1][1];
490 R[2][3] = 2.0*A[2][1]*A[2][2]; R[2][4] = 2.0*A[2][0]*A[2][2]; R[2][5] = 2.0*A[2][0]*A[2][1];
491
492 R[3][0] = A[1][0]*A[2][0]; R[3][1] = A[1][1]*A[2][1]; R[3][2] = A[1][2]*A[2][2];
493 R[4][0] = A[0][0]*A[2][0]; R[4][1] = A[0][1]*A[2][1]; R[4][2] = A[0][2]*A[2][2];
494 R[5][0] = A[0][0]*A[1][0]; R[5][1] = A[0][1]*A[1][1]; R[5][2] = A[0][2]*A[1][2];
495
496 R[3][3] = A[1][1]*A[2][2]+A[1][2]*A[2][1]; R[3][4] = A[1][0]*A[2][2]+A[1][2]*A[2][0]; R[3][5] = A[1][0]*A[2][1]+A[1][1]*A[2][0];
497 R[4][3] = A[0][1]*A[2][2]+A[2][1]*A[0][2]; R[4][4] = A[0][0]*A[2][2]+A[2][0]*A[0][2]; R[4][5] = A[2][0]*A[0][1]+A[2][1]*A[0][0];
498 R[5][3] = A[0][1]*A[1][2]+A[0][2]*A[1][1]; R[5][4] = A[0][0]*A[1][2]+A[0][2]*A[1][0]; R[5][5] = A[0][0]*A[1][1]+A[0][1]*A[1][0];
499 return R;
500 }
501
502 Dune::FieldMatrix<double,6,6> getMatrix(Dune::FieldMatrix<double,6,6> C, std::vector<double> theta) const{
503 double c = cos(theta[1]*M_PI/180.0);
504 double s = sin(theta[1]*M_PI/180.0);
505 Dune::FieldMatrix<double,3,3> A = {{c,0,-s},{0,1,0},{s,0,c}};
506
507 Dune::FieldMatrix<double,6,6> R = getR(A);
508
509 auto Ch = C;
510 Dune::FieldMatrix<double,6,6> RT(0.0);
511 Composites::transpose<Dune::FieldMatrix<double,6,6>,6>(R,RT);
512 Ch.rightmultiply(RT);
513 return Ch;
514 }
515
516 const GV& gv;
517 const PROBLEM& problem;
518 const DGF& dgf;
519 double& maxDisp;
520 };
521
522 }
523}
524
525#endif
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Sep 4, 22:38, 2025)