MFEM  v4.6.0
Finite element discretization library
nonlinearform.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 MFEM_NONLINEARFORM
13 #define MFEM_NONLINEARFORM
14 
15 #include "../config/config.hpp"
16 #include "nonlininteg.hpp"
17 #include "nonlinearform_ext.hpp"
18 #include "bilinearform.hpp"
19 #include "gridfunc.hpp"
20 
21 namespace mfem
22 {
23 
24 class NonlinearForm : public Operator
25 {
26 protected:
27  /// The assembly level.
29 
30  /// Extension for supporting different AssemblyLevel%s
31  /** For nonlinear operators, the "matrix" assembly levels usually do not make
32  sense, so only PARTIAL and NONE (matrix-free) are supported. */
34 
35  /// FE space on which the form lives.
36  FiniteElementSpace *fes; // not owned
37 
38  /// Set of Domain Integrators to be assembled (added).
40 
41  /// Set of interior face Integrators to be assembled (added).
43 
44  /// Set of boundary face Integrators to be assembled (added).
47 
48  mutable SparseMatrix *Grad, *cGrad; // owned
49  /// Gradient Operator when not assembled as a matrix.
50  mutable OperatorHandle hGrad; // has internal ownership flag
51 
52  /// A list of all essential true dofs
54 
55  /// Counter for updates propagated from the FiniteElementSpace.
56  long sequence;
57 
58  /// Auxiliary Vector%s
59  mutable Vector aux1, aux2;
60 
61  /// Pointer to the prolongation matrix of fes, may be NULL.
62  const Operator *P; // not owned
63  /// The result of dynamic-casting P to SparseMatrix pointer.
64  const SparseMatrix *cP; // not owned
65 
66  bool Serial() const { return (!P || cP); }
67  const Vector &Prolongate(const Vector &x) const;
68 
69 public:
70  /// Construct a NonlinearForm on the given FiniteElementSpace, @a f.
71  /** As an Operator, the NonlinearForm has input and output size equal to the
72  number of true degrees of freedom, i.e. f->GetTrueVSize(). */
74  : Operator(f->GetTrueVSize()), assembly(AssemblyLevel::LEGACY),
75  ext(NULL), fes(f), Grad(NULL), cGrad(NULL),
76  sequence(f->GetSequence()), P(f->GetProlongationMatrix()),
77  cP(dynamic_cast<const SparseMatrix*>(P))
78  { }
79 
80  /// Set the desired assembly level. The default is AssemblyLevel::LEGACY.
81  /** For nonlinear operators, the "matrix" assembly levels usually do not make
82  sense, so only LEGACY, NONE (matrix-free) and PARTIAL are supported.
83 
84  Currently, AssemblyLevel::LEGACY uses the standard nonlinear action
85  methods like AssembleElementVector of the NonlinearFormIntegrator class
86  which work only on CPU and do not utilize features such as fast
87  tensor-product basis evaluations. In this mode, the gradient operator is
88  constructed as a SparseMatrix (or, in parallel, format such as
89  HypreParMatrix).
90 
91  When using AssemblyLevel::PARTIAL, the action is performed using methods
92  like AddMultPA of the NonlinearFormIntegrator class which typically
93  support both CPU and GPU backends and utilize features such as fast
94  tensor-product basis evaluations. In this mode, the gradient operator
95  also uses partial assembly with support for CPU and GPU backends.
96 
97  When using AssemblyLevel::NONE, the action is performed using methods
98  like AddMultMF of the NonlinearFormIntegrator class which typically
99  support both CPU and GPU backends and utilize features such as fast
100  tensor-product basis evaluations. In this mode, the gradient operator
101  is currently not supported.
102 
103  This method must be called before "assembly" with Setup(). */
104  void SetAssemblyLevel(AssemblyLevel assembly_level);
105 
107  const FiniteElementSpace *FESpace() const { return fes; }
108 
109  /// Adds new Domain Integrator.
111  { dnfi.Append(nlfi); }
112 
113  /// Access all integrators added with AddDomainIntegrator().
115  const Array<NonlinearFormIntegrator*> *GetDNFI() const { return &dnfi; }
116 
117  /// Adds new Interior Face Integrator.
119  { fnfi.Append(nlfi); }
120 
121  /** @brief Access all interior face integrators added with
122  AddInteriorFaceIntegrator(). */
124  { return fnfi; }
125 
126  /// Adds new Boundary Face Integrator.
128  { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
129 
130  /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
131  attributes. */
133  Array<int> &bdr_marker)
134  { bfnfi.Append(nfi); bfnfi_marker.Append(&bdr_marker); }
135 
136  /** @brief Access all boundary face integrators added with
137  AddBdrFaceIntegrator(). */
139  { return bfnfi; }
140 
141  /// Specify essential boundary conditions.
142  /** This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
143  the result internally for use by other methods. If the @a rhs pointer is
144  not NULL, its essential true dofs will be set to zero. This makes it
145  "compatible" with the output vectors from the Mult() method which also
146  have zero entries at the essential true dofs. */
147  void SetEssentialBC(const Array<int> &bdr_attr_is_ess, Vector *rhs = NULL);
148 
149  /// Specify essential boundary conditions.
150  /** Use either SetEssentialBC() or SetEssentialTrueDofs() if possible. */
151  void SetEssentialVDofs(const Array<int> &ess_vdofs_list);
152 
153  /// Specify essential boundary conditions.
154  void SetEssentialTrueDofs(const Array<int> &ess_tdof_list_)
155  { ess_tdof_list_.Copy(this->ess_tdof_list); }
156 
157  /// Return a (read-only) list of all essential true dofs.
158  const Array<int> &GetEssentialTrueDofs() const { return ess_tdof_list; }
159 
160  /// Compute the energy corresponding to the state @a x.
161  /** In general, @a x may have non-homogeneous essential boundary values.
162 
163  The state @a x must be a "GridFunction size" vector, i.e. its size must
164  be fes->GetVSize(). */
165  double GetGridFunctionEnergy(const Vector &x) const;
166 
167  /// Compute the energy corresponding to the state @a x.
168  /** In general, @a x may have non-homogeneous essential boundary values.
169 
170  The state @a x must be a true-dof vector. */
171  virtual double GetEnergy(const Vector &x) const
172  { return GetGridFunctionEnergy(Prolongate(x)); }
173 
174  /// Evaluate the action of the NonlinearForm.
175  /** The input essential dofs in @a x will, generally, be non-zero. However,
176  the output essential dofs in @a y will always be set to zero.
177 
178  Both the input and the output vectors, @a x and @a y, must be true-dof
179  vectors, i.e. their size must be fes->GetTrueVSize(). */
180  virtual void Mult(const Vector &x, Vector &y) const;
181 
182  /** @brief Compute the gradient Operator of the NonlinearForm corresponding
183  to the state @a x. */
184  /** Any previously specified essential boundary conditions will be
185  automatically imposed on the gradient operator.
186 
187  The returned object is valid until the next call to this method or the
188  destruction of this object.
189 
190  In general, @a x may have non-homogeneous essential boundary values.
191 
192  The state @a x must be a true-dof vector. */
193  virtual Operator &GetGradient(const Vector &x) const;
194 
195  /// Update the NonlinearForm to propagate updates of the associated FE space.
196  /** After calling this method, the essential boundary conditions need to be
197  set again. */
198  virtual void Update();
199 
200  /** @brief Setup the NonlinearForm: based on the current AssemblyLevel and
201  the current mesh, optionally, precompute and store data that will be
202  reused in subsequent call to Mult(). */
203  /** Typically, this method has to be called before Mult() when using
204  AssemblyLevel::PARTIAL, after calling Update(), or after modifying the
205  mesh coordinates. */
206  virtual void Setup();
207 
208  /// Get the finite element space prolongation matrix
209  virtual const Operator *GetProlongation() const { return P; }
210  /// Get the finite element space restriction matrix
211  virtual const Operator *GetRestriction() const
212  { return fes->GetRestrictionMatrix(); }
213 
214  /** @brief Destroy the NonlinearForm including the owned
215  NonlinearFormIntegrator%s and gradient Operator. */
216  virtual ~NonlinearForm();
217 };
218 
219 
220 /** @brief A class representing a general block nonlinear operator defined on
221  the Cartesian product of multiple FiniteElementSpace%s. */
223 {
224 protected:
225  /// FE spaces on which the form lives.
227 
228  /// Set of Domain Integrators to be assembled (added).
230 
231  /// Set of interior face Integrators to be assembled (added).
233 
234  /// Set of Boundary Face Integrators to be assembled (added).
237 
238  /** Auxiliary block-vectors for wrapping input and output vectors or holding
239  GridFunction-like block-vector data (e.g. in parallel). */
240  mutable BlockVector xs, ys;
241 
244 
245  // A list of the offsets
248 
249  // Array of Arrays of tdofs for each space in 'fes'
251 
252  /// Array of pointers to the prolongation matrix of fes, may be NULL
254 
255  /// Array of results of dynamic-casting P to SparseMatrix pointer
257 
258  /// Indicator if the Operator is part of a parallel run
259  bool is_serial = true;
260 
261  /// Indicator if the Operator needs prolongation on assembly
262  bool needs_prolongation = false;
263 
265 
266  const BlockVector &Prolongate(const BlockVector &bx) const;
267 
268  /// Specialized version of GetEnergy() for BlockVectors
269  double GetEnergyBlocked(const BlockVector &bx) const;
270 
271  /// Specialized version of Mult() for BlockVector%s
272  /// Block L-Vector to Block L-Vector
273  void MultBlocked(const BlockVector &bx, BlockVector &by) const;
274 
275  /// Specialized version of GetGradient() for BlockVector
276  void ComputeGradientBlocked(const BlockVector &bx) const;
277 
278 public:
279  /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
281 
282  /// Construct a BlockNonlinearForm on the given set of FiniteElementSpace%s.
284 
285  /// Return the @a k-th FE space of the BlockNonlinearForm.
286  FiniteElementSpace *FESpace(int k) { return fes[k]; }
287  /// Return the @a k-th FE space of the BlockNonlinearForm (const version).
288  const FiniteElementSpace *FESpace(int k) const { return fes[k]; }
289 
290  /// (Re)initialize the BlockNonlinearForm.
291  /** After a call to SetSpaces(), the essential b.c. must be set again. */
293 
294  /// Return the regular dof offsets.
295  const Array<int> &GetBlockOffsets() const { return block_offsets; }
296  /// Return the true-dof offsets.
298 
299  /// Adds new Domain Integrator.
301  { dnfi.Append(nlfi); }
302 
303  /// Adds new Interior Face Integrator.
305  { fnfi.Append(nlfi); }
306 
307  /// Adds new Boundary Face Integrator.
309  { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
310 
311  /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
312  attributes. */
314  Array<int> &bdr_marker);
315 
316  virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
317  Array<Vector *> &rhs);
318 
319  virtual double GetEnergy(const Vector &x) const;
320 
321  /// Method is only called in serial, the parallel version calls MultBlocked
322  /// directly.
323  virtual void Mult(const Vector &x, Vector &y) const;
324 
325  /// Method is only called in serial, the parallel version calls
326  /// GetGradientBlocked directly.
327  virtual Operator &GetGradient(const Vector &x) const;
328 
329  /// Destructor.
330  virtual ~BlockNonlinearForm();
331 };
332 
333 
334 }
335 
336 #endif
const FiniteElementSpace * FESpace() const
Array< const Operator * > P
Array of pointers to the prolongation matrix of fes, may be NULL.
BlockNonlinearForm()
Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
bool is_serial
Indicator if the Operator is part of a parallel run.
const Vector & Prolongate(const Vector &x) const
Array< NonlinearFormIntegrator * > dnfi
Set of Domain Integrators to be assembled (added).
const Array< NonlinearFormIntegrator * > & GetBdrFaceIntegrators() const
Access all boundary face integrators added with AddBdrFaceIntegrator().
const Array< int > & GetBlockOffsets() const
Return the regular dof offsets.
Array< BlockNonlinearFormIntegrator * > dnfi
Set of Domain Integrators to be assembled (added).
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
void AddDomainIntegrator(BlockNonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
SparseMatrix * cGrad
Array2D< SparseMatrix * > cGrads
virtual void Mult(const Vector &x, Vector &y) const
virtual void Setup()
Setup the NonlinearForm: based on the current AssemblyLevel and the current mesh, optionally...
void ComputeGradientBlocked(const BlockVector &bx) const
Specialized version of GetGradient() for BlockVector.
void SetSpaces(Array< FiniteElementSpace *> &f)
(Re)initialize the BlockNonlinearForm.
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void AddBdrFaceIntegrator(NonlinearFormIntegrator *nfi, Array< int > &bdr_marker)
Adds new Boundary Face Integrator, restricted to specific boundary attributes.
A class representing a general block nonlinear operator defined on the Cartesian product of multiple ...
Array< NonlinearFormIntegrator * > bfnfi
Set of boundary face Integrators to be assembled (added).
const Array< int > & GetBlockTrueOffsets() const
Return the true-dof offsets.
Array< BlockNonlinearFormIntegrator * > bfnfi
Set of Boundary Face Integrators to be assembled (added).
Array< int > ess_tdof_list
A list of all essential true dofs.
FiniteElementSpace * fes
FE space on which the form lives.
void AddInteriorFaceIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Interior Face Integrator.
void AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Boundary Face Integrator.
Array< const SparseMatrix * > cP
Array of results of dynamic-casting P to SparseMatrix pointer.
Array2D< SparseMatrix * > Grads
virtual void Mult(const Vector &x, Vector &y) const
Evaluate the action of the NonlinearForm.
void MultBlocked(const BlockVector &bx, BlockVector &by) const
const Array< NonlinearFormIntegrator * > * GetDNFI() const
std::function< double(const Vector &)> f(double mass_coeff)
Definition: lor_mms.hpp:30
void AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
Adds new Boundary Face Integrator.
Array< FiniteElementSpace * > fes
FE spaces on which the form lives.
Data type sparse matrix.
Definition: sparsemat.hpp:50
void AddInteriorFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
Adds new Interior Face Integrator.
const Array< int > & GetEssentialTrueDofs() const
Return a (read-only) list of all essential true dofs.
Array< NonlinearFormIntegrator * > * GetDNFI()
Access all integrators added with AddDomainIntegrator().
Array< NonlinearFormIntegrator * > fnfi
Set of interior face Integrators to be assembled (added).
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:616
NonlinearForm(FiniteElementSpace *f)
Construct a NonlinearForm on the given FiniteElementSpace, f.
virtual Operator & GetGradient(const Vector &x) const
const SparseMatrix * cP
The result of dynamic-casting P to SparseMatrix pointer.
Vector aux1
Auxiliary Vectors.
virtual double GetEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
Array< BlockNonlinearFormIntegrator * > fnfi
Set of interior face Integrators to be assembled (added).
AssemblyLevel assembly
The assembly level.
Class extending the NonlinearForm class to support the different AssemblyLevels.
virtual void SetEssentialBC(const Array< Array< int > *> &bdr_attr_is_ess, Array< Vector *> &rhs)
const Operator * P
Pointer to the prolongation matrix of fes, may be NULL.
virtual double GetEnergy(const Vector &x) const
Array< Array< int > * > bfnfi_marker
double GetGridFunctionEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
Array< Array< int > * > bfnfi_marker
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACY.
Dynamic 2D array using row-major layout.
Definition: array.hpp:354
FiniteElementSpace * FESpace()
virtual void Update()
Update the NonlinearForm to propagate updates of the associated FE space.
virtual Operator & GetGradient(const Vector &x) const
Compute the gradient Operator of the NonlinearForm corresponding to the state x.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
NonlinearFormExtension * ext
Extension for supporting different AssemblyLevels.
const BlockVector & Prolongate(const BlockVector &bx) const
virtual ~BlockNonlinearForm()
Destructor.
const FiniteElementSpace * FESpace(int k) const
Return the k-th FE space of the BlockNonlinearForm (const version).
virtual ~NonlinearForm()
Destroy the NonlinearForm including the owned NonlinearFormIntegrators and gradient Operator...
double GetEnergyBlocked(const BlockVector &bx) const
Specialized version of GetEnergy() for BlockVectors.
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list_)
Specify essential boundary conditions.
virtual const Operator * GetRestriction() const
Get the finite element space restriction matrix.
This class is used to express the local action of a general nonlinear finite element operator...
Definition: nonlininteg.hpp:27
const Array< NonlinearFormIntegrator * > & GetInteriorFaceIntegrators() const
Access all interior face integrators added with AddInteriorFaceIntegrator().
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:864
long sequence
Counter for updates propagated from the FiniteElementSpace.
OperatorHandle hGrad
Gradient Operator when not assembled as a matrix.
bool needs_prolongation
Indicator if the Operator needs prolongation on assembly.
void SetEssentialBC(const Array< int > &bdr_attr_is_ess, Vector *rhs=NULL)
Specify essential boundary conditions.
Vector data type.
Definition: vector.hpp:58
SparseMatrix * Grad
virtual const Operator * GetProlongation() const
Get the finite element space prolongation matrix.
void SetEssentialVDofs(const Array< int > &ess_vdofs_list)
Specify essential boundary conditions.
Array< Array< int > * > ess_tdofs
Abstract operator.
Definition: operator.hpp:24
A class to handle Block systems in a matrix-free implementation.
FiniteElementSpace * FESpace(int k)
Return the k-th FE space of the BlockNonlinearForm.