MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
example.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#ifndef ADEXAMPLE_HPP
13#define ADEXAMPLE_HPP
14
15#include "mfem.hpp"
16#include "admfem.hpp"
17#include <memory>
18#include <iostream>
19#include <fstream>
20
21namespace mfem
22{
23
24/// Example: Implementation of the residual evaluation for p-Laplacian
25/// problem. The residual is evaluated at the integration points for PDE
26/// parameters vparam and state fields (derivatives with respect to x,y,z and
27/// primal field) stored in vector uu.
28template<typename TDataType, typename TParamVector, typename TStateVector,
29 int residual_size, int state_size, int param_size>
31{
32public:
33 /// The operator returns the first derivative of the energy with respect to
34 /// all state variables. These are set in vector uu and consist of the
35 /// derivatives with respect to x,y,z and the primal field. The derivative is
36 /// stored in vector rr with length equal to the length of vector uu.
37 void operator()(TParamVector &vparam, TStateVector &uu, TStateVector &rr)
38 {
39 MFEM_ASSERT(residual_size==4,
40 "PLaplacianResidual residual_size should be equal to 4!");
41 real_t pp = vparam[0];
42 real_t ee = vparam[1];
43 real_t ff = vparam[2];
44
45 // The vector rr holds the gradients of the following expression:
46 // (u_x^2+u_y^2+u_z^2+\varepsilon^2)^(p/2)-f.u,
47 // where u_x,u_y,u_z are the gradients of the scalar field u.
48 // The state vector is defined as uu=[u_x,u_y,u_z,u].
49
50 TDataType norm2 = uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2];
51 TDataType tvar = pow(ee * ee + norm2, (pp - 2.0) / 2.0);
52
53 rr[0] = tvar * uu[0];
54 rr[1] = tvar * uu[1];
55 rr[2] = tvar * uu[2];
56 rr[3] = -ff;
57 }
58};
59
60/// Defines template class (functor) for evaluating the energy of the
61/// p-Laplacian problem. The input parameters vparam are: vparam[0] - the
62/// p-Laplacian power, vparam[1] small value ensuring exciting of an unique
63/// solution, and vparam[2] - the distributed external input to the PDE. The
64/// template parameter TDataType will be replaced by the compiler with the
65/// appropriate AD type for automatic differentiation. The TParamVector
66/// represents the vector type used for the parameter vector, and TStateVector
67/// the vector type used for the state vector. The template parameters
68/// state_size and param_size provide information for the size of the state and
69/// the parameters vectors.
70template<typename TDataType, typename TParamVector, typename TStateVector
71 , int state_size, int param_size>
73{
74public:
75 /// Returns the energy of a p-Laplacian for state field input provided in
76 /// vector uu and parameters provided in vector vparam.
77 TDataType operator()(TParamVector &vparam, TStateVector &uu)
78 {
79 MFEM_ASSERT(state_size==4,"MyEnergyFunctor state_size should be equal to 4!");
80 MFEM_ASSERT(param_size==3,"MyEnergyFunctor param_size should be equal to 3!");
81 real_t pp = vparam[0];
82 real_t ee = vparam[1];
83 real_t ff = vparam[2];
84
85 TDataType u = uu[3];
86 TDataType norm2 = uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2];
87
88 TDataType rez = pow(ee * ee + norm2, pp / 2.0) / pp - ff * u;
89 return rez;
90 }
91};
92
93
94/// Implements integrator for a p-Laplacian problem. The integrator is based on
95/// a class QFunction utilized for evaluating the energy, the first derivative
96/// (residual) and the Hessian of the energy (the Jacobian of the residual).
97/// The template parameter CQVectAutoDiff represents the automatically
98/// differentiated energy or residual implemented by the user.
99/// CQVectAutoDiff::VectorFunc(Vector parameters, Vector state,Vector residual)
100/// evaluates the residual at an integration point.
101/// CQVectAutoDiff::Jacobian(Vector parameters, Vector state, Matrix hessian)
102/// evaluates the Hessian of the energy(the Jacobian of the residual).
103template<class CQVectAutoDiff>
105{
106protected:
110
111 CQVectAutoDiff rdf;
112
113public:
115 {
116 coeff = nullptr;
117 pp = nullptr;
118 load = nullptr;
119
120 vparam.SetSize(3);
121 vparam[0] = 2.0; // default power
122 vparam[1] = 1e-8; // default epsilon
123 vparam[2] = 1.0; // default load
124 }
125
126 pLaplaceAD(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr)
127 {
128 vparam.SetSize(3);
129 vparam[0] = 2.0; // default power
130 vparam[1] = 1e-8; // default epsilon
131 vparam[2] = 1.0; // default load
132
133 }
134
136 : pp(&pp_), coeff(&q), load(&ld_)
137 {
138 vparam.SetSize(3);
139 vparam[0] = 2.0; // default power
140 vparam[1] = 1e-8; // default epsilon
141 vparam[2] = 1.0; // default load
142 }
143
144 virtual ~pLaplaceAD() {}
145
148 const Vector &elfun)
149 {
150 real_t energy = 0.0;
151 const int ndof = el.GetDof();
152 const int ndim = el.GetDim();
153 const int spaceDim = trans.GetSpaceDim();
154 bool square = (ndim == spaceDim);
155 int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
156 const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
157
158 Vector shapef(ndof);
159 // derivatives in isoparametric coordinates
160 DenseMatrix dshape_iso(ndof, ndim);
161 // derivatives in physical space
162 DenseMatrix dshape_xyz(ndof, spaceDim);
163 Vector grad(spaceDim);
164
165 Vector uu(4); //[diff_x,diff_y,diff_z,u]
166
167 uu = 0.0;
168
169 // Calculates the functional/energy at an integration point.
171
172 real_t w;
173 real_t detJ;
174
175 for (int i = 0; i < ir.GetNPoints(); i++)
176 {
177 const IntegrationPoint &ip = ir.IntPoint(i);
178 trans.SetIntPoint(&ip);
179 w = trans.Weight();
180 detJ = (square ? w : w * w);
181 w = ip.weight * w;
182
183 el.CalcDShape(ip, dshape_iso);
184 el.CalcShape(ip, shapef);
185 // AdjugateJacobian = / adj(J), if J is square
186 // \ adj(J^t.J).J^t, otherwise
187 Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
188 // dshape_xyz should be divided by detJ for obtaining the real value
189 // calculate the gradient
190 dshape_xyz.MultTranspose(elfun, grad);
191
192 // set the power
193 if (pp != nullptr)
194 {
195 vparam[0] = pp->Eval(trans, ip);
196 }
197
198 // set the coefficient ensuring positiveness of the tangent matrix
199 if (coeff != nullptr)
200 {
201 vparam[1] = coeff->Eval(trans, ip);
202 }
203 // add the contribution from the load
204 if (load != nullptr)
205 {
206 vparam[2] = load->Eval(trans, ip);
207 }
208 // fill the values of vector uu
209 for (int jj = 0; jj < spaceDim; jj++)
210 {
211 uu[jj] = grad[jj] / detJ;
212 }
213 uu[3] = shapef * elfun;
214 // the energy is taken directly from the templated function
215 energy = energy + w * qfunc(vparam,uu);
216 }
217 return energy;
218 }
219
220 virtual void AssembleElementVector(const FiniteElement &el,
222 const Vector &elfun,
223 Vector &elvect)
224 {
225 MFEM_PERF_BEGIN("AssembleElementVector");
226 const int ndof = el.GetDof();
227 const int ndim = el.GetDim();
228 const int spaceDim = trans.GetSpaceDim();
229 int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
230 const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
231
232 Vector shapef(ndof);
233 DenseMatrix dshape_iso(ndof, ndim);
234 DenseMatrix dshape_xyz(ndof, spaceDim);
235 Vector lvec(ndof);
236 elvect.SetSize(ndof);
237 elvect = 0.0;
238
239 DenseMatrix B(ndof, 4); // [diff_x,diff_y,diff_z, shape]
240 Vector uu(4); // [diff_x,diff_y,diff_z,u]
241 Vector du(4);
242 B = 0.0;
243 uu = 0.0;
244 real_t w;
245
246 for (int i = 0; i < ir.GetNPoints(); i++)
247 {
248 const IntegrationPoint &ip = ir.IntPoint(i);
249 trans.SetIntPoint(&ip);
250 w = trans.Weight();
251 w = ip.weight * w;
252
253 el.CalcDShape(ip, dshape_iso);
254 el.CalcShape(ip, shapef);
255 Mult(dshape_iso, trans.InverseJacobian(), dshape_xyz);
256
257 // set the matrix B
258 for (int jj = 0; jj < spaceDim; jj++)
259 {
260 B.SetCol(jj, dshape_xyz.GetColumn(jj));
261 }
262 B.SetCol(3, shapef);
263
264 // set the power
265 if (pp != nullptr)
266 {
267 vparam[0] = pp->Eval(trans, ip);
268 }
269 // set the coefficient ensuring positiveness of the tangent matrix
270 if (coeff != nullptr)
271 {
272 vparam[1] = coeff->Eval(trans, ip);
273 }
274 // add the contribution from the load
275 if (load != nullptr)
276 {
277 vparam[2] = load->Eval(trans, ip);
278 }
279
280 // calculate uu
281 B.MultTranspose(elfun, uu);
282 // calculate derivative of the energy with respect to uu
283 rdf.VectorFunc(vparam,uu,du);
284 B.Mult(du, lvec);
285 elvect.Add(w, lvec);
286 } // end integration loop
287 MFEM_PERF_END("AssembleElementVector");
288 }
289
290 virtual void AssembleElementGrad(const FiniteElement &el,
292 const Vector &elfun,
293 DenseMatrix &elmat)
294 {
295 MFEM_PERF_BEGIN("AssembleElementGrad");
296 const int ndof = el.GetDof();
297 const int ndim = el.GetDim();
298 const int spaceDim = trans.GetSpaceDim();
299 int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
300 const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
301
302 Vector shapef(ndof);
303 DenseMatrix dshape_iso(ndof, ndim);
304 DenseMatrix dshape_xyz(ndof, spaceDim);
305 elmat.SetSize(ndof, ndof);
306 elmat = 0.0;
307
308 DenseMatrix B(ndof, 4); // [diff_x,diff_y,diff_z, shape]
309 DenseMatrix A(ndof, 4);
310 Vector uu(4); // [diff_x,diff_y,diff_z,u]
311 DenseMatrix duu(4, 4);
312 B = 0.0;
313 uu = 0.0;
314
315 real_t w;
316
317 for (int i = 0; i < ir.GetNPoints(); i++)
318 {
319 const IntegrationPoint &ip = ir.IntPoint(i);
320 trans.SetIntPoint(&ip);
321 w = trans.Weight();
322 w = ip.weight * w;
323
324 el.CalcDShape(ip, dshape_iso);
325 el.CalcShape(ip, shapef);
326 Mult(dshape_iso, trans.InverseJacobian(), dshape_xyz);
327
328 // set the matrix B
329 for (int jj = 0; jj < spaceDim; jj++)
330 {
331 B.SetCol(jj, dshape_xyz.GetColumn(jj));
332 }
333 B.SetCol(3, shapef);
334
335 // set the power
336 if (pp != nullptr)
337 {
338 vparam[0] = pp->Eval(trans, ip);
339 }
340 // set the coefficient ensuring positiveness of the tangent matrix
341 if (coeff != nullptr)
342 {
343 vparam[1] = coeff->Eval(trans, ip);
344 }
345 // add the contribution from the load
346 if (load != nullptr)
347 {
348 vparam[2] = load->Eval(trans, ip);
349 }
350
351 // calculate uu
352 B.MultTranspose(elfun, uu);
353 // calculate derivative of the energy with respect to uu
354 rdf.Jacobian(vparam,uu,duu);
355 Mult(B, duu, A);
356 AddMult_a_ABt(w, A, B, elmat);
357
358 } // end integration loop
359 MFEM_PERF_END("AssembleElementGrad");
360 }
361
362private:
363 Vector vparam; // [power, epsilon, load]
364
365};
366
367/// Implements hand-coded integrator for a p-Laplacian problem. Utilized as
368/// alternative for the pLaplaceAD class based on automatic differentiation.
370{
371protected:
375
376public:
378 {
379 coeff = nullptr;
380 pp = nullptr;
381 load = nullptr;
382 }
383
384 pLaplace(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr) {}
385
387 : pp(&pp_), coeff(&q), load(&ld_)
388 {}
389
390 virtual ~pLaplace() {}
391
394 const Vector &elfun)
395 {
396 real_t energy = 0.0;
397 const int ndof = el.GetDof();
398 const int ndim = el.GetDim();
399 const int spaceDim = trans.GetSpaceDim();
400 bool square = (ndim == spaceDim);
401 int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
402 const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
403
404 Vector shapef(ndof);
405 DenseMatrix dshape_iso(ndof, ndim);
406 DenseMatrix dshape_xyz(ndof, spaceDim);
407 Vector grad(spaceDim);
408
409 real_t w;
410 real_t detJ;
411 real_t nrgrad2;
412 real_t ppp = 2.0;
413 real_t eee = 0.0;
414
415 for (int i = 0; i < ir.GetNPoints(); i++)
416 {
417 const IntegrationPoint &ip = ir.IntPoint(i);
418 trans.SetIntPoint(&ip);
419 w = trans.Weight();
420 detJ = (square ? w : w * w);
421 w = ip.weight * w;
422
423 el.CalcDShape(ip, dshape_iso);
424 el.CalcShape(ip, shapef);
425 // AdjugateJacobian = / adj(J), if J is square
426 // \ adj(J^t.J).J^t, otherwise
427 Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
428 // dshape_xyz should be divided by detJ for obtaining the real value
429 // calculate the gradient
430 dshape_xyz.MultTranspose(elfun, grad);
431 nrgrad2 = grad * grad / (detJ * detJ);
432
433 // set the power
434 if (pp != nullptr)
435 {
436 ppp = pp->Eval(trans, ip);
437 }
438
439 // set the coefficient ensuring positiveness of the tangent matrix
440 if (coeff != nullptr)
441 {
442 eee = coeff->Eval(trans, ip);
443 }
444
445 energy = energy + w * std::pow(nrgrad2 + eee * eee, ppp / 2.0) / ppp;
446
447 // add the contribution from the load
448 if (load != nullptr)
449 {
450 energy = energy - w * (shapef * elfun) * load->Eval(trans, ip);
451 }
452 }
453 return energy;
454 }
455
456 virtual void AssembleElementVector(const FiniteElement &el,
458 const Vector &elfun,
459 Vector &elvect)
460 {
461 MFEM_PERF_BEGIN("AssembleElementVector");
462 const int ndof = el.GetDof();
463 const int ndim = el.GetDim();
464 const int spaceDim = trans.GetSpaceDim();
465 bool square = (ndim == spaceDim);
466 int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
467 const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
468
469 Vector shapef(ndof);
470 DenseMatrix dshape_iso(ndof, ndim);
471 DenseMatrix dshape_xyz(ndof, spaceDim);
472 Vector grad(spaceDim);
473 Vector lvec(ndof);
474 elvect.SetSize(ndof);
475 elvect = 0.0;
476
477 real_t w;
478 real_t detJ;
479 real_t nrgrad;
480 real_t aa;
481 real_t ppp = 2.0;
482 real_t eee = 0.0;
483
484 for (int i = 0; i < ir.GetNPoints(); i++)
485 {
486 const IntegrationPoint &ip = ir.IntPoint(i);
487 trans.SetIntPoint(&ip);
488 w = trans.Weight();
489 detJ = (square ? w : w * w);
490 w = ip.weight * w;
491
492 el.CalcDShape(ip, dshape_iso);
493 el.CalcShape(ip, shapef);
494 // AdjugateJacobian = / adj(J), if J is square
495 // \ adj(J^t.J).J^t, otherwise
496 Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
497 // dshape_xyz should be divided by detJ for obtaining the real value
498
499 // calculate the gradient
500 dshape_xyz.MultTranspose(elfun, grad);
501 nrgrad = grad.Norml2() / detJ;
502 // grad is not scaled so far, i.e., grad=grad/detJ
503
504 // set the power
505 if (pp != nullptr)
506 {
507 ppp = pp->Eval(trans, ip);
508 }
509
510 // set the coefficient ensuring positiveness of the tangent matrix
511 if (coeff != nullptr)
512 {
513 eee = coeff->Eval(trans, ip);
514 }
515 // compute (norm of the gradient)^2 + epsilon^2
516 aa = nrgrad * nrgrad + eee * eee;
517 aa = std::pow(aa, (ppp - 2.0) / 2.0);
518 dshape_xyz.Mult(grad, lvec);
519 elvect.Add(w * aa / (detJ * detJ), lvec);
520
521 // add loading
522 if (load != nullptr)
523 {
524 elvect.Add(-w * load->Eval(trans, ip), shapef);
525 }
526 } // end integration loop
527 MFEM_PERF_END("AssembleElementVector");
528 }
529
530 virtual void AssembleElementGrad(const FiniteElement &el,
532 const Vector &elfun,
533 DenseMatrix &elmat)
534 {
535 MFEM_PERF_BEGIN("AssembleElementGrad");
536 const int ndof = el.GetDof();
537 const int ndim = el.GetDim();
538 const int spaceDim = trans.GetSpaceDim();
539 bool square = (ndim == spaceDim);
540 int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
541 const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
542
543 DenseMatrix dshape_iso(ndof, ndim);
544 DenseMatrix dshape_xyz(ndof, spaceDim);
545 Vector grad(spaceDim);
546 Vector lvec(ndof);
547 // set the size of the element matrix
548 elmat.SetSize(ndof, ndof);
549 elmat = 0.0;
550
551 real_t w; // integration weight
552 real_t detJ;
553 real_t nrgrad; // norm of the gradient
554 real_t aa0; // original nonlinear diffusion coefficient
555 real_t aa1; // gradient of the above
556 real_t ppp = 2.0; // power in the P-Laplacian
557 real_t eee = 0.0; // regularization parameter
558
559 for (int i = 0; i < ir.GetNPoints(); i++)
560 {
561 const IntegrationPoint &ip = ir.IntPoint(i);
562 trans.SetIntPoint(&ip);
563 w = trans.Weight();
564 detJ = (square ? w : w * w);
565 w = ip.weight * w;
566
567 el.CalcDShape(ip, dshape_iso);
568 // AdjugateJacobian = / adj(J), if J is square
569 // \ adj(J^t.J).J^t, otherwise
570 Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
571 // dshape_xyz should be divided by detJ for obtaining the real value
572 // grad is not scaled so far,i.e., grad=grad/detJ
573
574 // set the power
575 if (pp != nullptr)
576 {
577 ppp = pp->Eval(trans, ip);
578 }
579 // set the coefficient ensuring positiveness of the tangent matrix
580 if (coeff != nullptr)
581 {
582 eee = coeff->Eval(trans, ip);
583 }
584
585 // calculate the gradient
586 dshape_xyz.MultTranspose(elfun, grad);
587 nrgrad = grad.Norml2() / detJ;
588 // (u_x^2+u_y^2+u_z^2+\varepsilon^2)
589 aa0 = nrgrad * nrgrad + eee * eee;
590 aa1 = std::pow(aa0, (ppp - 2.0) / 2.0);
591 aa0 = (ppp - 2.0) * std::pow(aa0, (ppp - 4.0) / 2.0);
592 dshape_xyz.Mult(grad, lvec);
593 w = w / (detJ * detJ);
594 AddMult_a_VVt(w * aa0 / (detJ * detJ), lvec, elmat);
595 AddMult_a_AAt(w * aa1, dshape_xyz, elmat);
596
597 } // end integration loop
598 MFEM_PERF_END("AssembleElementGrad");
599 }
600};
601
602/// Implements AD enabled integrator for a p-Laplacian problem. The tangent
603/// matrix is computed using the residual of the element. The template argument
604/// should be equal to the size of the residual vector (element vector), i.e.,
605/// the user should specify the size to match the exact vector size for the
606/// considered order of the shape functions.
607
608template<int sizeres=10>
610{
611protected:
615
616public:
618 {
619 coeff = nullptr;
620 pp = nullptr;
621 load = nullptr;
622 }
623
624 pLaplaceSL(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr) {}
625
627 : pp(&pp_), coeff(&q), load(&ld_)
628 {}
629
630 virtual ~pLaplaceSL() {}
631
634 const Vector &elfun)
635 {
636 real_t energy = 0.0;
637 const int ndof = el.GetDof();
638 const int ndim = el.GetDim();
639 const int spaceDim = trans.GetSpaceDim();
640 bool square = (ndim == spaceDim);
641 int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
642 const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
643
644 Vector shapef(ndof);
645 DenseMatrix dshape_iso(ndof, ndim);
646 DenseMatrix dshape_xyz(ndof, spaceDim);
647 Vector grad(spaceDim);
648
649 real_t w;
650 real_t detJ;
651 real_t nrgrad2;
652 real_t ppp = 2.0;
653 real_t eee = 0.0;
654
655 for (int i = 0; i < ir.GetNPoints(); i++)
656 {
657 const IntegrationPoint &ip = ir.IntPoint(i);
658 trans.SetIntPoint(&ip);
659 w = trans.Weight();
660 detJ = (square ? w : w * w);
661 w = ip.weight * w;
662
663 el.CalcDShape(ip, dshape_iso);
664 el.CalcShape(ip, shapef);
665 // AdjugateJacobian = / adj(J), if J is square
666 // \ adj(J^t.J).J^t, otherwise
667 Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
668 // dshape_xyz should be divided by detJ for obtaining the real value
669 // calculate the gradient
670 dshape_xyz.MultTranspose(elfun, grad);
671 nrgrad2 = grad * grad / (detJ * detJ);
672
673 // set the power
674 if (pp != nullptr)
675 {
676 ppp = pp->Eval(trans, ip);
677 }
678
679 // set the coefficient ensuring positiveness of the tangent matrix
680 if (coeff != nullptr)
681 {
682 eee = coeff->Eval(trans, ip);
683 }
684
685 energy = energy + w * std::pow(nrgrad2 + eee * eee, ppp / 2.0) / ppp;
686
687 // add the contribution from the load
688 if (load != nullptr)
689 {
690 energy = energy - w * (shapef * elfun) * load->Eval(trans, ip);
691 }
692 }
693 return energy;
694 }
695
696 virtual void AssembleElementVector(const FiniteElement &el,
698 const Vector &elfun,
699 Vector &elvect)
700 {
701 MFEM_PERF_BEGIN("AssembleElementVector");
702 const int ndof = el.GetDof();
703 const int ndim = el.GetDim();
704 const int spaceDim = trans.GetSpaceDim();
705 bool square = (ndim == spaceDim);
706 int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
707 const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
708
709 Vector shapef(ndof);
710 DenseMatrix dshape_iso(ndof, ndim);
711 DenseMatrix dshape_xyz(ndof, spaceDim);
712 Vector grad(spaceDim);
713 Vector lvec(ndof);
714 elvect.SetSize(ndof);
715 elvect = 0.0;
716
717 real_t w;
718 real_t detJ;
719 real_t nrgrad;
720 real_t aa;
721 real_t ppp = 2.0;
722 real_t eee = 0.0;
723
724 for (int i = 0; i < ir.GetNPoints(); i++)
725 {
726 const IntegrationPoint &ip = ir.IntPoint(i);
727 trans.SetIntPoint(&ip);
728 w = trans.Weight();
729 detJ = (square ? w : w * w);
730 w = ip.weight * w; //w;
731
732 el.CalcDShape(ip, dshape_iso);
733 el.CalcShape(ip, shapef);
734 // AdjugateJacobian = / adj(J), if J is square
735 // \ adj(J^t.J).J^t, otherwise
736 Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
737 // dshape_xyz should be divided by detJ for obtaining the real value
738
739 // calculate the gradient
740 dshape_xyz.MultTranspose(elfun, grad);
741 nrgrad = grad.Norml2() / detJ;
742 // grad is not scaled so far, i.e., grad=grad/detJ
743
744 // set the power
745 if (pp != nullptr)
746 {
747 ppp = pp->Eval(trans, ip);
748 }
749
750 // set the coefficient ensuring positiveness of the tangent matrix
751 if (coeff != nullptr)
752 {
753 eee = coeff->Eval(trans, ip);
754 }
755
756 aa = nrgrad * nrgrad + eee * eee;
757 aa = std::pow(aa, (ppp - 2.0) / 2.0);
758 dshape_xyz.Mult(grad, lvec);
759 elvect.Add(w * aa / (detJ * detJ), lvec);
760
761 // add loading
762 if (load != nullptr)
763 {
764 elvect.Add(-w * load->Eval(trans, ip), shapef);
765 }
766 } // end integration loop
767 MFEM_PERF_END("AssembleElementVector");
768 }
769
770 virtual void AssembleElementGrad(const FiniteElement &el,
772 const Vector &elfun,
773 DenseMatrix &elmat)
774 {
775 MFEM_PERF_BEGIN("AssembleElementGrad");
776 const int ndof = el.GetDof();
777 const int ndim = el.GetDim();
778 const int spaceDim = trans.GetSpaceDim();
779 bool square = (ndim == spaceDim);
780 int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
781 const IntegrationRule &ir(IntRules.Get(el.GetGeomType(), order));
782
783 DenseMatrix dshape_iso(ndof, ndim);
784 DenseMatrix dshape_xyz(ndof, spaceDim);
785 elmat.SetSize(ndof, ndof);
786 elmat = 0.0;
787
788 real_t w;
789 real_t detJ;
790 real_t ppp = 2.0;
791 real_t eee = 0.0;
792
793 mfem::Vector param(3); param=0.0;
794
795 // Computes the residual at an integration point. The implementation is a
796 // copy of the integration loop in AssembleElementVector.
797 auto resfun = [&](mfem::Vector& vparam, mfem::ad::ADVectorType& uu,
799 {
800
801 vres.SetSize(uu.Size()); vres=0.0;
802 mfem::ad::ADVectorType grad(spaceDim);
805 mfem::ad::ADVectorType lvec(ndof);
806
807 for (int q = 0; q < ir.GetNPoints(); q++)
808 {
809 lvec=0.0;
810
811 const IntegrationPoint &ip = ir.IntPoint(q);
812 trans.SetIntPoint(&ip);
813 w = trans.Weight();
814 detJ = (square ? w : w * w);
815 w = ip.weight * w;
816
817 el.CalcDShape(ip, dshape_iso);
818 // AdjugateJacobian = / adj(J), if J is square
819 // \ adj(J^t.J).J^t, otherwise
820 Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
821 // dshape_xyz should be divided by detJ for obtaining the real value
822 // grad is not scaled so far,i.e., grad=grad/detJ
823
824 // set the power
825 if (pp != nullptr)
826 {
827 ppp = pp->Eval(trans, ip);
828 }
829 // set the coefficient ensuring positiveness of the tangent matrix
830 if (coeff != nullptr)
831 {
832 eee = coeff->Eval(trans, ip);
833 }
834
835 grad=0.0;
836 // calculate the gradient
837 for (int i=0; i<spaceDim; i++)
838 {
839 for (int j=0; j<ndof; j++)
840 {
841 grad[i]= grad[i]+ dshape_xyz(j,i)*uu[j];
842 }
843 }
844
845
846 nrgrad= (grad*grad)/(detJ*detJ);
847
848 aa = nrgrad + eee * eee;
849 aa = pow(aa, (ppp - 2.0) / 2.0);
850
851 for (int i=0; i<spaceDim; i++)
852 {
853 for (int j=0; j<ndof; j++)
854 {
855 lvec[j] = lvec[j] + dshape_xyz(j,i) * grad[i];
856 }
857 }
858
859 for (int j=0; j<ndof; j++)
860 {
861 vres[j]=vres[j] + lvec[j] * (w*aa/(detJ*detJ));
862 }
863 }
864 };
865
866 mfem::Vector bla(elfun);
867 // calculate the gradient - only for a fixed ndof
869 fdr.Jacobian(param, bla, elmat);
870 MFEM_PERF_END("AssembleElementGrad");
871 }
872};
873
874} // namespace mfem
875
876#endif
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition densemat.cpp:220
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
Definition densemat.cpp:262
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
void SetCol(int c, const real_t *col)
void GetColumn(int c, Vector &col) const
Abstract class for all finite elements.
Definition fe_base.hpp:239
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition fe_base.hpp:333
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:316
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition fe_base.hpp:326
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:329
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
TDataType operator()(TParamVector &vparam, TStateVector &uu)
Definition example.hpp:77
void operator()(TParamVector &vparam, TStateVector &uu, TStateVector &rr)
Definition example.hpp:37
This class is used to express the local action of a general nonlinear finite element operator....
Templated vector data type.
Definition tadvector.hpp:41
void Jacobian(mfem::Vector &vparam, mfem::Vector &vstate, mfem::DenseMatrix &jac)
Definition admfem.hpp:70
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Definition vector.cpp:247
virtual real_t GetElementEnergy(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun)
Compute the local energy.
Definition example.hpp:146
Coefficient * pp
Definition example.hpp:107
virtual ~pLaplaceAD()
Definition example.hpp:144
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition example.hpp:290
pLaplaceAD(Coefficient &pp_)
Definition example.hpp:126
Coefficient * coeff
Definition example.hpp:108
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition example.hpp:220
CQVectAutoDiff rdf
Definition example.hpp:111
Coefficient * load
Definition example.hpp:109
pLaplaceAD(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
Definition example.hpp:135
Coefficient * pp
Definition example.hpp:612
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition example.hpp:770
pLaplaceSL(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
Definition example.hpp:626
virtual real_t GetElementEnergy(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun)
Compute the local energy.
Definition example.hpp:632
Coefficient * coeff
Definition example.hpp:613
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition example.hpp:696
virtual ~pLaplaceSL()
Definition example.hpp:630
pLaplaceSL(Coefficient &pp_)
Definition example.hpp:624
Coefficient * load
Definition example.hpp:614
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition example.hpp:530
virtual real_t GetElementEnergy(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun)
Compute the local energy.
Definition example.hpp:392
pLaplace(Coefficient &pp_)
Definition example.hpp:384
Coefficient * pp
Definition example.hpp:372
Coefficient * coeff
Definition example.hpp:373
pLaplace(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
Definition example.hpp:386
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition example.hpp:456
Coefficient * load
Definition example.hpp:374
virtual ~pLaplace()
Definition example.hpp:390
void trans(const Vector &u, Vector &x)
Definition ex27.cpp:412
codi::RealForward ADFloatType
Forward AD type declaration.
Definition admfem.hpp:30
void AddMult_a_ABt(real_t a, const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += a * A * B^t.
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:476
void AddMult_a_VVt(const real_t a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
void AddMult_a_AAt(real_t a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
float real_t
Definition config.hpp:43
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:486