MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
tmop_pa.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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 "../tmop.hpp"
13#include "../linearform.hpp"
14#include "../pgridfunc.hpp"
15#include "../tmop_tools.hpp"
19
20namespace mfem
21{
22
24 const FiniteElementSpace &fes)
25{
26 MFEM_VERIFY(PA.enabled, "PA extension setup has not been done!");
27 MFEM_VERIFY(PA.fes == &fes, "");
28 // TODO: we need a more robust way to check that the 'fes' used when
29 // AssemblePA() was called has not been modified or completely destroyed and
30 // a new object created at the same address.
31
32 // Form the Vector of node positions, depending on what's the input.
33 Vector xe(de.Size());
34 if (x_0)
35 {
36 // The input is the displacement.
37 add(PA.X0, de, xe);
38 }
39 else { xe = de; }
40
41 if (PA.Jtr_needs_update || targetC->UsesPhysicalCoordinates())
42 {
44 PA.Jtr_debug_grad = true;
45 }
46
47 if (PA.dim == 2)
48 {
51 }
52
53 if (PA.dim == 3)
54 {
57 }
58}
59
61{
62 const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
64 // Return immediately if limiting is not enabled
65 if (lim_coeff == nullptr) { return; }
66 MFEM_VERIFY(lim_nodes0, "internal error");
67
68 MFEM_VERIFY(PA.enabled, "AssemblePA_Limiting but PA is not enabled!");
69 MFEM_VERIFY(lim_func, "No TMOP_LimiterFunction specification!")
70 MFEM_VERIFY(dynamic_cast<TMOP_QuadraticLimiter*>(lim_func) ||
71 dynamic_cast<TMOP_ExponentialLimiter*>(lim_func),
72 "Only TMOP_QuadraticLimiter and TMOP_ExponentialLimiter are supported");
73
74 const FiniteElementSpace *fes = PA.fes;
75 const int NE = PA.ne;
76 if (NE == 0) { return; } // Quick return for empty processors
77 const IntegrationRule &ir = *PA.ir;
78
80
81 // H0 for lim_coeff, (dim x dim) Q-vector
82 PA.H0.UseDevice(true);
83 PA.H0.SetSize(PA.dim * PA.dim * PA.nq * NE, mt);
84
85 // lim_coeff -> PA.C0 (Q-vector)
86 PA.C0.UseDevice(true);
87 if (ConstantCoefficient* cQ =
88 dynamic_cast<ConstantCoefficient*>(lim_coeff))
89 {
90 PA.C0.SetSize(1, Device::GetMemoryType());
91 PA.C0.HostWrite();
92 PA.C0(0) = cQ->constant;
93 }
94 else
95 {
96 PA.C0.SetSize(PA.nq * PA.ne, Device::GetMemoryType());
97 auto C0 = Reshape(PA.C0.HostWrite(), PA.nq, PA.ne);
98 for (int e = 0; e < NE; ++e)
99 {
101 for (int q = 0; q < ir.GetNPoints(); ++q)
102 {
103 C0(q,e) = lim_coeff->Eval(T, ir.IntPoint(q));
104 }
105 }
106 }
107
108 // lim_nodes0 -> PA.XL (E-vector)
109 MFEM_VERIFY(lim_nodes0->FESpace() == fes, "");
110 const Operator *n0_R = fes->GetElementRestriction(ordering);
111 PA.XL.SetSize(n0_R->Height(), Device::GetMemoryType());
112 PA.XL.UseDevice(true);
113 n0_R->Mult(*lim_nodes0, PA.XL);
114
115 // Limiting distances: lim_dist -> PA.LD (E-vector)
116 // TODO: remove the hack for the case lim_dist == NULL.
117 const FiniteElementSpace *limfes = (lim_dist) ? lim_dist->FESpace() : fes;
118 const FiniteElement &lim_fe = *limfes->GetTypicalFE();
119 PA.maps_lim = &lim_fe.GetDofToQuad(ir, DofToQuad::TENSOR);
120 PA.LD.SetSize(NE*lim_fe.GetDof(), Device::GetMemoryType());
121 PA.LD.UseDevice(true);
122 if (lim_dist)
123 {
124 const Operator *ld_R = limfes->GetElementRestriction(ordering);
125 ld_R->Mult(*lim_dist, PA.LD);
126 }
127 else { PA.LD = 1.0; }
128}
129
131 const IntegrationRule &ir,
132 const Vector &xe,
133 DenseTensor &Jtr) const
134{
135 MFEM_VERIFY(Jtr.SizeI() == Jtr.SizeJ() && Jtr.SizeI() > 1, "");
136 const int dim = Jtr.SizeI();
137 bool done = false;
138 if (dim == 2) { done = ComputeAllElementTargets<2>(fes, ir, xe, Jtr); }
139 if (dim == 3) { done = ComputeAllElementTargets<3>(fes, ir, xe, Jtr); }
140
141 if (!done) { ComputeAllElementTargets_Fallback(fes, ir, xe, Jtr); }
142}
143
145 const IntegrationRule &ir,
146 const Vector &xe,
147 DenseTensor &Jtr) const
148{
149 ComputeAllElementTargets_Fallback(fes, ir, xe, Jtr);
150}
151
152
153// Code paths leading to ComputeElementTargets:
154// - GetElementEnergy(elfun) which is done through GetGridFunctionEnergyPA(x)
155// - AssembleElementVectorExact(elfun)
156// - AssembleElementGradExact(elfun)
157// - EnableNormalization(x) -> ComputeNormalizationEnergies(x)
158// - (AssembleElementVectorFD(elfun))
159// - (AssembleElementGradFD(elfun))
160// ============================================================================
161// - TargetConstructor():
162// - IDEAL_SHAPE_UNIT_SIZE: Wideal
163// - IDEAL_SHAPE_EQUAL_SIZE: α * Wideal
164// - IDEAL_SHAPE_GIVEN_SIZE: β * Wideal
165// - GIVEN_SHAPE_AND_SIZE: β * Wideal
166// - AnalyticAdaptTC(elfun):
167// - GIVEN_FULL: matrix_tspec->Eval(Jtr(elfun))
168// - DiscreteAdaptTC():
169// - IDEAL_SHAPE_GIVEN_SIZE: size^{1.0/dim} * Jtr(i) (size)
170// - GIVEN_SHAPE_AND_SIZE: Jtr(i) *= D_rho (ratio)
171// Jtr(i) *= Q_phi (skew)
172// Jtr(i) *= R_theta (orientation)
174{
175 PA.Jtr_needs_update = false;
176 PA.Jtr_debug_grad = false;
177 const FiniteElementSpace *fes = PA.fes;
178 if (PA.ne == 0) { return; } // Quick return for empty processors
179 const IntegrationRule &ir = *PA.ir;
180
181 // Compute PA.Jtr for all elements
183}
184
186{
187 Vector x_loc;
188 if (periodic)
189 {
190 GetPeriodicPositions(*x_0, d_loc, *x_0->FESpace(), *PA.fes, x_loc);
191 }
192 else
193 {
194 x_loc.SetSize(x_0->Size());
195 add(*x_0, d_loc, x_loc);
196 }
197
198 // Both are constant or not specified.
199 if (PA.MC.Size() == 1 && PA.C0.Size() == 1) { return; }
200
201 // Coefficients are always evaluated on the CPU for now.
202 PA.MC.HostWrite();
203 PA.C0.HostWrite();
204
205 const IntegrationRule &ir = *PA.ir;
206 auto T = new IsoparametricTransformation;
207 for (int e = 0; e < PA.ne; ++e)
208 {
209 // Uses the node positions in x_loc.
210 PA.fes->GetMesh()->GetElementTransformation(e, x_loc, T);
211
212 if (PA.MC.Size() > 1)
213 {
214 for (int q = 0; q < PA.nq; ++q)
215 {
216 PA.MC(q + e * PA.nq) = metric_coeff->Eval(*T, ir.IntPoint(q));
217 }
218 }
219
220 if (PA.C0.Size() > 1)
221 {
222 for (int q = 0; q < PA.nq; ++q)
223 {
224 PA.C0(q + e * PA.nq) = lim_coeff->Eval(*T, ir.IntPoint(q));
225 }
226 }
227 }
228
229 delete T;
230}
231
233{
234 const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
236 PA.enabled = true;
237 PA.fes = &fes;
238 Mesh *mesh = fes.GetMesh();
239 const int ne = PA.ne = mesh->GetNE();
240 const int dim = PA.dim = mesh->Dimension();
241
242 MFEM_VERIFY(PA.dim == 2 || PA.dim == 3, "Not yet implemented!");
243 MFEM_VERIFY(mesh->GetNumGeometries(dim) <= 1,
244 "TMOP+PA does not support mixed meshes.");
245 MFEM_VERIFY(mesh->HasGeometry(Geometry::SQUARE) ||
247 "TMOP+PA only supports squares and cubes.");
248 MFEM_VERIFY(!fes.IsVariableOrder(), "variable orders are not supported");
249 MFEM_VERIFY(fes.GetOrdering() == Ordering::byNODES,
250 "TMOP+PAP only supports Ordering::byNODES!");
251
252 const FiniteElement &fe = *fes.GetTypicalFE();
253 PA.ir = &EnergyIntegrationRule(fe);
254 const IntegrationRule &ir = *PA.ir;
255 const int nq = PA.nq = ir.GetNPoints();
257 PA.maps = &fe.GetDofToQuad(ir, mode);
258 // Note - initial mesh. TODO delete this?
260
261 // Initial node positions.
262 // PA.X0 is also updated in TMOP_Integrator::SetInitialMeshPos.
263 if (x_0 != nullptr)
264 {
266 auto n0_R = x_0->FESpace()->GetElementRestriction(ord);
267 PA.X0.UseDevice(true);
268 PA.X0.SetSize(n0_R->Height(), Device::GetMemoryType());
269 n0_R->Mult(*x_0, PA.X0);
270 }
271
272 // Energy vector, scalar Q-vector
273 PA.E.UseDevice(true);
274 PA.E.SetSize(ne*nq, Device::GetDeviceMemoryType());
275
276 // H for Grad, (dim x dim) Q-vector
277 PA.H.UseDevice(true);
278 PA.H.SetSize(dim*dim * dim*dim * nq*ne, mt);
279
280 // Scalar Q-vector of '1', used to compute sums via dot product
281 PA.O.SetSize(ne*nq, Device::GetDeviceMemoryType());
282 PA.O = 1.0;
283
284 if (metric_coeff)
285 {
286 if (auto cc = dynamic_cast<ConstantCoefficient *>(metric_coeff))
287 {
288 PA.MC.SetSize(1, Device::GetMemoryType());
289 PA.MC.HostWrite();
290 PA.MC(0) = cc->constant;
291 }
292 else
293 {
294 PA.MC.SetSize(PA.nq * PA.ne, Device::GetMemoryType());
295 auto M0 = Reshape(PA.MC.HostWrite(), PA.nq, PA.ne);
296 for (int e = 0; e < PA.ne; ++e)
297 {
298 ElementTransformation& T = *PA.fes->GetElementTransformation(e);
299 for (int q = 0; q < ir.GetNPoints(); ++q)
300 {
301 // Note that this is always on the initial mesh.
302 M0(q,e) = metric_coeff->Eval(T, ir.IntPoint(q));
303 }
304 }
305 }
306 }
307 else
308 {
309 PA.MC.SetSize(1, Device::GetMemoryType());
310 PA.MC.HostWrite();
311 PA.MC(0) = 1.0;
312 }
313
314 // Setup ref->target Jacobians, PA.Jtr, (dim x dim) Q-vector, DenseTensor
315 PA.Jtr.SetSize(dim, dim, PA.ne*PA.nq, mt);
316 PA.Jtr_needs_update = true;
317 PA.Jtr_debug_grad = false;
318
319 // Limiting: lim_coeff -> PA.C0, lim_nodes0 -> PA.XL, lim_dist -> PA.LD, PA.H0
321}
322
324{
325 // This method must be called after AssembleGradPA().
326
327 MFEM_VERIFY(PA.Jtr_needs_update == false, "");
328
330 {
331 MFEM_VERIFY(PA.Jtr_debug_grad == true, "AssembleGradPA() was not called"
332 " or Jtr was overwritten by another method!");
333 }
334
335 if (PA.dim == 2)
336 {
339 }
340
341 if (PA.dim == 3)
342 {
345 }
346}
347
348void TMOP_Integrator::AddMultPA(const Vector &de, Vector &ye) const
349{
350 // This method must be called after AssemblePA().
351
352 // Form the Vector of node positions, depending on what's the input.
353 Vector xe(de.Size());
354 if (x_0)
355 {
356 // The input is the displacement.
357 add(PA.X0, de, xe);
358 }
359 else { xe = de; }
360
361 if (PA.Jtr_needs_update || targetC->UsesPhysicalCoordinates())
362 {
364 }
365
366 if (PA.dim == 2)
367 {
368 AddMultPA_2D(xe,ye);
369 if (lim_coeff) { AddMultPA_C0_2D(xe,ye); }
370 }
371
372 if (PA.dim == 3)
373 {
374 AddMultPA_3D(xe,ye);
375 if (lim_coeff) { AddMultPA_C0_3D(xe,ye); }
376 }
377}
378
380{
381 // This method must be called after AssembleGradPA().
382
383 MFEM_VERIFY(PA.Jtr_needs_update == false, "");
384
386 {
387 MFEM_VERIFY(PA.Jtr_debug_grad == true, "AssembleGradPA() was not called or "
388 "Jtr was overwritten by another method!");
389 }
390
391 if (PA.dim == 2)
392 {
393 AddMultGradPA_2D(re,ce);
394 if (lim_coeff) { AddMultGradPA_C0_2D(re,ce); }
395 }
396
397 if (PA.dim == 3)
398 {
399 AddMultGradPA_3D(re,ce);
400 if (lim_coeff) { AddMultGradPA_C0_3D(re,ce); }
401 }
402}
403
405{
406 // This method must be called after AssemblePA().
407
408 // Form the Vector of node positions, depending on what's the input.
409 Vector xe(de.Size());
410 if (x_0)
411 {
412 // The input is the displacement.
413 add(PA.X0, de, xe);
414 }
415 else { xe = de; }
416
417 real_t energy = 0.0;
418
419 if (PA.Jtr_needs_update || targetC->UsesPhysicalCoordinates())
420 {
422 }
423
424 if (PA.dim == 2)
425 {
426 energy = GetLocalStateEnergyPA_2D(xe);
427 if (lim_coeff) { energy += GetLocalStateEnergyPA_C0_2D(xe); }
428 }
429
430 if (PA.dim == 3)
431 {
432 energy = GetLocalStateEnergyPA_3D(xe);
433 if (lim_coeff) { energy += GetLocalStateEnergyPA_C0_3D(xe); }
434 }
435
436 return energy;
437}
438
439} // namespace mfem
void ComputeAllElementTargets(const FiniteElementSpace &fes, const IntegrationRule &ir, const Vector &xe, DenseTensor &Jtr) const override
Computes reference-to-target transformation Jacobians for all quadrature points in all elements.
Definition tmop_pa.cpp:144
bool UseDevice() const
Return the device flag of the Memory object used by the Array.
Definition array.hpp:132
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
A coefficient that is constant across space and time.
Rank 3 tensor (array of matrices)
int SizeJ() const
int SizeI() const
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition device.hpp:278
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition device.hpp:274
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
Definition fe_base.hpp:154
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:165
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition fespace.hpp:709
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition fespace.hpp:924
Ordering::Type GetOrdering() const
Return the ordering method.
Definition fespace.hpp:876
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition fespace.cpp:1510
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
Definition fespace.cpp:3871
Abstract class for all finite elements.
Definition fe_base.hpp:244
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:365
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:334
FiniteElementSpace * FESpace()
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
A standard isoparametric element transformation.
Definition eltrans.hpp:629
Mesh data type.
Definition mesh.hpp:64
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
Definition mesh.cpp:880
bool HasGeometry(Geometry::Type geom) const
Return true iff the given geom is encountered in the mesh. Geometries of dimensions lower than Dimens...
Definition mesh.hpp:1250
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition mesh.cpp:7243
Abstract operator.
Definition operator.hpp:25
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Exponential limiter function in TMOP_Integrator.
Definition tmop.hpp:1349
void AddMultGradPA(const Vector &, Vector &) const override
Method for partially assembled gradient action.
Definition tmop_pa.cpp:379
real_t GetLocalStateEnergyPA_C0_2D(const Vector &) const
TMOP_LimiterFunction * lim_func
Definition tmop.hpp:1920
void ComputeAllElementTargets(const Vector &xe=Vector()) const
Definition tmop_pa.cpp:173
void AddMultPA_C0_3D(const Vector &, Vector &) const
real_t GetLocalStateEnergyPA_2D(const Vector &) const
void AssembleDiagonalPA_2D(Vector &) const
Coefficient * lim_coeff
Definition tmop.hpp:1916
const GridFunction * lim_dist
Definition tmop.hpp:1918
void AssembleDiagonalPA_C0_3D(Vector &) const
void AddMultGradPA_2D(const Vector &, Vector &) const
void AssembleGradPA_C0_2D(const Vector &) const
const TargetConstructor * targetC
Definition tmop.hpp:1900
const GridFunction * lim_nodes0
Definition tmop.hpp:1915
const IntegrationRule * ir
Definition tmop.hpp:2024
void AssemblePA(const FiniteElementSpace &) override
Method defining partial assembly.
Definition tmop_pa.cpp:232
void AssembleGradDiagonalPA(Vector &) const override
Method for computing the diagonal of the gradient with partial assembly.
Definition tmop_pa.cpp:323
void AssembleDiagonalPA_C0_2D(Vector &) const
real_t GetLocalStateEnergyPA_C0_3D(const Vector &) const
void AddMultGradPA_3D(const Vector &, Vector &) const
const FiniteElementSpace * fes
Definition tmop.hpp:2023
const IntegrationRule & EnergyIntegrationRule(const FiniteElement &el) const
Definition tmop.hpp:2085
void AddMultPA_2D(const Vector &, Vector &) const
void AddMultPA_3D(const Vector &, Vector &) const
void AssembleGradPA_C0_3D(const Vector &) const
void AssembleGradPA_2D(const Vector &) const
struct mfem::TMOP_Integrator::@26 PA
Coefficient * metric_coeff
Definition tmop.hpp:1908
const GridFunction * x_0
Definition tmop.hpp:1895
void AddMultPA(const Vector &, Vector &) const override
Method for partially assembled action.
Definition tmop_pa.cpp:348
void AddMultGradPA_C0_2D(const Vector &, Vector &) const
void AssembleDiagonalPA_3D(Vector &) const
void UpdateCoefficientsPA(const Vector &d_loc)
Definition tmop_pa.cpp:185
void AssembleGradPA_3D(const Vector &) const
void AssembleGradPA(const Vector &, const FiniteElementSpace &) override
Prepare the integrator for partial assembly (PA) gradient evaluations on the given FE space fes at th...
Definition tmop_pa.cpp:23
real_t GetLocalStateEnergyPA_3D(const Vector &) const
void AddMultPA_C0_2D(const Vector &, Vector &) const
real_t GetLocalStateEnergyPA(const Vector &) const override
Compute the local (to the MPI rank) energy with partial assembly.
Definition tmop_pa.cpp:404
void AddMultGradPA_C0_3D(const Vector &, Vector &) const
Default limiter function in TMOP_Integrator.
Definition tmop.hpp:1318
bool UsesPhysicalCoordinates() const
Return true if the methods ComputeElementTargets(), ComputeAllElementTargets(), and ComputeElementTar...
Definition tmop.hpp:1573
void ComputeAllElementTargets_Fallback(const FiniteElementSpace &fes, const IntegrationRule &ir, const Vector &xe, DenseTensor &Jtr) const
Definition tmop.cpp:2187
bool ComputeAllElementTargets(const FiniteElementSpace &fes, const IntegrationRule &ir, const Vector &xe, DenseTensor &Jtr) const
Vector data type.
Definition vector.hpp:82
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
int dim
Definition ex24.cpp:53
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:391
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition dtensor.hpp:131
void GetPeriodicPositions(const Vector &x_0, const Vector &dx, const FiniteElementSpace &fesL2, const FiniteElementSpace &fesH1, Vector &x)
float real_t
Definition config.hpp:43
MemoryType
Memory types supported by MFEM.
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition fespace.hpp:83