MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mtop_integrators.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 #include "mtop_integrators.hpp"
13 
14 namespace mfem
15 {
16 
21  const Array<const Vector *> &elfun,
22  const Array<const Vector *> &pelfun)
23 {
24  int dof_u0 = el[0]->GetDof();
25  int dof_r0 = pel[0]->GetDof();
26 
27  int dim = el[0]->GetDim();
28  int spaceDim = Tr.GetSpaceDim();
29  if (dim != spaceDim)
30  {
31  mfem::mfem_error("ParametricLinearDiffusion::GetElementEnergy"
32  " is not defined on manifold meshes");
33  }
34 
35  // shape functions
36  Vector shu0(dof_u0);
37  Vector shr0(dof_r0);
38  DenseMatrix dsu0(dof_u0,dim);
39  DenseMatrix B(dof_u0, 4);
40  B=0.0;
41 
42  double w;
43 
44  Vector param(1); param=0.0;
45  Vector uu(4); uu=0.0;
46 
47  double energy =0.0;
48 
49  const IntegrationRule *ir;
50  {
51  int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
52  +pel[0]->GetOrder();
53  ir=&IntRules.Get(Tr.GetGeometryType(),order);
54  }
55 
56  for (int i = 0; i < ir->GetNPoints(); i++)
57  {
58  const IntegrationPoint &ip = ir->IntPoint(i);
59  Tr.SetIntPoint(&ip);
60  w=Tr.Weight();
61  w = ip.weight * w;
62 
63  el[0]->CalcPhysDShape(Tr,dsu0);
64  el[0]->CalcPhysShape(Tr,shu0);
65  pel[0]->CalcPhysShape(Tr,shr0);
66 
67  param[0]=shr0*(*pelfun[0]);
68 
69  // set the matrix B
70  for (int jj=0; jj<dim; jj++)
71  {
72  B.SetCol(jj,dsu0.GetColumn(jj));
73  }
74  B.SetCol(3,shu0);
75  B.MultTranspose(*elfun[0],uu);
76  energy=energy+w * qfun.QEnergy(Tr,ip,param,uu);
77  }
78  return energy;
79 }
80 
81 
86  const Array<const Vector *> &elfun,
87  const Array<const Vector *> &pelfun,
88  const Array<Vector *> &elvec)
89 {
90  int dof_u0 = el[0]->GetDof();
91  int dof_r0 = pel[0]->GetDof();
92 
93  int dim = el[0]->GetDim();
94 
95  elvec[0]->SetSize(dof_u0);
96  *elvec[0]=0.0;
97  int spaceDim = Tr.GetSpaceDim();
98  if (dim != spaceDim)
99  {
100  mfem::mfem_error("ParametricLinearDiffusion::AssembleElementVector"
101  " is not defined on manifold meshes");
102  }
103 
104  // shape functions
105  Vector shu0(dof_u0);
106  Vector shr0(dof_r0);
107  DenseMatrix dsu0(dof_u0,dim);
108  DenseMatrix B(dof_u0, 4);
109  B=0.0;
110 
111  double w;
112 
113  Vector param(1); param=0.0;
114  Vector uu(4); uu=0.0;
115  Vector rr(4);
116  Vector lvec; lvec.SetSize(dof_u0);
117 
118  const IntegrationRule *ir = nullptr;
119  int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
120  +pel[0]->GetOrder();
121  ir=&IntRules.Get(Tr.GetGeometryType(),order);
122 
123 
124  for (int i = 0; i < ir->GetNPoints(); i++)
125  {
126  const IntegrationPoint &ip = ir->IntPoint(i);
127  Tr.SetIntPoint(&ip);
128  w=Tr.Weight();
129  w = ip.weight * w;
130 
131  el[0]->CalcPhysDShape(Tr,dsu0);
132  el[0]->CalcPhysShape(Tr,shu0);
133  pel[0]->CalcPhysShape(Tr,shr0);
134 
135  param[0]=shr0*(*pelfun[0]);
136 
137  // set the matrix B
138  for (int jj=0; jj<dim; jj++)
139  {
140  B.SetCol(jj,dsu0.GetColumn(jj));
141  }
142  B.SetCol(3,shu0);
143  B.MultTranspose(*elfun[0],uu);
144  qfun.QResidual(Tr,ip,param, uu, rr);
145 
146  B.Mult(rr,lvec);
147  elvec[0]->Add(w,lvec);
148  }
149 }
150 
153  const Array<const FiniteElement *> &pel,
155  const Array<const Vector *> &elfun,
156  const Array<const Vector *> &pelfun,
157  const Array2D<DenseMatrix *> &elmats)
158 {
159  int dof_u0 = el[0]->GetDof();
160  int dof_r0 = pel[0]->GetDof();
161 
162  int dim = el[0]->GetDim();
163 
164  DenseMatrix* K=elmats(0,0);
165  K->SetSize(dof_u0,dof_u0);
166  (*K)=0.0;
167 
168  int spaceDim = Tr.GetSpaceDim();
169  if (dim != spaceDim)
170  {
171  mfem::mfem_error("ParametricLinearDiffusion::AssembleElementGrad"
172  " is not defined on manifold meshes");
173  }
174 
175  // shape functions
176  Vector shu0(dof_u0);
177  Vector shr0(dof_r0);
178  DenseMatrix dsu0(dof_u0,dim);
179  DenseMatrix B(dof_u0, 4);
180  DenseMatrix A(dof_u0, 4);
181  B=0.0;
182  double w;
183 
184  Vector param(1); param=0.0;
185  Vector uu(4); uu=0.0;
186  DenseMatrix hh(4,4);
187  Vector lvec; lvec.SetSize(dof_u0);
188 
189  const IntegrationRule *ir = nullptr;
190  int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
191  +pel[0]->GetOrder();
192  ir=&IntRules.Get(Tr.GetGeometryType(),order);
193 
194  for (int i = 0; i < ir->GetNPoints(); i++)
195  {
196  const IntegrationPoint &ip = ir->IntPoint(i);
197  Tr.SetIntPoint(&ip);
198  w = Tr.Weight();
199  w = ip.weight * w;
200 
201  el[0]->CalcPhysDShape(Tr,dsu0);
202  el[0]->CalcPhysShape(Tr,shu0);
203  pel[0]->CalcPhysShape(Tr,shr0);
204 
205  param[0]=shr0*(*pelfun[0]);
206 
207  // set the matrix B
208  for (int jj=0; jj<dim; jj++)
209  {
210  B.SetCol(jj,dsu0.GetColumn(jj));
211  }
212  B.SetCol(3,shu0);
213  B.MultTranspose(*elfun[0],uu);
214  qfun.QGradResidual(Tr,ip,param,uu,hh);
215  Mult(B,hh,A);
216  AddMult_a_ABt(w,A,B,*K);
217  }
218 }
219 
222  const Array<const FiniteElement *> &pel,
224  const Array<const Vector *> &elfun,
225  const Array<const Vector *> &alfun,
226  const Array<const Vector *> &pelfun,
227  const Array<Vector *> &elvec)
228 {
229  int dof_u0 = el[0]->GetDof();
230  int dof_r0 = pel[0]->GetDof();
231 
232  int dim = el[0]->GetDim();
233  Vector& e0 = *(elvec[0]);
234 
235  e0.SetSize(dof_r0);
236  e0=0.0;
237 
238  int spaceDim = Tr.GetSpaceDim();
239  if (dim != spaceDim)
240  {
241  mfem::mfem_error("ParametricLinearDiffusion::AssemblePrmElementVector"
242  " is not defined on manifold meshes");
243  }
244 
245  // shape functions
246  Vector shu0(dof_u0);
247  Vector shr0(dof_r0);
248  DenseMatrix dsu0(dof_u0,dim);
249  DenseMatrix B(dof_u0, 4);
250  B=0.0;
251 
252  double w;
253 
254  Vector param(1); param=0.0;
255  Vector uu(4); uu=0.0;
256  Vector aa(4); aa=0.0;
257  Vector rr(1);
258  Vector lvec0; lvec0.SetSize(dof_r0);
259 
260  const IntegrationRule *ir;
261  {
262  int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
263  +pel[0]->GetOrder();
264  ir=&IntRules.Get(Tr.GetGeometryType(),order);
265  }
266 
267  for (int i = 0; i < ir->GetNPoints(); i++)
268  {
269  const IntegrationPoint &ip = ir->IntPoint(i);
270  Tr.SetIntPoint(&ip);
271  w=Tr.Weight();
272  w = ip.weight * w;
273 
274  el[0]->CalcPhysDShape(Tr,dsu0);
275  el[0]->CalcPhysShape(Tr,shu0);
276  pel[0]->CalcPhysShape(Tr,shr0);
277 
278  param[0]=shr0*(*pelfun[0]);
279 
280  // set the matrix B
281  for (int jj=0; jj<dim; jj++)
282  {
283  B.SetCol(jj,dsu0.GetColumn(jj));
284  }
285  B.SetCol(3,shu0);
286  B.MultTranspose(*elfun[0],uu);
287  B.MultTranspose(*alfun[0],aa);
288 
289  qfun.AQResidual(Tr, ip, param, uu, aa, rr);
290 
291  lvec0=shr0;
292  lvec0*=rr[0];
293 
294  e0.Add(w,lvec0);
295  }
296 }
297 
301  const Array<const Vector *> &elfun)
302 {
303  int dof_u0 = el[0]->GetDof();
304  int dim = el[0]->GetDim();
305  int spaceDim = Tr.GetSpaceDim();
306  if (dim != spaceDim)
307  {
308  mfem::mfem_error("DiffusionObjIntegrator::GetElementEnergy"
309  " is not defined on manifold meshes");
310  }
311 
312  // shape functions
313  Vector shu0(dof_u0);
314 
315  double w;
316  double val;
317 
318  double energy = 0.0;
319 
320  const IntegrationRule *ir;
321  {
322  int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
323  ir=&IntRules.Get(Tr.GetGeometryType(),order);
324  }
325 
326  for (int i = 0; i < ir->GetNPoints(); i++)
327  {
328  const IntegrationPoint &ip = ir->IntPoint(i);
329  Tr.SetIntPoint(&ip);
330  w=Tr.Weight();
331 
332  w = ip.weight * w;
333 
334  el[0]->CalcPhysShape(Tr,shu0);
335 
336  val=shu0*(*elfun[0]);
337  energy=energy + w * val * val;
338  }
339  return 0.5*energy;
340 }
341 
345  const Array<const Vector *> &elfun,
346  const Array<Vector *> &elvec)
347 {
348  int dof_u0 = el[0]->GetDof();
349  int dim = el[0]->GetDim();
350  int spaceDim = Tr.GetSpaceDim();
351 
352  elvec[0]->SetSize(dof_u0);
353  *elvec[0]=0.0;
354 
355  if (dim != spaceDim)
356  {
357  mfem::mfem_error("DiffusionObjIntegrator::GetElementEnergy"
358  " is not defined on manifold meshes");
359  }
360 
361  // shape functions
362  Vector shu0(dof_u0);
363 
364  double w;
365  double val;
366 
367  const IntegrationRule *ir;
368  {
369  int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
370  ir=&IntRules.Get(Tr.GetGeometryType(),order);
371  }
372 
373  for (int i = 0; i < ir->GetNPoints(); i++)
374  {
375  const IntegrationPoint &ip = ir->IntPoint(i);
376  Tr.SetIntPoint(&ip);
377  w=Tr.Weight();
378 
379  w = ip.weight * w;
380 
381  el[0]->CalcPhysShape(Tr,shu0);
382 
383  val=shu0*(*elfun[0]);
384 
385  elvec[0]->Add(w*val,shu0);
386  }
387 }
388 
389 
390 } // end mfem namespace
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual int OrderGrad(const FiniteElement *fe) const =0
Return the order of .
virtual void QResidual(ElementTransformation &T, const IntegrationPoint &ip, mfem::Vector &dd, mfem::Vector &uu, mfem::Vector &rr)=0
void SetCol(int c, const double *col)
Definition: densemat.cpp:2154
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:162
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += a * A * B^t.
Definition: densemat.cpp:2940
virtual void AQResidual(ElementTransformation &T, const IntegrationPoint &ip, mfem::Vector &dd, mfem::Vector &uu, mfem::Vector &aa, mfem::Vector &rr)=0
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 double QEnergy(ElementTransformation &T, const IntegrationPoint &ip, mfem::Vector &dd, mfem::Vector &uu)
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:201
virtual double GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun) override
Returns the objective contribution at element level.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1292
Dynamic 2D array using row-major layout.
Definition: array.hpp:351
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
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.
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:252
Class for integration point with weight.
Definition: intrules.hpp:25
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 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
int dim
Definition: ex24.cpp:53
Vector data type.
Definition: vector.hpp:60
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
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.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
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.