1 // Copyright (c) 2010-2021, 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_PRMNONLINEARFORM 13 #define MFEM_PRMNONLINEARFORM 14 15 #include "mfem.hpp" 16 17 namespace mfem 18 { 19 20 /** The abstract base class ParametricBNLFormIntegrator is a generalization of 21 the BlockNonlinearFormIntegrator class suitable for block state and 22 parameter vectors. */ 23 class ParametricBNLFormIntegrator 24 { 25 public: 26 /// Compute the local energy 27 virtual double GetElementEnergy(const Array<const FiniteElement *>&el, 28 const Array<const FiniteElement *>&pel, 29 ElementTransformation &Tr, 30 const Array<const Vector *>&elfun, 31 const Array<const Vector *>&pelfun); 32 33 /// Perform the local action of the BlockNonlinearFormIntegrator 34 virtual void AssembleElementVector(const Array<const FiniteElement *> &el, 35 const Array<const FiniteElement *>&pel, 36 ElementTransformation &Tr, 37 const Array<const Vector *> &elfun, 38 const Array<const Vector *>&pelfun, 39 const Array<Vector *> &elvec); 40 41 /// Perform the local action of the BlockNonlinearFormIntegrator on element 42 /// faces 43 virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1, 44 const Array<const FiniteElement *> &el2, 45 const Array<const FiniteElement *> &pel1, 46 const Array<const FiniteElement *> &pel2, 47 FaceElementTransformations &Tr, 48 const Array<const Vector *> &elfun, 49 const Array<const Vector *>&pelfun, 50 const Array<Vector *> &elvect); 51 52 /// Perform the local action on the parameters of the BNLFormIntegrator 53 virtual void AssemblePrmElementVector(const Array<const FiniteElement *> &el, 54 const Array<const FiniteElement *>&pel, 55 ElementTransformation &Tr, 56 const Array<const Vector *> &elfun, 57 const Array<const Vector *> &alfun, 58 const Array<const Vector *>&pelfun, 59 const Array<Vector *> &pelvec); 60 61 /// Perform the local action on the parameters of the BNLFormIntegrator on 62 /// faces 63 virtual void AssemblePrmFaceVector(const Array<const FiniteElement *> &el1, 64 const Array<const FiniteElement *> &el2, 65 const Array<const FiniteElement *> &pel1, 66 const Array<const FiniteElement *> &pel2, 67 FaceElementTransformations &Tr, 68 const Array<const Vector *> &elfun, 69 const Array<const Vector *> &alfun, 70 const Array<const Vector *>&pelfun, 71 const Array<Vector *> &pelvect); 72 73 /// Assemble the local gradient matrix 74 virtual void AssembleElementGrad(const Array<const FiniteElement*> &el, 75 const Array<const FiniteElement *>&pel, 76 ElementTransformation &Tr, 77 const Array<const Vector *> &elfun, 78 const Array<const Vector *>&pelfun, 79 const Array2D<DenseMatrix *> &elmats); 80 81 /// Assemble the local gradient matrix on faces of the elements 82 virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1, 83 const Array<const FiniteElement *>&el2, 84 const Array<const FiniteElement *> &pel1, 85 const Array<const FiniteElement *> &pel2, 86 FaceElementTransformations &Tr, 87 const Array<const Vector *> &elfun, 88 const Array<const Vector *>&pelfun, 89 const Array2D<DenseMatrix *> &elmats); 90 91 ~ParametricBNLFormIntegrator()92 virtual ~ParametricBNLFormIntegrator() { } 93 }; 94 95 96 /** @brief A class representing a general parametric block nonlinear operator 97 defined on the Cartesian product of multiple FiniteElementSpace%s. */ 98 class ParametricBNLForm : public Operator 99 { 100 protected: 101 /// FE spaces on which the form lives. 102 Array<FiniteElementSpace*> fes; 103 104 /// FE spaces for the parametric fields 105 Array<FiniteElementSpace*> paramfes; 106 107 int paramheight; 108 int paramwidth; 109 110 /// Set of Domain Integrators to be assembled (added). 111 Array<ParametricBNLFormIntegrator*> dnfi; 112 113 /// Set of interior face Integrators to be assembled (added). 114 Array<ParametricBNLFormIntegrator*> fnfi; 115 116 /// Set of Boundary Face Integrators to be assembled (added). 117 Array<ParametricBNLFormIntegrator*> bfnfi; 118 Array<Array<int>*> bfnfi_marker; 119 120 /** Auxiliary block-vectors for wrapping input and output vectors or holding 121 GridFunction-like block-vector data (e.g. in parallel). */ 122 mutable BlockVector xs, ys; 123 mutable BlockVector prmxs, prmys; 124 125 /** Auxiliary block-vectors for holding GridFunction-like block-vector data 126 (e.g. in parallel). */ 127 mutable BlockVector xsv; 128 129 /** Auxiliary block-vectors for holding GridFunction-like block-vector data 130 for the parameter fields (e.g. in parallel). */ 131 mutable BlockVector xdv; 132 /** Auxiliary block-vectors for holding GridFunction-like block-vector data 133 for the adjoint fields (e.g. in parallel). */ 134 mutable BlockVector adv; 135 136 mutable Array2D<SparseMatrix*> Grads, cGrads; 137 mutable BlockOperator *BlockGrad; 138 139 // A list of the offsets 140 Array<int> block_offsets; 141 Array<int> block_trueOffsets; 142 // A list with the offsets for the parametric fields 143 Array<int> paramblock_offsets; 144 Array<int> paramblock_trueOffsets; 145 146 // Array of Arrays of tdofs for each space in 'fes' 147 Array<Array<int> *> ess_tdofs; 148 149 // Array of Arrays of tdofs for each space in 'paramfes' 150 Array<Array<int> *> paramess_tdofs; 151 152 /// Array of pointers to the prolongation matrix of fes, may be NULL 153 Array<const Operator *> P; 154 155 /// Array of pointers to the prolongation matrix of paramfes, may be NULL 156 Array<const Operator *> Pparam; 157 158 /// Array of results of dynamic-casting P to SparseMatrix pointer 159 Array<const SparseMatrix *> cP; 160 161 /// Array of results of dynamic-casting Pparam to SparseMatrix pointer 162 Array<const SparseMatrix *> cPparam; 163 164 /// Indicator if the Operator is part of a parallel run 165 bool is_serial = true; 166 167 /// Indicator if the Operator needs prolongation on assembly 168 bool needs_prolongation = false; 169 170 /// Indicator if the Operator needs prolongation on assembly 171 bool prmneeds_prolongation = false; 172 173 mutable BlockVector aux1, aux2; 174 175 mutable BlockVector prmaux1, prmaux2; 176 177 const BlockVector &Prolongate(const BlockVector &bx) const; 178 179 const BlockVector &ParamProlongate(const BlockVector &bx) const; 180 181 double GetEnergyBlocked(const BlockVector &bx, const BlockVector &dx) const; 182 183 184 /// Specialized version of Mult() for BlockVector%s 185 /// Block L-Vector to Block L-Vector 186 void MultBlocked(const BlockVector &bx, const BlockVector &dx, 187 BlockVector &by) const; 188 189 /// Specialized version of Mult() for BlockVector%s 190 /// Block L-Vector to Block L-Vector 191 /// bx - state vector, ax - adjoint vector, dx - parametric fields 192 /// dy = ax' d(residual(bx))/d(dx) 193 void MultParamBlocked(const BlockVector &bx, const BlockVector & ax, 194 const BlockVector &dx, BlockVector &dy) const; 195 196 197 /// Specialized version of GetGradient() for BlockVector 198 void ComputeGradientBlocked(const BlockVector &bx, const BlockVector &dx) const; 199 200 public: 201 /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces(). 202 ParametricBNLForm(); 203 204 /// Construct a BlockNonlinearForm on the given set of FiniteElementSpace%s. 205 ParametricBNLForm(Array<FiniteElementSpace *> &statef, 206 Array<FiniteElementSpace *> ¶mf); 207 208 /// Return the @a k-th FE space of the ParametricBNLForm. FESpace(int k)209 FiniteElementSpace *FESpace(int k) { return fes[k]; } 210 211 /// Return the @a k-th parametric FE space of the ParametricBNLForm. ParamFESpace(int k)212 FiniteElementSpace *ParamFESpace(int k) { return paramfes[k]; } 213 214 215 /// Return the @a k-th FE space of the BlockNonlinearForm (const version). FESpace(int k) const216 const FiniteElementSpace *FESpace(int k) const { return fes[k]; } 217 218 /// Return the @a k-th parametric FE space of the BlockNonlinearForm (const 219 /// version). ParamFESpace(int k) const220 const FiniteElementSpace *ParamFESpace(int k) const { return paramfes[k]; } 221 222 /// Return the integrators GetDNFI()223 Array<ParametricBNLFormIntegrator*>& GetDNFI() { return dnfi;} 224 225 226 /// (Re)initialize the ParametricBNLForm. 227 /** After a call to SetSpaces(), the essential b.c. must be set again. */ 228 void SetSpaces(Array<FiniteElementSpace *> &statef, 229 Array<FiniteElementSpace *> ¶mf); 230 231 /// Return the regular dof offsets. GetBlockOffsets() const232 const Array<int> &GetBlockOffsets() const { return block_offsets; } 233 234 /// Return the true-dof offsets. GetBlockTrueOffsets() const235 const Array<int> &GetBlockTrueOffsets() const { return block_trueOffsets; } 236 237 /// Return the regular dof offsets for the parameters. ParamGetBlockOffsets() const238 const Array<int> &ParamGetBlockOffsets() const { return paramblock_offsets; } 239 240 /// Return the true-dof offsets for the parameters. ParamGetBlockTrueOffsets() const241 const Array<int> &ParamGetBlockTrueOffsets() const { return paramblock_trueOffsets; } 242 243 /// Adds new Domain Integrator. AddDomainIntegrator(ParametricBNLFormIntegrator * nlfi)244 void AddDomainIntegrator(ParametricBNLFormIntegrator *nlfi) 245 { dnfi.Append(nlfi); } 246 247 /// Adds new Interior Face Integrator. AddInteriorFaceIntegrator(ParametricBNLFormIntegrator * nlfi)248 void AddInteriorFaceIntegrator(ParametricBNLFormIntegrator *nlfi) 249 { fnfi.Append(nlfi); } 250 251 /// Adds new Boundary Face Integrator. AddBdrFaceIntegrator(ParametricBNLFormIntegrator * nlfi)252 void AddBdrFaceIntegrator(ParametricBNLFormIntegrator *nlfi) 253 { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); } 254 255 /** @brief Adds new Boundary Face Integrator, restricted to specific boundary 256 attributes. */ 257 void AddBdrFaceIntegrator(ParametricBNLFormIntegrator *nlfi, 258 Array<int> &bdr_marker); 259 260 /// Set the essential boundary conditions. 261 virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess, 262 Array<Vector *> &rhs); 263 264 /// Set the essential boundary conditions on the parametric fields. 265 virtual void SetParamEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess, 266 Array<Vector *> &rhs); 267 268 269 /// Computes the energy for a state vector x. 270 virtual double GetEnergy(const Vector &x) const; 271 272 /// Method is only called in serial, the parallel version calls MultBlocked 273 /// directly. 274 virtual void Mult(const Vector &x, Vector &y) const; 275 276 /// Method is only called in serial, the parallel version calls MultBlocked 277 /// directly. 278 virtual void ParamMult(const Vector &x, Vector &y) const; 279 280 /// Method is only called in serial, the parallel version calls 281 /// GetGradientBlocked directly. 282 virtual BlockOperator &GetGradient(const Vector &x) const; 283 284 /// Set the state fields 285 virtual void SetStateFields(const Vector &xv) const; 286 287 /// Set the adjoint fields 288 virtual void SetAdjointFields(const Vector &av) const; 289 290 /// Set the parameters/design fields 291 virtual void SetParamFields(const Vector &dv) const; 292 293 /// Destructor. 294 virtual ~ParametricBNLForm(); 295 296 }; 297 298 } 299 300 #endif 301