1 // -*- tab-width: 4; indent-tabs-mode: nil -*- 2 #ifndef DUNE_PDELAB_LOCALOPERATOR_DIFFUSIONMIXED_HH 3 #define DUNE_PDELAB_LOCALOPERATOR_DIFFUSIONMIXED_HH 4 5 #include <cstddef> 6 #include<vector> 7 8 #include<dune/common/exceptions.hh> 9 #include<dune/common/fvector.hh> 10 #include<dune/common/fmatrix.hh> 11 12 #include<dune/geometry/type.hh> 13 #include<dune/geometry/referenceelements.hh> 14 15 #include<dune/pdelab/common/quadraturerules.hh> 16 17 #include"defaultimp.hh" 18 #include"pattern.hh" 19 #include"flags.hh" 20 #include "convectiondiffusionparameter.hh" 21 22 namespace Dune { 23 namespace PDELab { 24 25 // a local operator for solving the Poisson equation 26 // div sigma +a_0 u = f in \Omega, 27 // sigma = -K grad u in \Omega, 28 // u = g on \partial\Omega_D 29 // sigma \cdot \nu = j on \partial\Omega_N 30 // with H(div) conforming (mixed) finite elements 31 // param.A : diffusion tensor dependent on position 32 // param.c : Helmholtz term 33 // param.f : grid function type giving f 34 // param.bctype : grid function type selecting boundary condition 35 // param.g : grid function type giving g 36 template<typename PARAM> 37 class DiffusionMixed : public NumericalJacobianApplyVolume<DiffusionMixed<PARAM> >, 38 public NumericalJacobianVolume<DiffusionMixed<PARAM> >, 39 public FullVolumePattern, 40 public LocalOperatorDefaultFlags 41 { 42 43 using BCType = typename ConvectionDiffusionBoundaryConditions::Type; 44 45 public: 46 // pattern assembly flags 47 enum { doPatternVolume = true }; 48 49 // residual assembly flags 50 enum { doAlphaVolume = true }; 51 enum { doLambdaVolume = true }; 52 enum { doLambdaBoundary = true }; 53 DiffusionMixed(const PARAM & param_,int qorder_v_=2,int qorder_p_=1)54 DiffusionMixed ( const PARAM& param_, 55 int qorder_v_=2, 56 int qorder_p_=1 ) 57 : param(param_), 58 qorder_v(qorder_v_), 59 qorder_p(qorder_p_) 60 { 61 } 62 63 // volume integral depending on test and ansatz functions 64 template<typename EG, typename LFSU, typename X, typename LFSV, typename R> alpha_volume(const EG & eg,const LFSU & lfsu,const X & x,const LFSV & lfsv,R & r) const65 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const 66 { 67 // Define types 68 using VelocitySpace = typename LFSU::template Child<0>::Type; 69 using PressureSpace = typename LFSU::template Child<1>::Type; 70 using DF = typename VelocitySpace::Traits::FiniteElementType:: 71 Traits::LocalBasisType::Traits::DomainFieldType; 72 using RF = typename VelocitySpace::Traits::FiniteElementType:: 73 Traits::LocalBasisType::Traits::RangeFieldType; 74 using VelocityJacobianType = typename VelocitySpace::Traits::FiniteElementType:: 75 Traits::LocalBasisType::Traits::JacobianType; 76 using VelocityRangeType = typename VelocitySpace::Traits::FiniteElementType:: 77 Traits::LocalBasisType::Traits::RangeType; 78 using PressureRangeType = typename PressureSpace::Traits::FiniteElementType:: 79 Traits::LocalBasisType::Traits::RangeType; 80 81 // select the two components 82 using namespace Indices; 83 const auto& velocityspace = child(lfsu,_0); 84 const auto& pressurespace = lfsu.template child<1>(); 85 86 // dimensions 87 const int dim = EG::Geometry::mydimension; 88 89 // References to the cell 90 const auto& cell = eg.entity(); 91 92 // Get geometry 93 auto geo = eg.geometry(); 94 95 // evaluate transformation which must be linear 96 Dune::FieldVector<DF,dim> pos; 97 pos=0.0; 98 auto jac = geo.jacobianInverseTransposed(pos); 99 jac.invert(); 100 auto det = geo.integrationElement(pos); 101 102 // evaluate diffusion tensor at cell center, assume it is constant over elements 103 auto ref_el = referenceElement(geo); 104 auto localcenter = ref_el.position(0,0); 105 auto tensor = param.A(cell,localcenter); 106 tensor.invert(); // need iverse for mixed method 107 108 // Initialize vectors outside for loop 109 std::vector<VelocityRangeType> vbasis(velocityspace.size()); 110 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size()); 111 VelocityRangeType sigma; 112 VelocityRangeType Kinvsigma; 113 std::vector<VelocityJacobianType> vjacbasis(velocityspace.size()); 114 std::vector<PressureRangeType> pbasis(pressurespace.size()); 115 std::vector<RF> divergence(velocityspace.size(),0.0); 116 117 118 // \sigma\cdot v term 119 // loop over quadrature points 120 for (const auto& ip : quadratureRule(geo,qorder_v)) 121 { 122 // evaluate shape functions at ip (this is a Galerkin method) 123 velocityspace.finiteElement().localBasis().evaluateFunction(ip.position(),vbasis); 124 125 // transform basis vectors 126 for (std::size_t i=0; i<velocityspace.size(); i++) 127 { 128 vtransformedbasis[i] = 0.0; 129 jac.umtv(vbasis[i],vtransformedbasis[i]); 130 } 131 132 // compute sigma 133 sigma=0.0; 134 for (std::size_t i=0; i<velocityspace.size(); i++) 135 sigma.axpy(x(velocityspace,i),vtransformedbasis[i]); 136 137 // K^{-1} * sigma 138 tensor.mv(sigma,Kinvsigma); 139 140 // integrate (K^{-1}*sigma) * phi_i 141 auto factor = ip.weight() / det; 142 for (std::size_t i=0; i<velocityspace.size(); i++) 143 r.accumulate(velocityspace,i,(Kinvsigma*vtransformedbasis[i])*factor); 144 } 145 146 // u div v term, div sigma q term, a0*u term 147 // loop over quadrature points 148 for (const auto& ip : quadratureRule(geo,qorder_p)) 149 { 150 // evaluate shape functions at ip (this is a Galerkin method) 151 velocityspace.finiteElement().localBasis().evaluateJacobian(ip.position(),vjacbasis); 152 pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis); 153 154 // compute u 155 PressureRangeType u; 156 u=0.0; 157 for (std::size_t i=0; i<pressurespace.size(); i++) 158 u.axpy(x(pressurespace,i),pbasis[i]); 159 160 // evaluate Helmholtz term (reaction term) 161 auto a0value = param.c(cell,ip.position()); 162 163 // integrate a0 * u * q 164 RF factor = ip.weight(); 165 for (std::size_t i=0; i<pressurespace.size(); i++) 166 r.accumulate(pressurespace,i,-a0value*u*pbasis[i]*factor); 167 168 // compute divergence of velocity basis functions on reference element 169 for (std::size_t i=0; i<velocityspace.size(); i++){ 170 divergence[i] = 0; 171 for (int j=0; j<dim; j++) 172 divergence[i] += vjacbasis[i][j][j]; 173 } 174 175 // integrate sigma * phi_i 176 for (std::size_t i=0; i<velocityspace.size(); i++) 177 r.accumulate(velocityspace,i,-u*divergence[i]*factor); 178 179 // compute divergence of sigma 180 RF divergencesigma = 0.0; 181 for (std::size_t i=0; i<velocityspace.size(); i++) 182 divergencesigma += x(velocityspace,i)*divergence[i]; 183 184 // integrate div sigma * q 185 for (std::size_t i=0; i<pressurespace.size(); i++) 186 r.accumulate(pressurespace,i,-divergencesigma*pbasis[i]*factor); 187 } 188 } 189 190 // volume integral depending only on test functions 191 template<typename EG, typename LFSV, typename R> lambda_volume(const EG & eg,const LFSV & lfsv,R & r) const192 void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const 193 { 194 // Define types 195 using PressureSpace = typename LFSV::template Child<1>::Type; 196 using PressureRangeType = typename PressureSpace::Traits::FiniteElementType:: 197 Traits::LocalBasisType::Traits::RangeType; 198 199 // select the pressure component 200 using namespace Indices; 201 const auto& pressurespace = child(lfsv,_1); 202 203 // References to the cell 204 const auto& cell = eg.entity(); 205 206 // Get geometry 207 auto geo = eg.geometry(); 208 209 // Initialize vectors outside for loop 210 std::vector<PressureRangeType> pbasis(pressurespace.size()); 211 212 // loop over quadrature points 213 for (const auto& ip : quadratureRule(geo,qorder_p)) 214 { 215 // evaluate shape functions 216 pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis); 217 218 // evaluate right hand side parameter function 219 auto y = param.f(cell,ip.position()); 220 221 // integrate f 222 auto factor = ip.weight() * geo.integrationElement(ip.position()); 223 for (std::size_t i=0; i<pressurespace.size(); i++) 224 r.accumulate(pressurespace,i,y*pbasis[i]*factor); 225 } 226 } 227 228 // boundary integral independen of ansatz functions 229 template<typename IG, typename LFSV, typename R> lambda_boundary(const IG & ig,const LFSV & lfsv,R & r) const230 void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const 231 { 232 // Define types 233 using VelocitySpace = typename LFSV::template Child<0>::Type; 234 using DF = typename VelocitySpace::Traits::FiniteElementType:: 235 Traits::LocalBasisType::Traits::DomainFieldType; 236 using VelocityRangeType = typename VelocitySpace::Traits::FiniteElementType:: 237 Traits::LocalBasisType::Traits::RangeType; 238 239 // select the two velocity component 240 using namespace Indices; 241 const auto& velocityspace = child(lfsv,_0); 242 243 // dimensions 244 const int dim = IG::Entity::dimension; 245 246 // References to the inside cell 247 const auto& cell_inside = ig.inside(); 248 249 // Get geometry 250 auto geo = ig.geometry(); 251 auto geo_inside = cell_inside.geometry(); 252 253 // Get geometry of intersection in local coordinates of cell_inside 254 auto geo_in_inside = ig.geometryInInside(); 255 256 // evaluate transformation which must be linear 257 Dune::FieldVector<DF,dim> pos; 258 pos = 0.0; 259 typename IG::Entity::Geometry::JacobianInverseTransposed jac; 260 jac = geo_inside.jacobianInverseTransposed(pos); 261 jac.invert(); 262 auto det = geo_inside.integrationElement(pos); 263 264 // Declare vectors outside for loop 265 std::vector<VelocityRangeType> vbasis(velocityspace.size()); 266 std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size()); 267 268 // loop over quadrature points and integrate normal flux 269 for (const auto& ip : quadratureRule(geo,qorder_v)) 270 { 271 // evaluate boundary condition type 272 auto bctype = param.bctype(ig.intersection(),ip.position()); 273 274 // skip rest if we are on Neumann boundary 275 if (bctype == ConvectionDiffusionBoundaryConditions::Neumann) 276 continue; 277 278 // position of quadrature point in local coordinates of element 279 auto local = geo_in_inside.global(ip.position()); 280 281 // evaluate test shape functions 282 velocityspace.finiteElement().localBasis().evaluateFunction(local,vbasis); 283 284 // transform basis vectors 285 for (std::size_t i=0; i<velocityspace.size(); i++) 286 { 287 vtransformedbasis[i] = 0.0; 288 jac.umtv(vbasis[i],vtransformedbasis[i]); 289 } 290 291 // evaluate Dirichlet boundary condition 292 auto y = param.g(cell_inside,local); 293 294 // integrate g v*normal 295 auto factor = ip.weight()*geo.integrationElement(ip.position())/det; 296 for (std::size_t i=0; i<velocityspace.size(); i++) 297 r.accumulate(velocityspace,i,y*(vtransformedbasis[i]*ig.unitOuterNormal(ip.position()))*factor); 298 } 299 } 300 301 private: 302 const PARAM& param; 303 int qorder_v; 304 int qorder_p; 305 }; 306 307 //! \} group GridFunctionSpace 308 } // namespace PDELab 309 } // namespace Dune 310 311 #endif // DUNE_PDELAB_LOCALOPERATOR_DIFFUSIONMIXED_HH 312