dune-composites (2.5.1)

linearelasticity.hh
1#ifndef DUNE_PDELAB_LINEARELASTIC_HH
2#define DUNE_PDELAB_LINEARELASTIC_HH
3
4#include "../../MathFunctions/transpose.hh"
5
6#include "defaultimp.hh"
7#include "pattern.hh"
8#include "flags.hh"
9#include "idefault.hh"
10
11namespace Dune {
12 namespace PDELab {
13 template<typename GV, class MODEL, int dofel>
14 class linearelasticity :
15 public NumericalJacobianApplyVolume<linearelasticity<GV,MODEL,dofel>>,
16 public NumericalJacobianApplyBoundary<linearelasticity<GV,MODEL,dofel>>,
17 public FullVolumePattern,
18 public LocalOperatorDefaultFlags,
19 public InstationaryLocalOperatorDefaultMethods<double>,
20 public Dune::PDELab::NumericalJacobianVolume<linearelasticity<GV,MODEL,dofel> >,
21 public Dune::PDELab::NumericalJacobianBoundary<linearelasticity<GV,MODEL,dofel> >
22
23 {
24 public:
25 // pattern assembly flags
26 enum { doPatternVolume = true };
27
28 // residual assembly flags
29 enum { doAlphaVolume = true };
30 enum { doLambdaVolume = true };
31 enum { doLambdaBoundary = true };
32
33 //constructor
34 linearelasticity (const GV gv_, const MODEL& myModel_, int intorder_=2, double eps_=100) :
35 myModel(myModel_),gv(gv_),intorder(intorder_),eps(eps_)
36 {
37 my_rank = gv.comm().rank();
38 }
39
40 //Write the local stiffness matrix
41 template<typename DuneMatrix1,typename DuneMatrix2, typename DuneVector>
42 void BCB(const unsigned int nodes_per_element, DuneVector &gradphi, DuneMatrix1 &C, DuneMatrix2 &K) const
43 {
44 Dune::FieldMatrix<double,6,dofel> tmp(0.0), B(0.0);
45 for (unsigned int i = 0; i < nodes_per_element; i++)
46 {
47 B[0][i] = gradphi[i][0]; // E11
48 B[1][i + nodes_per_element] = gradphi[i][1]; // E22
49 B[2][i + 2 * nodes_per_element] = gradphi[i][2]; // E22
50 B[3][i + nodes_per_element] = gradphi[i][2]; B[3][i + 2 * nodes_per_element] = gradphi[i][1]; // E23
51 B[4][i] = gradphi[i][2]; B[4][i + 2 * nodes_per_element] = gradphi[i][0]; // E13
52 B[5][i] = gradphi[i][1]; B[5][i + nodes_per_element] = gradphi[i][0]; // E12
53 }
54
55 //tmp = C * B
56 mm<DuneMatrix1,Dune::FieldMatrix<double,6,dofel>,Dune::FieldMatrix<double,6,dofel>>(C,B,tmp);
57
58 // K = B' * tmp = B' * C * B
59 mtm<Dune::FieldMatrix<double,6,dofel>,Dune::FieldMatrix<double,6,dofel>,DuneMatrix2>(B,tmp,K);
60 } //end BCB
61
62 //Write the local stiffness matrix
63 template<typename DuneMatrix, typename DuneVector1, typename DuneVector2, typename DuneVector>
64 void Bdelta(const unsigned int nodes_per_element, DuneVector1 &gradphi, DuneMatrix &C, DuneVector2 &delta, DuneVector &result) const
65 {
66 Dune::FieldMatrix<double,6,dofel> B(0.0);
67 for (unsigned int i = 0; i < nodes_per_element; i++)
68 {
69 B[0][i] = gradphi[i][0]; // E11
70 B[1][i + nodes_per_element] = gradphi[i][1]; // E22
71 B[2][i + 2 * nodes_per_element] = gradphi[i][2]; // E22
72 B[3][i + nodes_per_element] = gradphi[i][2]; B[3][i + 2 * nodes_per_element] = gradphi[i][1]; // E23
73 B[4][i] = gradphi[i][2]; B[4][i + 2 * nodes_per_element] = gradphi[i][0]; // E13
74 B[5][i] = gradphi[i][1]; B[5][i + nodes_per_element] = gradphi[i][0]; // E12
75 }
76 DuneVector2 tmp(0.0);
77 C.mv(delta,tmp);
78 B.mtv(tmp,result);
79 } //end Bdelta
80
81 template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
82 void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, M & mat) const
83 {
84 // dimensions
85 const int dim = 3;
86
87 // extract local function spaces
88 typedef typename LFSU::template Child<0>::Type LFSU_U1;
89 typedef typename LFSU::template Child<1>::Type LFSU_U2;
90 typedef typename LFSU::template Child<2>::Type LFSU_U3;
91
92 const LFSU_U1& lfsu_u1 = lfsu.template child<0>();
93 const LFSU_U2& lfsu_u2 = lfsu.template child<1>();
94 const LFSU_U3& lfsu_u3 = lfsu.template child<2>();
95
96 const unsigned int npe = lfsu_u1.size();
97 assert(npe == dofel/dim); // ensure dof per element are set correctly
98
99 // domain and range field type
100 typedef typename LFSU_U1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::DomainFieldType DF;
101 typedef typename M::value_type RF;
102 typedef typename LFSU_U1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::JacobianType JacobianType;
103
104 // select quadrature rule
105 GeometryType gt = eg.geometry().type();
106 const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,myModel.intOrder);
107
108 // Evaluate elasticity tensor
109 int id = gv.indexSet().index(eg.entity());
110 Dune::FieldMatrix<double,6,6> C = myModel.getElasticTensor(id);
111
112
113 Dune::FieldMatrix<double,dofel,dofel> Kh(0.0);
114 for (typename Dune::QuadratureRule<DF,dim>::const_iterator it = rule.begin(),endit = rule.end(); it != endit; ++it){
115 //Rotate material tensor
116 auto theta = myModel.returnTheta(it->position());
117 //auto Ch = rotate(C, theta, 2); //rotate the stacking sequence first
118 auto Ch = myModel.TensorRotation(C, theta[1], 1); //then geometry rotation
119 //Ch = rotate(C, theta, 0);
120
121 Dune::FieldMatrix<double,dofel,dofel> K(0.0);
122 // Evaluate Jacobian
123 std::vector<JacobianType> js(npe);
124 lfsu_u1.finiteElement().localBasis().evaluateJacobian(it->position(),js);
125 // Transform gradient to real element
126 const typename EG::Geometry::JacobianInverseTransposed jac = eg.geometry().jacobianInverseTransposed(it->position());
127 std::vector<Dune::FieldVector<RF,dim> > gradphi(npe);
128 for (unsigned int i=0; i < npe; i++){
129 gradphi[i] = 0.0;
130 jac.umv(js[i][0],gradphi[i]);
131 }
132 Dune::FieldMatrix<double,6,dofel> B(0.0);
133 //Fill in stiffness matrix
134 BCB(npe,gradphi,Ch, K);
135
136 // Integration Point Weight
137 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
138 K *= factor;
139 Kh += K;
140 } // end for-each quadrature point
141 //if(intorder == 1) hourglass_control<Dune::FieldMatrix<double,dofel,dofel>>(Kh);
142
143 Dune::FieldMatrix<double,dim,dim> Kij(0.0);
144 for (unsigned int i = 0; i < npe; i++){
145 for (unsigned int j = 0; j < npe; j++){
146 Kij = {{Kh[i][j], Kh[i][npe + j], Kh[i][2*npe + j]},
147 {Kh[npe + i][j], Kh[npe + i][npe + j], Kh[npe + i][2*npe + j]},
148 {Kh[2*npe + i][j],Kh[2*npe + i][npe + j], Kh[2*npe + i][2*npe + j]}};
149 mat.accumulate(lfsu_u1,i,lfsu_u1,j,Kij[0][0] );
150 mat.accumulate(lfsu_u1,i,lfsu_u2,j,Kij[0][1] );
151 mat.accumulate(lfsu_u1,i,lfsu_u3,j,Kij[0][2] );
152 mat.accumulate(lfsu_u2,i,lfsu_u1,j,Kij[1][0] );
153 mat.accumulate(lfsu_u2,i,lfsu_u2,j,Kij[1][1]);
154 mat.accumulate(lfsu_u2,i,lfsu_u3,j,Kij[1][2] );
155 mat.accumulate(lfsu_u3,i,lfsu_u1,j,Kij[2][0] );
156 mat.accumulate(lfsu_u3,i,lfsu_u2,j,Kij[2][1] );
157 mat.accumulate(lfsu_u3,i,lfsu_u3,j,Kij[2][2] );
158 } // end for j
159 } // end for i
160 } // end jacobian_volume
161
162
163 // volume integral depending on test and ansatz functions
164 template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
165 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
166 {
167 // dimensions
168 const int dim = 3;
169
170 // extract local function spaces
171 typedef typename LFSU::template Child<0>::Type LFSU_U1;
172 typedef typename LFSU::template Child<1>::Type LFSU_U2;
173 typedef typename LFSU::template Child<2>::Type LFSU_U3;
174
175 const LFSU_U1& lfsu_u1 = lfsu.template child<0>();
176 const LFSU_U2& lfsu_u2 = lfsu.template child<1>();
177 const LFSU_U3& lfsu_u3 = lfsu.template child<2>();
178
179 const unsigned int npe = lfsu_u1.size();
180
181 // domain and range field type
182 typedef typename LFSU_U1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::DomainFieldType DF;
183 typedef typename R::value_type RF;
184 typedef typename LFSU_U1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::JacobianType JacobianType;
185
186 // select quadrature rule
187 GeometryType gt = eg.geometry().type();
188 const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,myModel.intOrder);
189
190 Dune::FieldMatrix<double,dofel,dofel> K(0.0);
191
192 // Unwrap u
193 FieldVector<double,dofel> u(0.0);
194 for (unsigned int i=0; i<npe; i++){
195 u[i] = x(lfsu_u1,i);
196 u[i + npe] = x(lfsu_u2,i);
197 u[i + 2*npe] = x(lfsu_u3,i);
198 }
199
200 // Evaluate elasticity tensor
201 int id = gv.indexSet().index(eg.entity());
202 Dune::FieldMatrix<double,6,6> C = myModel.getElasticTensor(id);
203
204 // Loop over quadrature points
205 for (typename Dune::QuadratureRule<DF,dim>::const_iterator it = rule.begin(),endit = rule.end(); it != endit; ++it)
206 {
207 //Rotate material tensor
208 auto theta = myModel.returnTheta(it->position());
209 //auto Ch = rotate(C, theta, 2); //rotate the stacking sequence first
210 auto Ch = myModel.TensorRotation(C, theta[1], 1); //then geometry rotation
211 //Ch = rotate(C, theta, 0);
212
213 // Evaluate Jacobian
214 std::vector<JacobianType> js(npe);
215 lfsu_u1.finiteElement().localBasis().evaluateJacobian(it->position(),js);
216
217 // Transform gradient to real element
218 const typename EG::Geometry::JacobianInverseTransposed jac = eg.geometry().jacobianInverseTransposed(it->position());
219 std::vector<Dune::FieldVector<RF,dim> > gradphi(npe);
220 for (unsigned int i=0; i < npe; i++){
221 gradphi[i] = 0.0;
222 jac.umv(js[i][0],gradphi[i]);
223 }
224
225 //Fill in stiffness matrix
226 BCB(npe,gradphi,Ch, K);
227
228 // if(intorder == 1) hourglass_control<Dune::FieldMatrix<double,dofel,dofel>>(K);
229
230 // Compute residual = K * U - checked
231 FieldVector<double,dofel> residual(0.0);
232
233
234 for (int i = 0; i < dofel; i++){
235 for (int j = 0; j < dofel; j++){
236 residual[i] += K[i][j] * u[j];
237 }
238 }
239
240 // geometric weight
241 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
242
243 for (unsigned int i=0; i < npe; i++){
244 r.accumulate(lfsu_u1,i,residual[i] * factor);
245 r.accumulate(lfsu_u2,i,residual[i + npe] * factor);
246 r.accumulate(lfsu_u3,i,residual[i + 2*npe] * factor);
247 }
248
249 } // end for each quadrature point
250
251 } // end alpha_volume
252
253 // volume integral depending only on test functions
254 template<typename EG, typename LFSV, typename R>
255 void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
256 {
257 // dimensions
258 const int dim = 3;
259
260 // extract local function spaces
261 typedef typename LFSV::template Child<0>::Type LFSV_U1;
262 typedef typename LFSV::template Child<1>::Type LFSV_U2;
263 typedef typename LFSV::template Child<2>::Type LFSV_U3;
264
265 const LFSV_U1& lfsv_u1 = lfsv.template child<0>();
266 const LFSV_U2& lfsv_u2 = lfsv.template child<1>();
267 const LFSV_U3& lfsv_u3 = lfsv.template child<2>();
268
269 const unsigned int npe = lfsv_u1.size();
270
271 // domain and range field type
272 typedef typename LFSV_U1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::DomainFieldType DF;
273 typedef typename R::value_type RF;
274 typedef typename LFSV_U1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::JacobianType JacobianType;
275 typedef typename LFSV_U1::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType RT_V;
276
277
278 int ele_id = gv.indexSet().index(eg.entity());
279
280 //Get input heat strains
281 Dune::FieldVector<double,6> delta(0.0);
282
283 myModel.evaluateHeat(delta,ele_id);
284
285 //Get input forcing term
286 Dune::FieldVector<double,3> f1(0.0);
287 myModel.evaluateWeight(f1,ele_id);
288
289 // select quadrature rule
290 GeometryType gt = eg.geometry().type();
291 const QuadratureRule<DF,dim>& rule = QuadratureRules<DF,dim>::rule(gt,myModel.intOrder);
292
293 // Evaluate elasticity tensor
294 int id = gv.indexSet().index(eg.entity());
295 Dune::FieldMatrix<double,6,6> C = myModel.getElasticTensor(id);
296
297 // Loop over quadrature points
298 for (typename Dune::QuadratureRule<DF,dim>::const_iterator it = rule.begin(),endit = rule.end(); it != endit; ++it)
299 {
300 //Rotate material tensor
301 auto theta = myModel.returnTheta(it->position());
302 //auto Ch = rotate(C, theta, 2); //rotate the stacking sequence first
303 auto Ch = myModel.TensorRotation(C, theta[1], 1); //then geometry rotation
304 //Ch = rotate(C, theta, 0);
305
306 //Shape functions
307 std::vector<RT_V> phi(npe);
308 lfsv_u1.finiteElement().localBasis().evaluateFunction(it->position(),phi);
309
310 // Evaluate Jacobian
311 std::vector<JacobianType> js(npe);
312 lfsv_u1.finiteElement().localBasis().evaluateJacobian(it->position(),js);
313
314 // Transform gradient to real element
315 const typename EG::Geometry::JacobianInverseTransposed jac = eg.geometry().jacobianInverseTransposed(it->position());
316 std::vector<Dune::FieldVector<RF,dim> > gradphi(npe);
317 for (unsigned int i=0; i < npe; i++){
318 gradphi[i] = 0.0;
319 jac.umv(js[i][0],gradphi[i]);
320 }
321
322 Dune::FieldVector<double,dofel> result;
323 Bdelta(npe,gradphi,Ch, delta, result);
324
325 // geometric weight
326 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
327
328 for (unsigned int i=0; i < npe; i++){
329 r.accumulate(lfsv_u1,i, -(f1[0]*phi[i] + result[i]) * factor);
330 r.accumulate(lfsv_u2,i, -(f1[1]*phi[i] + result[i+npe]) * factor);
331 r.accumulate(lfsv_u3,i, -(f1[2]*phi[i] + result[i+2*npe]) * factor);
332 }
333 } // end loop over quadrature points
334 } // end lambda_volume
335
336
337 // boundary integral
338 template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
339 void alpha_boundary (const IG& ig,
340 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
341 R& r_s) const
342 {
343 }
344
345 // jacobian of boundary term
346 template<typename IG, typename LFSV_HAT, typename R>
347 void lambda_boundary (const IG& ig, const LFSV_HAT& lfsv_hat, R& r) const
348 {
349 // dimensions
350 const int dim = 3;
351 const int dimw = 3;
352
353 // extract local function spaces
354 typedef typename LFSV_HAT::template Child<0>::Type LFSV_U1;
355 typedef typename LFSV_HAT::template Child<1>::Type LFSV_U2;
356 typedef typename LFSV_HAT::template Child<2>::Type LFSV_U3;
357
358 const LFSV_U1& lfsv_u1 = lfsv_hat.template child<0>();
359 const LFSV_U2& lfsv_u2 = lfsv_hat.template child<1>();
360 const LFSV_U3& lfsv_u3 = lfsv_hat.template child<2>();
361
362 const unsigned int npe = lfsv_u1.size();
363
364 // domain and range field type
365 typedef typename LFSV_U1::Traits::FiniteElementType::
366 Traits::LocalBasisType::Traits::DomainFieldType DF;
367 typedef typename R::value_type RF;
368 typedef typename LFSV_U1::Traits::FiniteElementType::
369 Traits::LocalBasisType::Traits::RangeType RT_V;
370
371 // select quadrature rule
372 GeometryType gtface = ig.geometry().type();
373 const QuadratureRule<DF,dim-1>& rule = QuadratureRules<DF,dim-1>::rule(gtface,intorder);
374
375 // loop over quadrature points and integrate normal flux
376 for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it = rule.begin(), endit = rule.end();it != endit;++it)
377 {
378 // position of quadrature point in local coordinates of element
379 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
380
381 // evaluate basis functions
382 std::vector<RT_V> phi(npe);
383 lfsv_u1.finiteElement().localBasis().evaluateFunction(local,phi);
384
385 // Compute Weight Factor
386 const RF factor = it->weight() * ig.geometry().integrationElement(it->position());
387 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(it->position());
388 // evaluate flux boundary condition - needs updating
389 Dune::FieldVector<DF,dimw> neumann_stress(0.0);
390 myModel.evaluateNeumann(ig.geometry().global(it->position()),neumann_stress,normal);
391 for (size_t i=0; i < npe; i++){
392 r.accumulate(lfsv_u1,i, neumann_stress[0] * phi[i] * factor);
393 r.accumulate(lfsv_u2,i, neumann_stress[1] * phi[i] * factor);
394 r.accumulate(lfsv_u3,i, neumann_stress[2] * phi[i] * factor);
395
396 }
397 }
398 }
399
400 private:
401 template<typename Matrix1, typename Matrix2, typename Matrix3>
402 void mm(Matrix1& A, Matrix2&B, Matrix3& C) const
403 {
404 int a1 = A.N(); int a2 = A.M();
405 int b1 = B.N(); int b2 = B.M();
406 assert(a2 == b1);
407 C = 0.0;
408 for (int i = 0; i < a1; i++){
409 for (int j = 0; j < b2; j++){
410 for (int k = 0; k < a2; k++){
411 C[i][j] += A[i][k] * B[k][j];
412 }
413 }
414 }
415 }
416
417 template<typename Matrix1, typename Matrix2, typename Matrix3>
418 void mtm(Matrix1& A, Matrix2&B, Matrix3& C) const
419 {
420
421 int a1 = A.N(); int a2 = A.M();
422 int b1 = B.N(); int b2 = B.M();
423 assert(a1 == b1);
424 C = 0.0;
425 for (int i = 0; i < a2; i++){
426 for (int j = 0; j < b2; j++){
427 for (int k = 0; k < a1; k++){
428 C[i][j] += A[k][i] * B[k][j];
429 }
430 }
431 }
432 }
433
434
435
436 template<typename Matrix>
437 void hourglass_control(Matrix & Kh) const{
438 if(dofel != 24){std::cout << "Warning: hourglass control not implemented for this case";}
439 Dune::FieldMatrix<double,dofel,dofel> v(0.0),vT(0.0);
440 v = {
441 { 1, 0, 0, 0, -1, 1, -1, 0, 0, -1, -1, 0, 1, 0, 0, 1, 0, 0, -1, 0, 0, -1, 0, 0},
442 { 1, 0, 0, 0, -1, 1, 1, 0, 0, -1, -1, 0, 1, 0, 0, -1, 0, 0, 1, 0, 0, 1, 0, 0},
443 { 1, 0, 0, 0, -1, -1, -1, 0, 0, 1, -1, 0, -1, 0, 0, -1, 0, 0, 1, 0, 0, -1, 0, 0},
444 { 1, 0, 0, 0, -1, -1, 1, 0, 0, 1, -1, 0, -1, 0, 0, 1, 0, 0, -1, 0, 0, 1, 0, 0},
445 { 1, 0, 0, 0, 1, 1, -1, 0, 0, -1, 1, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, 1, 0, 0},
446 { 1, 0, 0, 0, 1, 1, 1, 0, 0, -1, 1, 0, -1, 0, 0, 1, 0, 0, 1, 0, 0, -1, 0, 0},
447 { 1, 0, 0, 0, 1, -1, -1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0},
448 { 1, 0, 0, 0, 1, -1, 1, 0, 0, 1, 1, 0, 1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0},
449 { 0, 1, 0, 1, 0, -1, 0, -1, 0, -1, 0, -1, 0, 1, 0, 0, 1, 0, 0, -1, 0, 0, -1, 0},
450 { 0, 1, 0, 1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, 0, -1, 0, 0, 1, 0, 0, 1, 0},
451 { 0, 1, 0, 1, 0, -1, 0, 1, 0, -1, 0, -1, 0, -1, 0, 0, -1, 0, 0, 1, 0, 0, -1, 0},
452 { 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, -1, 0, -1, 0, 0, 1, 0, 0, -1, 0, 0, 1, 0},
453 { 0, 1, 0, -1, 0, -1, 0, -1, 0, -1, 0, 1, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, 1, 0},
454 { 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 0, 1, 0, 0, -1, 0},
455 { 0, 1, 0, -1, 0, -1, 0, 1, 0, -1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0},
456 { 0, 1, 0, -1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0},
457 { 0, 0, 1, -1, 1, 0, 0, 0, -1, 0, -1, -1, 0, 0, 1, 0, 0, 1, 0, 0, -1, 0, 0, -1},
458 { 0, 0, 1, -1, -1, 0, 0, 0, -1, 0, 1, -1, 0, 0, 1, 0, 0, -1, 0, 0, 1, 0, 0, 1},
459 { 0, 0, 1, 1, 1, 0, 0, 0, -1, 0, -1, 1, 0, 0, -1, 0, 0, -1, 0, 0, 1, 0, 0, -1},
460 { 0, 0, 1, 1, -1, 0, 0, 0, -1, 0, 1, 1, 0, 0, -1, 0, 0, 1, 0, 0, -1, 0, 0, 1},
461 { 0, 0, 1, -1, 1, 0, 0, 0, 1, 0, -1, -1, 0, 0, -1, 0, 0, -1, 0, 0, -1, 0, 0, 1},
462 { 0, 0, 1, -1, -1, 0, 0, 0, 1, 0, 1, -1, 0, 0, -1, 0, 0, 1, 0, 0, 1, 0, 0, -1},
463 { 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, -1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1},
464 { 0, 0, 1, 1, -1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, -1, 0, 0, -1, 0, 0, -1}
465 };
466 //Normalise
467 for(int m=0; m< dofel; m++){
468 double sum = 0;
469 for(int j=0; j< dofel; j++){
470 sum+=v[j][m]*v[j][m];
471 }
472 for(int j=0; j< dofel; j++){
473 v[j][m]*=sqrt(1./sum);
474 }
475
476 }
477 Composites::transpose<Dune::FieldMatrix<double,dofel,dofel>, dofel>(v,vT);
478 //Diagonalise matrix
479 Kh.rightmultiply(v);
480 Kh.leftmultiply(vT);
481 //Remove improper nullspace
482 for(int i = 12; i < dofel; i++){
483 Kh[i][i] = eps;
484 }
485 //Return to initial basis
486 Kh.rightmultiply(vT);
487 Kh.leftmultiply(v);
488 }
489
490 const MODEL& myModel;
491 const GV gv;
492 //const typename GV::IndexSet& is;
493 int intorder;
494 double eps;
495 int my_rank;
496 };
497
498 }
499}
500#endif
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Sep 5, 22:35, 2025)