1 // @HEADER 2 // ************************************************************************ 3 // 4 // Rapid Optimization Library (ROL) Package 5 // Copyright (2014) Sandia Corporation 6 // 7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 8 // license for use of this work by or on behalf of the U.S. Government. 9 // 10 // Redistribution and use in source and binary forms, with or without 11 // modification, are permitted provided that the following conditions are 12 // met: 13 // 14 // 1. Redistributions of source code must retain the above copyright 15 // notice, this list of conditions and the following disclaimer. 16 // 17 // 2. Redistributions in binary form must reproduce the above copyright 18 // notice, this list of conditions and the following disclaimer in the 19 // documentation and/or other materials provided with the distribution. 20 // 21 // 3. Neither the name of the Corporation nor the names of the 22 // contributors may be used to endorse or promote products derived from 23 // this software without specific prior written permission. 24 // 25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 36 // 37 // Questions? Contact lead developers: 38 // Drew Kouri (dpkouri@sandia.gov) and 39 // Denis Ridzal (dridzal@sandia.gov) 40 // 41 // ************************************************************************ 42 // @HEADER 43 44 /*! \file pde_poisson_topOpt.hpp 45 \brief Implements the local PDE interface for the Poisson 46 topology optimization problem. 47 */ 48 49 #ifndef PDE_POISSON_TOPOPT_HPP 50 #define PDE_POISSON_TOPOPT_HPP 51 52 #include "../../TOOLS/pde.hpp" 53 #include "../../TOOLS/fe.hpp" 54 55 #include "Intrepid_HGRAD_QUAD_C1_FEM.hpp" 56 #include "Intrepid_HGRAD_QUAD_C2_FEM.hpp" 57 #include "Intrepid_DefaultCubatureFactory.hpp" 58 #include "Intrepid_FunctionSpaceTools.hpp" 59 #include "Intrepid_RealSpaceTools.hpp" 60 #include "Intrepid_CellTools.hpp" 61 62 #include "ROL_Ptr.hpp" 63 64 template <class Real> 65 class PDE_Poisson_TopOpt : public PDE<Real> { 66 private: 67 // Finite element basis information 68 ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > basisPtr_; 69 std::vector<ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > > basisPtrs_; 70 // Cell cubature information 71 ROL::Ptr<Intrepid::Cubature<Real> > cellCub_; 72 // Cell node information 73 ROL::Ptr<Intrepid::FieldContainer<Real> > volCellNodes_; 74 std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real> > > > bdryCellNodes_; 75 std::vector<std::vector<std::vector<int> > > bdryCellLocIds_; 76 // Finite element definition 77 ROL::Ptr<FE<Real> > fe_vol_; 78 // Local degrees of freedom on boundary, for each side of the reference cell (first index). 79 std::vector<std::vector<int> > fidx_; 80 // Coordinates of degrees freedom on boundary cells. 81 // Indexing: [sideset number][local side id](cell number, value at dof) 82 std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real> > > > bdryCellDofValues_; 83 // Force function evaluated at cubature points 84 ROL::Ptr<Intrepid::FieldContainer<Real> > force_eval_; 85 ROL::Ptr<Intrepid::FieldContainer<Real> > nforce_eval_; 86 // Inputs 87 Real minConductivity_; 88 Real SIMPpower_; 89 ForceFunc(const std::vector<Real> & x) const90 Real ForceFunc(const std::vector<Real> &x) const { 91 return static_cast<Real>(0.01); 92 } 93 dirichletFunc(const std::vector<Real> & coords,const int sideset,const int locSideId) const94 Real dirichletFunc(const std::vector<Real> & coords, const int sideset, const int locSideId) const { 95 return static_cast<Real>(0); 96 } 97 98 public: PDE_Poisson_TopOpt(Teuchos::ParameterList & parlist)99 PDE_Poisson_TopOpt(Teuchos::ParameterList &parlist) { 100 // Finite element fields. 101 int basisOrder = parlist.sublist("Problem").get("Order of FE discretization",1); 102 if (basisOrder == 1) { 103 basisPtr_ = ROL::makePtr<Intrepid::Basis_HGRAD_QUAD_C1_FEM<Real, Intrepid::FieldContainer<Real> >>(); 104 } 105 else if (basisOrder == 2) { 106 basisPtr_ = ROL::makePtr<Intrepid::Basis_HGRAD_QUAD_C2_FEM<Real, Intrepid::FieldContainer<Real> >>(); 107 } 108 basisPtrs_.clear(); basisPtrs_.push_back(basisPtr_); 109 // Quadrature rules. 110 shards::CellTopology cellType = basisPtr_->getBaseCellTopology(); // get the cell type from any basis 111 Intrepid::DefaultCubatureFactory<Real> cubFactory; // create cubature factory 112 int cubDegree = parlist.sublist("Problem").get("Cubature Degree",2); // set cubature degree, e.g., 2 113 cellCub_ = cubFactory.create(cellType, cubDegree); // create default cubature 114 115 minConductivity_ = parlist.sublist("Problem").get("Minimum Conductivity",1.e-3); 116 SIMPpower_ = parlist.sublist("Problem").get("SIMP Power",3.0); 117 } 118 residual(ROL::Ptr<Intrepid::FieldContainer<Real>> & res,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)119 void residual(ROL::Ptr<Intrepid::FieldContainer<Real> > & res, 120 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 121 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 122 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 123 const Real one(1); 124 // Get dimensions 125 int c = fe_vol_->gradN()->dimension(0); 126 int f = fe_vol_->gradN()->dimension(1); 127 int p = fe_vol_->gradN()->dimension(2); 128 int d = fe_vol_->gradN()->dimension(3); 129 // Initialize residual storage 130 res = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f); 131 // Evaluate density at cubature points 132 ROL::Ptr<Intrepid::FieldContainer<Real> > valZ_eval = 133 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p); 134 fe_vol_->evaluateValue(valZ_eval, z_coeff); 135 // Build SIMP density at cubature points 136 for (int i = 0; i < c; ++i) { 137 for (int j = 0; j < p; ++j) { 138 (*valZ_eval)(i,j) = minConductivity_ 139 + (one - minConductivity_) * std::pow((*valZ_eval)(i,j),SIMPpower_); 140 } 141 } 142 // Build flux function 143 ROL::Ptr<Intrepid::FieldContainer<Real> > gradU_eval = 144 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p, d); 145 fe_vol_->evaluateGradient(gradU_eval, u_coeff); 146 Intrepid::FieldContainer<Real> KgradU(c,p,d); 147 Intrepid::FunctionSpaceTools::scalarMultiplyDataData<Real>(KgradU, 148 *valZ_eval, 149 *gradU_eval); 150 // Integrate stiffness term 151 Intrepid::FunctionSpaceTools::integrate<Real>(*res, 152 KgradU, 153 *(fe_vol_->gradNdetJ()), 154 Intrepid::COMP_CPP, false); 155 // Add force term 156 Intrepid::FunctionSpaceTools::integrate<Real>(*res, 157 *nforce_eval_, 158 *(fe_vol_->NdetJ()), 159 Intrepid::COMP_CPP, true); 160 // Apply Dirichlet boundary conditions 161 int numSideSets = bdryCellLocIds_.size(); 162 if (numSideSets > 0) { 163 for (int i = 0; i < numSideSets; ++i) { 164 if ((i==2)) { 165 int numLocalSideIds = bdryCellLocIds_[i].size(); 166 for (int j = 0; j < numLocalSideIds; ++j) { 167 int numCellsSide = bdryCellLocIds_[i][j].size(); 168 int numBdryDofs = fidx_[j].size(); 169 for (int k = 0; k < numCellsSide; ++k) { 170 int cidx = bdryCellLocIds_[i][j][k]; 171 for (int l = 0; l < numBdryDofs; ++l) { 172 (*res)(cidx,fidx_[j][l]) = (*u_coeff)(cidx,fidx_[j][l]) - (*bdryCellDofValues_[i][j])(k,fidx_[j][l]); 173 } 174 } 175 } 176 } 177 } 178 } 179 } 180 Jacobian_1(ROL::Ptr<Intrepid::FieldContainer<Real>> & jac,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)181 void Jacobian_1(ROL::Ptr<Intrepid::FieldContainer<Real> > & jac, 182 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 183 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 184 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 185 const Real one(1); 186 // Get dimensions 187 int c = fe_vol_->gradN()->dimension(0); 188 int f = fe_vol_->gradN()->dimension(1); 189 int p = fe_vol_->gradN()->dimension(2); 190 int d = fe_vol_->gradN()->dimension(3); 191 // Initialize Jacobian storage 192 jac = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, f); 193 // Build density-dependent conductivity function 194 ROL::Ptr<Intrepid::FieldContainer<Real> > valZ_eval = 195 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p); 196 fe_vol_->evaluateValue(valZ_eval, z_coeff); 197 for (int i = 0; i < c; ++i) { 198 for (int j = 0; j < p; ++j) { 199 (*valZ_eval)(i,j) = minConductivity_ 200 + (one - minConductivity_) * std::pow((*valZ_eval)(i,j),SIMPpower_); 201 } 202 } 203 // Build flux function 204 Intrepid::FieldContainer<Real> KgradN(c,f,p,d); 205 Intrepid::FunctionSpaceTools::scalarMultiplyDataField<Real>(KgradN, 206 *valZ_eval, 207 *(fe_vol_->gradN())); 208 // Integrate stiffness term 209 Intrepid::FunctionSpaceTools::integrate<Real>(*jac, 210 KgradN, 211 *(fe_vol_->gradNdetJ()), 212 Intrepid::COMP_CPP, false); 213 214 // Apply Dirichlet boundary conditions 215 int numSideSets = bdryCellLocIds_.size(); 216 if (numSideSets > 0) { 217 for (int i = 0; i < numSideSets; ++i) { 218 if ((i==2)) { 219 int numLocalSideIds = bdryCellLocIds_[i].size(); 220 for (int j = 0; j < numLocalSideIds; ++j) { 221 int numCellsSide = bdryCellLocIds_[i][j].size(); 222 int numBdryDofs = fidx_[j].size(); 223 for (int k = 0; k < numCellsSide; ++k) { 224 int cidx = bdryCellLocIds_[i][j][k]; 225 for (int l = 0; l < numBdryDofs; ++l) { 226 for (int m = 0; m < f; ++m) { 227 (*jac)(cidx,fidx_[j][l],m) = static_cast<Real>(0); 228 } 229 (*jac)(cidx,fidx_[j][l],fidx_[j][l]) = static_cast<Real>(1); 230 } 231 } 232 } 233 } 234 } 235 } 236 } 237 Jacobian_2(ROL::Ptr<Intrepid::FieldContainer<Real>> & jac,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)238 void Jacobian_2(ROL::Ptr<Intrepid::FieldContainer<Real> > & jac, 239 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 240 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 241 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 242 const Real one(1); 243 // Get dimensions 244 int c = fe_vol_->gradN()->dimension(0); 245 int f = fe_vol_->gradN()->dimension(1); 246 int p = fe_vol_->gradN()->dimension(2); 247 int d = fe_vol_->gradN()->dimension(3); 248 // Initialize Jacobian storage 249 jac = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, f); 250 // Build density-dependent conductivity function 251 ROL::Ptr<Intrepid::FieldContainer<Real> > valZ_eval = 252 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p); 253 fe_vol_->evaluateValue(valZ_eval, z_coeff); 254 for (int i = 0; i < c; ++i) { 255 for (int j = 0; j < p; ++j) { 256 (*valZ_eval)(i,j) = SIMPpower_*(one - minConductivity_) 257 * std::pow((*valZ_eval)(i,j),SIMPpower_-one); 258 } 259 } 260 // Build derivative of conductivity function 261 Intrepid::FieldContainer<Real> dKN(c,f,p); 262 Intrepid::FunctionSpaceTools::scalarMultiplyDataField<Real>(dKN, 263 *valZ_eval, 264 *(fe_vol_->N())); 265 // Integrate stiffness term 266 ROL::Ptr<Intrepid::FieldContainer<Real> > gradU_eval = 267 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p, d); 268 fe_vol_->evaluateGradient(gradU_eval, u_coeff); 269 Intrepid::FieldContainer<Real> gradUgradN(c,f,p); 270 Intrepid::FunctionSpaceTools::dotMultiplyDataField<Real>(gradUgradN, 271 *gradU_eval, 272 (*fe_vol_->gradNdetJ())); 273 Intrepid::FunctionSpaceTools::integrate<Real>(*jac, 274 gradUgradN, 275 dKN, 276 Intrepid::COMP_CPP, false); 277 278 // Apply Dirichlet boundary conditions 279 int numSideSets = bdryCellLocIds_.size(); 280 if (numSideSets > 0) { 281 for (int i = 0; i < numSideSets; ++i) { 282 if ((i==2)) { 283 int numLocalSideIds = bdryCellLocIds_[i].size(); 284 for (int j = 0; j < numLocalSideIds; ++j) { 285 int numCellsSide = bdryCellLocIds_[i][j].size(); 286 int numBdryDofs = fidx_[j].size(); 287 for (int k = 0; k < numCellsSide; ++k) { 288 int cidx = bdryCellLocIds_[i][j][k]; 289 for (int l = 0; l < numBdryDofs; ++l) { 290 for (int m=0; m < f; ++m) { 291 (*jac)(cidx,fidx_[j][l],m) = static_cast<Real>(0); 292 } 293 } 294 } 295 } 296 } 297 } 298 } 299 } 300 Hessian_11(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)301 void Hessian_11(ROL::Ptr<Intrepid::FieldContainer<Real> > & hess, 302 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & l_coeff, 303 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 304 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 305 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 306 throw Exception::Zero(">>> (PDE_Poisson_TopOpt::Hessian_11): Zero Hessian."); 307 } 308 Hessian_12(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)309 void Hessian_12(ROL::Ptr<Intrepid::FieldContainer<Real> > & hess, 310 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & l_coeff, 311 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 312 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 313 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 314 const Real one(1); 315 // Get dimensions 316 int c = fe_vol_->gradN()->dimension(0); 317 int f = fe_vol_->gradN()->dimension(1); 318 int p = fe_vol_->gradN()->dimension(2); 319 int d = fe_vol_->gradN()->dimension(3); 320 // Initialize Hessian storage 321 hess = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, f); 322 // Build density-dependent conductivity function 323 ROL::Ptr<Intrepid::FieldContainer<Real> > valZ_eval = 324 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p); 325 fe_vol_->evaluateValue(valZ_eval, z_coeff); 326 for (int i = 0; i < c; ++i) { 327 for (int j = 0; j < p; ++j) { 328 (*valZ_eval)(i,j) = SIMPpower_*(one - minConductivity_) 329 * std::pow((*valZ_eval)(i,j),SIMPpower_-one); 330 } 331 } 332 // Apply Dirichlet conditions to the multipliers. 333 ROL::Ptr<Intrepid::FieldContainer<Real> > l0_coeff 334 = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f); 335 *l0_coeff = *l_coeff; 336 int numSideSets = bdryCellLocIds_.size(); 337 if (numSideSets > 0) { 338 for (int i = 0; i < numSideSets; ++i) { 339 if ((i==2)) { 340 int numLocalSideIds = bdryCellLocIds_[i].size(); 341 for (int j = 0; j < numLocalSideIds; ++j) { 342 int numCellsSide = bdryCellLocIds_[i][j].size(); 343 int numBdryDofs = fidx_[j].size(); 344 for (int k = 0; k < numCellsSide; ++k) { 345 int cidx = bdryCellLocIds_[i][j][k]; 346 for (int l = 0; l < numBdryDofs; ++l) { 347 (*l0_coeff)(cidx,fidx_[j][l]) = static_cast<Real>(0); 348 } 349 } 350 } 351 } 352 } 353 } 354 // Build derivative of conductivity function 355 Intrepid::FieldContainer<Real> dKN(c,f,p); 356 Intrepid::FunctionSpaceTools::scalarMultiplyDataField<Real>(dKN, 357 *valZ_eval, 358 *(fe_vol_->N())); 359 // Integrate stiffness term 360 ROL::Ptr<Intrepid::FieldContainer<Real> > gradL_eval = 361 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p, d); 362 fe_vol_->evaluateGradient(gradL_eval, l0_coeff); 363 Intrepid::FieldContainer<Real> gradLgradN(c,f,p); 364 Intrepid::FunctionSpaceTools::dotMultiplyDataField<Real>(gradLgradN, 365 *gradL_eval, 366 (*fe_vol_->gradNdetJ())); 367 Intrepid::FunctionSpaceTools::integrate<Real>(*hess, 368 dKN, 369 gradLgradN, 370 Intrepid::COMP_CPP, false); 371 } 372 Hessian_21(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)373 void Hessian_21(ROL::Ptr<Intrepid::FieldContainer<Real> > & hess, 374 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & l_coeff, 375 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 376 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 377 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 378 const Real one(1); 379 // Get dimensions 380 int c = fe_vol_->gradN()->dimension(0); 381 int f = fe_vol_->gradN()->dimension(1); 382 int p = fe_vol_->gradN()->dimension(2); 383 int d = fe_vol_->gradN()->dimension(3); 384 // Initialize Hessian storage 385 hess = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, f); 386 // Build density-dependent conductivity function 387 ROL::Ptr<Intrepid::FieldContainer<Real> > valZ_eval = 388 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p); 389 fe_vol_->evaluateValue(valZ_eval, z_coeff); 390 for (int i = 0; i < c; ++i) { 391 for (int j = 0; j < p; ++j) { 392 (*valZ_eval)(i,j) = SIMPpower_*(one - minConductivity_) 393 * std::pow((*valZ_eval)(i,j),SIMPpower_-one); 394 } 395 } 396 // Apply Dirichlet conditions to the multipliers. 397 ROL::Ptr<Intrepid::FieldContainer<Real> > l0_coeff 398 = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f); 399 *l0_coeff = *l_coeff; 400 int numSideSets = bdryCellLocIds_.size(); 401 if (numSideSets > 0) { 402 for (int i = 0; i < numSideSets; ++i) { 403 if ((i==2)) { 404 int numLocalSideIds = bdryCellLocIds_[i].size(); 405 for (int j = 0; j < numLocalSideIds; ++j) { 406 int numCellsSide = bdryCellLocIds_[i][j].size(); 407 int numBdryDofs = fidx_[j].size(); 408 for (int k = 0; k < numCellsSide; ++k) { 409 int cidx = bdryCellLocIds_[i][j][k]; 410 for (int l = 0; l < numBdryDofs; ++l) { 411 (*l0_coeff)(cidx,fidx_[j][l]) = static_cast<Real>(0); 412 } 413 } 414 } 415 } 416 } 417 } 418 // Build derivative of conductivity function 419 Intrepid::FieldContainer<Real> dKN(c,f,p); 420 Intrepid::FunctionSpaceTools::scalarMultiplyDataField<Real>(dKN, 421 *valZ_eval, 422 *(fe_vol_->N())); 423 // Integrate stiffness term 424 ROL::Ptr<Intrepid::FieldContainer<Real> > gradL_eval = 425 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p, d); 426 fe_vol_->evaluateGradient(gradL_eval, l0_coeff); 427 Intrepid::FieldContainer<Real> gradLgradN(c,f,p); 428 Intrepid::FunctionSpaceTools::dotMultiplyDataField<Real>(gradLgradN, 429 *gradL_eval, 430 (*fe_vol_->gradNdetJ())); 431 Intrepid::FunctionSpaceTools::integrate<Real>(*hess, 432 gradLgradN, 433 dKN, 434 Intrepid::COMP_CPP, false); 435 } 436 Hessian_22(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)437 void Hessian_22(ROL::Ptr<Intrepid::FieldContainer<Real> > & hess, 438 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & l_coeff, 439 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 440 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 441 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 442 const Real one(1), two(2); 443 // Get dimensions 444 int c = fe_vol_->gradN()->dimension(0); 445 int f = fe_vol_->gradN()->dimension(1); 446 int p = fe_vol_->gradN()->dimension(2); 447 int d = fe_vol_->gradN()->dimension(3); 448 // Initialize Hessian storage 449 hess = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, f); 450 // Build density-dependent conductivity function 451 ROL::Ptr<Intrepid::FieldContainer<Real> > valZ_eval = 452 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p); 453 fe_vol_->evaluateValue(valZ_eval, z_coeff); 454 for (int i = 0; i < c; ++i) { 455 for (int j = 0; j < p; ++j) { 456 (*valZ_eval)(i,j) = SIMPpower_ * (SIMPpower_ - one)*(one - minConductivity_) 457 * std::pow((*valZ_eval)(i,j),SIMPpower_-two); 458 } 459 } 460 // Apply Dirichlet conditions to the multipliers. 461 ROL::Ptr<Intrepid::FieldContainer<Real> > l0_coeff 462 = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f); 463 *l0_coeff = *l_coeff; 464 int numSideSets = bdryCellLocIds_.size(); 465 if (numSideSets > 0) { 466 for (int i = 0; i < numSideSets; ++i) { 467 if ((i==2)) { 468 int numLocalSideIds = bdryCellLocIds_[i].size(); 469 for (int j = 0; j < numLocalSideIds; ++j) { 470 int numCellsSide = bdryCellLocIds_[i][j].size(); 471 int numBdryDofs = fidx_[j].size(); 472 for (int k = 0; k < numCellsSide; ++k) { 473 int cidx = bdryCellLocIds_[i][j][k]; 474 for (int l = 0; l < numBdryDofs; ++l) { 475 (*l0_coeff)(cidx,fidx_[j][l]) = static_cast<Real>(0); 476 } 477 } 478 } 479 } 480 } 481 } 482 // Build derivative of conductivity function 483 Intrepid::FieldContainer<Real> dKN(c,f,p); 484 Intrepid::FunctionSpaceTools::scalarMultiplyDataField<Real>(dKN, 485 *valZ_eval, 486 *(fe_vol_->N())); 487 // Integrate stiffness term 488 ROL::Ptr<Intrepid::FieldContainer<Real> > gradU_eval = 489 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p, d); 490 fe_vol_->evaluateGradient(gradU_eval, u_coeff); 491 ROL::Ptr<Intrepid::FieldContainer<Real> > gradL_eval = 492 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p, d); 493 fe_vol_->evaluateGradient(gradL_eval, l0_coeff); 494 Intrepid::FieldContainer<Real> gradUgradL(c,p); 495 Intrepid::FunctionSpaceTools::dotMultiplyDataData<Real>(gradUgradL, 496 *gradU_eval, 497 *gradL_eval); 498 Intrepid::FieldContainer<Real> NgradUgradL(c,f,p); 499 Intrepid::FunctionSpaceTools::scalarMultiplyDataField<Real>(NgradUgradL, 500 gradUgradL, 501 *(fe_vol_->NdetJ())); 502 Intrepid::FunctionSpaceTools::integrate<Real>(*hess, 503 dKN, 504 NgradUgradL, 505 Intrepid::COMP_CPP, false); 506 } 507 RieszMap_1(ROL::Ptr<Intrepid::FieldContainer<Real>> & riesz)508 void RieszMap_1(ROL::Ptr<Intrepid::FieldContainer<Real> > & riesz) { 509 // GET DIMENSIONS 510 int c = fe_vol_->N()->dimension(0); 511 int f = fe_vol_->N()->dimension(1); 512 // INITIALIZE RIESZ 513 riesz = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, f); 514 *riesz = *fe_vol_->stiffMat(); 515 Intrepid::RealSpaceTools<Real>::add(*riesz,*(fe_vol_->massMat())); 516 } 517 RieszMap_2(ROL::Ptr<Intrepid::FieldContainer<Real>> & riesz)518 void RieszMap_2(ROL::Ptr<Intrepid::FieldContainer<Real> > & riesz) { 519 // GET DIMENSIONS 520 int c = fe_vol_->N()->dimension(0); 521 int f = fe_vol_->N()->dimension(1); 522 // INITIALIZE RIESZ 523 riesz = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, f); 524 *riesz = *fe_vol_->massMat(); 525 } 526 getFields()527 std::vector<ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > > getFields() { 528 return basisPtrs_; 529 } 530 setCellNodes(const ROL::Ptr<Intrepid::FieldContainer<Real>> & volCellNodes,const std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>>> & bdryCellNodes,const std::vector<std::vector<std::vector<int>>> & bdryCellLocIds)531 void setCellNodes(const ROL::Ptr<Intrepid::FieldContainer<Real> > &volCellNodes, 532 const std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real> > > > &bdryCellNodes, 533 const std::vector<std::vector<std::vector<int> > > &bdryCellLocIds) { 534 volCellNodes_ = volCellNodes; 535 bdryCellNodes_ = bdryCellNodes; 536 bdryCellLocIds_ = bdryCellLocIds; 537 // Finite element definition. 538 fe_vol_ = ROL::makePtr<FE<Real>>(volCellNodes_,basisPtr_,cellCub_); 539 fidx_ = fe_vol_->getBoundaryDofs(); 540 computeForce(); 541 computeDirichlet(); 542 } 543 computeForce(void)544 void computeForce(void) { 545 int c = fe_vol_->cubPts()->dimension(0); 546 int p = fe_vol_->cubPts()->dimension(1); 547 int d = fe_vol_->cubPts()->dimension(2); 548 std::vector<Real> pt(d,0); 549 force_eval_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(c,p); 550 nforce_eval_ = ROL::makePtr<Intrepid::FieldContainer<Real>>(c,p); 551 for (int i = 0; i < c; ++i) { 552 for (int j = 0; j < p; ++j) { 553 for (int k = 0; k < d; ++k) { 554 pt[k] = (*fe_vol_->cubPts())(i,j,k); 555 } 556 (*force_eval_)(i,j) = ForceFunc(pt); 557 (*nforce_eval_)(i,j) = -(*force_eval_)(i,j); 558 } 559 } 560 } 561 computeDirichlet(void)562 void computeDirichlet(void) { 563 // Compute Dirichlet values at DOFs. 564 int d = basisPtr_->getBaseCellTopology().getDimension(); 565 int numSidesets = bdryCellLocIds_.size(); 566 bdryCellDofValues_.resize(numSidesets); 567 for (int i=0; i<numSidesets; ++i) { 568 int numLocSides = bdryCellLocIds_[i].size(); 569 bdryCellDofValues_[i].resize(numLocSides); 570 for (int j=0; j<numLocSides; ++j) { 571 int c = bdryCellLocIds_[i][j].size(); 572 int f = basisPtr_->getCardinality(); 573 bdryCellDofValues_[i][j] = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f); 574 ROL::Ptr<Intrepid::FieldContainer<Real> > coords = 575 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, d); 576 if (c > 0) { 577 fe_vol_->computeDofCoords(coords, bdryCellNodes_[i][j]); 578 } 579 for (int k=0; k<c; ++k) { 580 for (int l=0; l<f; ++l) { 581 std::vector<Real> dofpoint(d); 582 for (int m=0; m<d; ++m) { 583 dofpoint[m] = (*coords)(k, l, m); 584 } 585 (*bdryCellDofValues_[i][j])(k, l) = dirichletFunc(dofpoint, i, j); 586 } 587 } 588 } 589 } 590 } 591 getFE(void) const592 const ROL::Ptr<FE<Real> > getFE(void) const { 593 return fe_vol_; 594 } 595 getForce(void) const596 const ROL::Ptr<Intrepid::FieldContainer<Real> > getForce(void) const { 597 return force_eval_; 598 } 599 600 }; // PDE_Poisson 601 602 603 604 template <class Real> 605 class PDE_Filter : public PDE<Real> { 606 private: 607 // Finite element basis information 608 ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > basisPtr_; 609 std::vector<ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > > basisPtrs_; 610 // Cell cubature information 611 ROL::Ptr<Intrepid::Cubature<Real> > cellCub_; 612 // Cell node information 613 ROL::Ptr<Intrepid::FieldContainer<Real> > volCellNodes_; 614 // Finite element definition 615 ROL::Ptr<FE<Real> > fe_; 616 // Problem parameters. 617 Real lengthScale_; 618 619 public: PDE_Filter(Teuchos::ParameterList & parlist)620 PDE_Filter(Teuchos::ParameterList &parlist) { 621 // Finite element fields. 622 basisPtr_ = ROL::makePtr<Intrepid::Basis_HGRAD_QUAD_C1_FEM<Real, Intrepid::FieldContainer<Real> >>(); 623 // Quadrature rules. 624 shards::CellTopology cellType = basisPtr_->getBaseCellTopology(); // get the cell type from the basis 625 Intrepid::DefaultCubatureFactory<Real> cubFactory; // create cubature factory 626 int cubDegree = parlist.sublist("Problem").get("Cubature Degree", 2); // set cubature degree, e.g., 2 627 cellCub_ = cubFactory.create(cellType, cubDegree); // create default cubature 628 629 basisPtrs_.clear(); 630 basisPtrs_.push_back(basisPtr_); // Filter components; there is only one, but we need d because of the infrastructure. 631 632 // Other problem parameters. 633 Real filterRadius = parlist.sublist("Problem").get("Filter Radius", 0.1); 634 lengthScale_ = std::pow(filterRadius, 2)/12.0; 635 } 636 residual(ROL::Ptr<Intrepid::FieldContainer<Real>> & res,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)637 void residual(ROL::Ptr<Intrepid::FieldContainer<Real> > & res, 638 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 639 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 640 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 641 // Retrieve dimensions. 642 int c = fe_->gradN()->dimension(0); 643 int f = fe_->gradN()->dimension(1); 644 int p = fe_->gradN()->dimension(2); 645 int d = fe_->gradN()->dimension(3); 646 647 // Initialize residuals. 648 res = ROL::makePtr<Intrepid::FieldContainer<Real>>(c,f); 649 650 // Evaluate/interpolate finite element fields on cells. 651 ROL::Ptr<Intrepid::FieldContainer<Real> > valU_eval 652 = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p); 653 ROL::Ptr<Intrepid::FieldContainer<Real> > gradU_eval 654 = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p, d); 655 ROL::Ptr<Intrepid::FieldContainer<Real> > valZ_eval 656 = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p); 657 fe_->evaluateValue(valU_eval, u_coeff); 658 fe_->evaluateValue(valZ_eval, z_coeff); 659 fe_->evaluateGradient(gradU_eval, u_coeff); 660 661 Intrepid::RealSpaceTools<Real>::scale(*gradU_eval, lengthScale_); 662 Intrepid::RealSpaceTools<Real>::scale(*valZ_eval, static_cast<Real>(-1)); 663 664 /*** Evaluate weak form of the residual. ***/ 665 Intrepid::FunctionSpaceTools::integrate<Real>(*res, 666 *gradU_eval, // R*gradU 667 *fe_->gradNdetJ(), // gradN 668 Intrepid::COMP_CPP, 669 false); 670 Intrepid::FunctionSpaceTools::integrate<Real>(*res, 671 *valU_eval, // U 672 *fe_->NdetJ(), // N 673 Intrepid::COMP_CPP, 674 true); 675 Intrepid::FunctionSpaceTools::integrate<Real>(*res, 676 *valZ_eval, // -Z 677 *fe_->NdetJ(), // N 678 Intrepid::COMP_CPP, 679 true); 680 } 681 Jacobian_1(ROL::Ptr<Intrepid::FieldContainer<Real>> & jac,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)682 void Jacobian_1(ROL::Ptr<Intrepid::FieldContainer<Real> > & jac, 683 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 684 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 685 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 686 // Retrieve dimensions. 687 int c = fe_->gradN()->dimension(0); 688 int f = fe_->gradN()->dimension(1); 689 690 // Initialize Jacobians. 691 jac = ROL::makePtr<Intrepid::FieldContainer<Real>>(c,f,f); 692 693 /*** Evaluate weak form of the Jacobian. ***/ 694 *jac = *(fe_->stiffMat()); 695 Intrepid::RealSpaceTools<Real>::scale(*jac, lengthScale_); // ls*gradN1 . gradN2 696 Intrepid::RealSpaceTools<Real>::add(*jac,*(fe_->massMat())); // + N1 * N2 697 } 698 699 Jacobian_2(ROL::Ptr<Intrepid::FieldContainer<Real>> & jac,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)700 void Jacobian_2(ROL::Ptr<Intrepid::FieldContainer<Real> > & jac, 701 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 702 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 703 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 704 // Retrieve dimensions. 705 int c = fe_->gradN()->dimension(0); 706 int f = fe_->gradN()->dimension(1); 707 708 // Initialize Jacobians. 709 jac = ROL::makePtr<Intrepid::FieldContainer<Real>>(c,f,f); 710 711 /*** Evaluate weak form of the Jacobian. ***/ 712 *jac = *(fe_->massMat()); 713 Intrepid::RealSpaceTools<Real>::scale(*jac, static_cast<Real>(-1)); // -N1 * N2 714 } 715 Hessian_11(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)716 void Hessian_11(ROL::Ptr<Intrepid::FieldContainer<Real> > & hess, 717 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & l_coeff, 718 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 719 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 720 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 721 throw Exception::Zero(">>> (PDE_Filter::Hessian_11): Hessian is zero."); 722 } 723 Hessian_12(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)724 void Hessian_12(ROL::Ptr<Intrepid::FieldContainer<Real> > & hess, 725 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & l_coeff, 726 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 727 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 728 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 729 throw Exception::Zero(">>> (PDE_Filter::Hessian_12): Hessian is zero."); 730 } 731 Hessian_21(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)732 void Hessian_21(ROL::Ptr<Intrepid::FieldContainer<Real> > & hess, 733 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & l_coeff, 734 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 735 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 736 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 737 throw Exception::Zero(">>> (PDE_Filter::Hessian_21): Hessian is zero."); 738 } 739 Hessian_22(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & u_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)740 void Hessian_22(ROL::Ptr<Intrepid::FieldContainer<Real> > & hess, 741 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & l_coeff, 742 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & u_coeff, 743 const ROL::Ptr<const Intrepid::FieldContainer<Real> > & z_coeff = ROL::nullPtr, 744 const ROL::Ptr<const std::vector<Real> > & z_param = ROL::nullPtr) { 745 throw Exception::Zero(">>> (PDE_Filter::Hessian_22): Hessian is zero."); 746 } 747 RieszMap_1(ROL::Ptr<Intrepid::FieldContainer<Real>> & riesz)748 void RieszMap_1(ROL::Ptr<Intrepid::FieldContainer<Real> > & riesz) { 749 //throw Exception::NotImplemented(">>> (PDE_Filter::RieszMap_1): Not implemented."); 750 // Retrieve dimensions. 751 int c = fe_->gradN()->dimension(0); 752 int f = fe_->gradN()->dimension(1); 753 754 // Initialize Riesz map. 755 riesz = ROL::makePtr<Intrepid::FieldContainer<Real>>(c,f,f); 756 757 *riesz = *(fe_->stiffMat()); 758 Intrepid::RealSpaceTools<Real>::add(*riesz,*(fe_->massMat())); 759 } 760 RieszMap_2(ROL::Ptr<Intrepid::FieldContainer<Real>> & riesz)761 void RieszMap_2(ROL::Ptr<Intrepid::FieldContainer<Real> > & riesz) { 762 //throw Exception::NotImplemented(">>> (PDE_Filter::RieszMap_2): Not implemented."); 763 int c = fe_->gradN()->dimension(0); 764 int f = fe_->gradN()->dimension(1); 765 766 // Initialize Riesz map. 767 riesz = ROL::makePtr<Intrepid::FieldContainer<Real>>(c,f,f); 768 769 *riesz = *(fe_->massMat()); 770 } 771 getFields()772 std::vector<ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real> > > > getFields() { 773 return basisPtrs_; 774 } 775 setCellNodes(const ROL::Ptr<Intrepid::FieldContainer<Real>> & volCellNodes,const std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>>> & bdryCellNodes,const std::vector<std::vector<std::vector<int>>> & bdryCellLocIds)776 void setCellNodes(const ROL::Ptr<Intrepid::FieldContainer<Real> > &volCellNodes, 777 const std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real> > > > &bdryCellNodes, 778 const std::vector<std::vector<std::vector<int> > > &bdryCellLocIds) { 779 volCellNodes_ = volCellNodes; 780 // Finite element definition. 781 fe_ = ROL::makePtr<FE<Real>>(volCellNodes_,basisPtr_,cellCub_); 782 } 783 setFieldPattern(const std::vector<std::vector<int>> & fieldPattern)784 void setFieldPattern(const std::vector<std::vector<int> > & fieldPattern) {} 785 786 }; // PDE_Filter 787 788 #endif 789