MFEM  v4.6.0 Finite element discretization library
mtop_integrators.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
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 MTOPINTEGRATORS_HPP
13 #define MTOPINTEGRATORS_HPP
14
15 #include "mfem.hpp"
16 #include "paramnonlinearform.hpp"
17
18 #include <map>
19
20 namespace mfem
21 {
22
23 /// Base class for representing function at integration points.
25 {
26 public:
27  virtual ~BaseQFunction() {}
28
29  /// Returns a user defined string identifying the function.
30  virtual std::string GetType()=0;
31
32  // Returns the energy at an integration point.
33  virtual
35  mfem::Vector &dd, mfem::Vector &uu)
36  {
37  return 0.0;
38  }
39
40  // Returns the residual at an integration point.
41  virtual
43  mfem::Vector &dd, mfem::Vector &uu, mfem::Vector &rr)=0;
44
45  /// Returns the gradient of the residual at a integration point.
46  virtual
49
50  /// Returns the gradient of the residual with respect to the design
51  /// parameters, multiplied by the adjoint.
52  virtual
54  mfem::Vector &dd, mfem::Vector &uu,
55  mfem::Vector &aa, mfem::Vector &rr)=0;
56
57 };
58
59 /* QLinearDiffusion implements methods for computing the energy, the residual,
60  * gradient of the residual and the product of the adjoint fields with the
61  * derivative of the residual with respect to the parameters. All computations
62  * are performed at a integration point. Therefore the vectors (vv,uu,aa,rr ..)
63  * hold the fields' values and the fields' derivatives at the integration
64  * point. For example for a single scalar parametric field representing the
65  * density in topology optimization the vector dd will have size one and the
66  * element will be the density at the integration point. The map between state
67  * and parameter is not fixed and depends on the implementation of the QFunction
68  * class. */
70 {
71 public:
73  double pp=1.0, double minrho=1e-7, double betac=4.0, double etac=0.5):
74  diff(diffco),load(hsrco), powerc(pp), rhomin(minrho), beta(betac), eta(etac)
75  {
76
77  }
78
79  virtual std::string GetType() override
80  {
81  return "QLinearDiffusion";
82  }
83
84  virtual
86  Vector &dd, Vector &uu) override
87  {
88  // dd[0] - density
92  // uu[3] - temperature/scalar field
93
94  double di=diff.Eval(T,ip);
96  // Computes the physical density using projection.
97  double rz=0.5+0.5*std::tanh(beta*(dd[0]-eta)); //projection
98  // Computes the diffusion coefficient at the integration point.
99  double fd=di*(std::pow(rz,powerc)+rhomin);
100  // Computes the sum of the energy and the product of the temperature and
101  // the external input at the integration point.
102  double rez = 0.5*(uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2])*fd-uu[3]*ll;
103  return rez;
104  }
105
106  /// Returns the derivative of QEnergy with respect to the state vector uu.
107  virtual
109  Vector &dd, Vector &uu, Vector &rr) override
110  {
111  double di=diff.Eval(T,ip);
113  double rz=0.5+0.5*std::tanh(beta*(dd[0]-eta));
114  double fd=di*(std::pow(rz,powerc)+rhomin);
115
116  rr[0]=uu[0]*fd;
117  rr[1]=uu[1]*fd;
118  rr[2]=uu[2]*fd;
119  rr[3]=-ll;
120  }
121
122
123  // Returns the derivative, with respect to the density, of the product of
124  // the adjoint field with the residual at the integration point ip.
125  virtual
127  Vector &dd, Vector &uu, Vector &aa, Vector &rr) override
128  {
129  double di=diff.Eval(T,ip);
130  double tt=std::tanh(beta*(dd[0]-eta));
131  double rz=0.5+0.5*tt;
132  double fd=di*powerc*std::pow(rz,powerc-1.0)*0.5*(1.0-tt*tt)*beta;
133
134  rr[0] = -(aa[0]*uu[0]+aa[1]*uu[1]+aa[2]*uu[2])*fd;
135  }
136
137  // Returns the gradient of the residual with respect to the state vector at
138  // the integration point ip.
139  virtual
141  Vector &dd, Vector &uu, DenseMatrix &hh) override
142  {
143  double di=diff.Eval(T,ip);
144  double tt=std::tanh(beta*(dd[0]-eta));
145  double rz=0.5+0.5*tt;
146  double fd=di*(std::pow(rz,powerc)+rhomin);
147  hh=0.0;
148
149  hh(0,0)=fd;
150  hh(1,1)=fd;
151  hh(2,2)=fd;
152  hh(3,3)=0.0;
153  }
154
155 private:
156  mfem::Coefficient& diff; //diffusion coefficient
158  double powerc; //penalization coefficient
159  double rhomin; //lower bound for the density
160  double beta; //controls the sharpness of the projection
161  double eta; //projection threshold for tanh
162 };
163
164 /// Provides implementation of an integrator for linear diffusion with
165 /// parametrization provided by a density field. The setup is standard for
166 /// topology optimization problems.
168 {
169 public:
171  {
172
173  }
174
175  /// Computes the local energy.
176  virtual
178  const Array<const FiniteElement *> &pel,
180  const Array<const Vector *> &elfun,
181  const Array<const Vector *> &pelfun) override;
182
183  /// Computes the element's residual.
184  virtual
186  const Array<const FiniteElement *> &pel,
188  const Array<const Vector *> &elfun,
189  const Array<const Vector *> &pelfun,
190  const Array<Vector *> &elvec) override;
191
192  /// Computes the stiffness/tangent matrix.
193  virtual
195  const Array<const FiniteElement *> &pel,
197  const Array<const Vector *> &elfun,
198  const Array<const Vector *> &pelfun,
199  const Array2D<DenseMatrix *> &elmats) override;
200
201  /// Computes the product of the adjoint solution and the derivative of the
202  /// residual with respect to the parametric fields.
203  virtual
205  const Array<const FiniteElement *> &pel,
207  const Array<const Vector *> &elfun,
208  const Array<const Vector *> &alfun,
209  const Array<const Vector *> &pelfun,
210  const Array<Vector *> &elvec) override;
211 private:
212  BaseQFunction& qfun;
213 };
214
215
216 /// Computes an example of nonlinear objective
217 /// \f$\int \rm{field}*\rm{field}*\rm{weight})\rm{d}\Omega_e\f$.
219 {
220 public:
221
223  {
224
225  }
226
227  /// Returns the objective contribution at element level.
228  virtual
231  const Array<const Vector *> &elfun) override;
232
233  /// Returns the gradient of the objective contribution at element level.
234  virtual
237  const Array<const Vector *> &elfun,
238  const Array<Vector *> &elvec) override;
239 };
240
241 }
242
243 #endif
virtual void QResidual(ElementTransformation &T, const IntegrationPoint &ip, mfem::Vector &dd, mfem::Vector &uu, mfem::Vector &rr)=0
virtual void AssembleElementVector(const Array< const FiniteElement *> &el, ElementTransformation &Tr, const Array< const Vector *> &elfun, const Array< Vector *> &elvec) override
Returns the gradient of the objective contribution at element level.
virtual double GetElementEnergy(const Array< const FiniteElement *> &el, const Array< const FiniteElement *> &pel, ElementTransformation &Tr, const Array< const Vector *> &elfun, const Array< const Vector *> &pelfun) override
Computes the local energy.
virtual void QResidual(ElementTransformation &T, const IntegrationPoint &ip, Vector &dd, Vector &uu, Vector &rr) override
Returns the derivative of QEnergy with respect to the state vector uu.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
ParametricLinearDiffusion(BaseQFunction &qfunm)
virtual void AQResidual(ElementTransformation &T, const IntegrationPoint &ip, mfem::Vector &dd, mfem::Vector &uu, mfem::Vector &aa, mfem::Vector &rr)=0
virtual double QEnergy(ElementTransformation &T, const IntegrationPoint &ip, mfem::Vector &dd, mfem::Vector &uu)
Base class for representing function at integration points.
virtual std::string GetType() override
Returns a user defined string identifying the function.
virtual void AQResidual(ElementTransformation &T, const IntegrationPoint &ip, Vector &dd, Vector &uu, Vector &aa, Vector &rr) override
virtual double QEnergy(ElementTransformation &T, const IntegrationPoint &ip, Vector &dd, Vector &uu) override
virtual void AssembleElementVector(const Array< const FiniteElement *> &el, const Array< const FiniteElement *> &pel, ElementTransformation &Tr, const Array< const Vector *> &elfun, const Array< const Vector *> &pelfun, const Array< Vector *> &elvec) override
Computes the element&#39;s residual.
virtual void QGradResidual(ElementTransformation &T, const IntegrationPoint &ip, Vector &dd, Vector &uu, DenseMatrix &hh) override
Returns the gradient of the residual at a integration point.
Dynamic 2D array using row-major layout.
Definition: array.hpp:354
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
virtual void QGradResidual(ElementTransformation &T, const IntegrationPoint &ip, mfem::Vector &dd, mfem::Vector &uu, mfem::DenseMatrix &hh)=0
Returns the gradient of the residual at a integration point.
QLinearDiffusion(mfem::Coefficient &diffco, mfem::Coefficient &hsrco, double pp=1.0, double minrho=1e-7, double betac=4.0, double etac=0.5)
Class for integration point with weight.
Definition: intrules.hpp:31
virtual double GetElementEnergy(const Array< const FiniteElement *> &el, ElementTransformation &Tr, const Array< const Vector *> &elfun) override
Returns the objective contribution at element level.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
Vector data type.
Definition: vector.hpp:58
virtual std::string GetType()=0
Returns a user defined string identifying the function.
virtual void AssembleElementGrad(const Array< const FiniteElement *> &el, const Array< const FiniteElement *> &pel, ElementTransformation &Tr, const Array< const Vector *> &elfun, const Array< const Vector *> &pelfun, const Array2D< DenseMatrix *> &elmats) override
Computes the stiffness/tangent matrix.
virtual void AssemblePrmElementVector(const Array< const FiniteElement *> &el, const Array< const FiniteElement *> &pel, ElementTransformation &Tr, const Array< const Vector *> &elfun, const Array< const Vector *> &alfun, const Array< const Vector *> &pelfun, const Array< Vector *> &elvec) override