1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 2 // vi: set et ts=4 sw=2 sts=2: 3 #ifndef DUNE_REFINED_SIMPLEX_LOCALBASIS_HH 4 #define DUNE_REFINED_SIMPLEX_LOCALBASIS_HH 5 6 /** \file 7 \brief Contains a base class for LocalBasis classes based on uniform refinement 8 */ 9 10 #include <dune/common/fvector.hh> 11 #include <dune/common/exceptions.hh> 12 #include <dune/localfunctions/common/localbasis.hh> 13 14 namespace Dune 15 { 16 template<class D, int dim> 17 class RefinedSimplexLocalBasis 18 { 19 protected: RefinedSimplexLocalBasis()20 RefinedSimplexLocalBasis() 21 { 22 DUNE_THROW(Dune::NotImplemented,"RefinedSimplexLocalBasis not implemented for dim > 3."); 23 } 24 }; 25 26 /**@ingroup LocalBasisImplementation 27 \brief Base class for LocalBasis classes based on uniform refinement in 1D; provides numbering and local coordinates of subelements 28 29 \tparam D Type to represent the field in the domain. 30 31 \nosubgrouping 32 */ 33 template<class D> 34 class RefinedSimplexLocalBasis<D,1> 35 { 36 protected: 37 38 /** \brief Protected default constructor so this class can only be instantiated as a base class. */ RefinedSimplexLocalBasis()39 RefinedSimplexLocalBasis() {} 40 41 /** \brief Get the number of the subelement containing a given point. 42 * 43 * The subelements are ordered according to 44 * 45 * 0 1 46 * |-------:-------| 47 * 48 * \param[in] global Coordinates in the reference element 49 * \returns Number of the subtriangle containing <tt>global</tt> 50 */ getSubElement(const FieldVector<D,1> & global)51 static int getSubElement(const FieldVector<D,1>& global) 52 { 53 if (global[0] <= 0.5) 54 return 0; 55 else if (global[0] <= 1.0) 56 return 1; 57 58 DUNE_THROW(InvalidStateException, "no subelement defined"); 59 } 60 61 /** \brief Get local coordinates in the subelement 62 63 \param[in] global Coordinates in the reference element 64 \param[out] subElement Number of the subelement containing <tt>global</tt> 65 \param[out] local The local coordinates in the subelement 66 */ getSubElement(const FieldVector<D,1> & global,int & subElement,FieldVector<D,1> & local)67 static void getSubElement(const FieldVector<D,1>& global, 68 int& subElement, 69 FieldVector<D,1>& local) 70 { 71 if (global[0] <= 0.5) { 72 subElement = 0; 73 local[0] = 2.0 * global[0]; 74 return; 75 } 76 77 subElement = 1; 78 local[0] = 2.0 * global[0] - 1.0; 79 } 80 81 }; 82 83 84 /**@ingroup LocalBasisImplementation 85 \brief Base class for LocalBasis classes based on uniform refinement in 2D; provides numbering and local coordinates of subelements 86 87 Shape functions like these are necessary for hierarchical error estimators 88 for certain nonlinear problems. 89 90 \tparam D Type to represent the field in the domain. 91 92 \nosubgrouping 93 */ 94 template<class D> 95 class RefinedSimplexLocalBasis<D,2> 96 { 97 protected: 98 99 /** \brief Protected default constructor so this class can only be instantiated as a base class. */ RefinedSimplexLocalBasis()100 RefinedSimplexLocalBasis() {} 101 102 /** \brief Get the number of the subtriangle containing a given point. 103 * 104 * The triangles are ordered according to 105 * \verbatim 106 |\ 107 |2\ 108 |--\ 109 |\3|\ 110 |0\|1\ 111 ------ 112 \endverbatim 113 * 114 * \param[in] global Coordinates in the reference triangle 115 * \returns Number of the subtriangle containing <tt>global</tt> 116 */ getSubElement(const FieldVector<D,2> & global)117 static int getSubElement(const FieldVector<D,2>& global) 118 { 119 if (global[0] + global[1] <= 0.5) 120 return 0; 121 else if (global[0] >= 0.5) 122 return 1; 123 else if (global[1] >= 0.5) 124 return 2; 125 126 return 3; 127 } 128 129 /** \brief Get local coordinates in the subtriangle 130 131 \param[in] global Coordinates in the reference triangle 132 \param[out] subElement Number of the subtriangle containing <tt>global</tt> 133 \param[out] local The local coordinates in the subtriangle 134 */ getSubElement(const FieldVector<D,2> & global,int & subElement,FieldVector<D,2> & local)135 static void getSubElement(const FieldVector<D,2>& global, 136 int& subElement, 137 FieldVector<D,2>& local) 138 { 139 if (global[0] + global[1] <= 0.5) { 140 subElement = 0; 141 local[0] = 2*global[0]; 142 local[1] = 2*global[1]; 143 return; 144 } else if (global[0] >= 0.5) { 145 subElement = 1; 146 local[0] = 2*global[0]-1; 147 local[1] = 2*global[1]; 148 return; 149 } else if (global[1] >= 0.5) { 150 subElement = 2; 151 local[0] = 2*global[0]; 152 local[1] = 2*global[1]-1; 153 return; 154 } 155 156 subElement = 3; 157 local[0] = -2 * global[0] + 1; 158 local[1] = -2 * global[1] + 1; 159 160 } 161 162 163 }; 164 165 /**@ingroup LocalBasisImplementation 166 \brief Base class for LocalBasis classes based on uniform refinement in 3D; provides numbering and local coordinates of subelements 167 168 Shape functions like these are necessary for hierarchical error estimators 169 for certain nonlinear problems. 170 171 \tparam D Type to represent the field in the domain. 172 173 \nosubgrouping 174 */ 175 template<class D> 176 class RefinedSimplexLocalBasis<D,3> 177 { 178 protected: 179 180 /** \brief Protected default constructor so this class can only be instantiated as a base class. */ RefinedSimplexLocalBasis()181 RefinedSimplexLocalBasis() {} 182 183 /** \brief Get the number of the subsimplex containing a given point in the reference element 184 * 185 * Defining the following points in the reference simplex 186 * 187 * 0: (0.0, 0.0, 0.0) 188 * 1: (1.0, 0.0, 0.0) 189 * 2: (0.0, 1.0, 0.0) 190 * 3: (0.0, 0.0, 1.0) 191 * 4: (0.5, 0.0, 0.0) 192 * 5: (0.5, 0.5, 0.0) 193 * 6: (0.0, 0.5, 0.0) 194 * 7: (0.0, 0.0, 0.5) 195 * 8: (0.5, 0.0, 0.5) 196 * 9: (0.0, 0.5, 0.5) 197 * 198 * The subsimplices are numbered according to 199 * 200 * 0: 0467 - 201 * 1: 4158 |_ "cut off" vertices 202 * 2: 6529 | 203 * 3: 7893 - 204 * 205 * 4: 6487 - 206 * 5: 4568 |_ octahedron partition 207 * 6: 6897 | 208 * 7: 6895 - 209 * 210 * \param[in] global Coordinates in the reference simplex 211 * \returns Number of the subsimplex containing <tt>global</tt> 212 */ getSubElement(const FieldVector<D,3> & global)213 static int getSubElement(const FieldVector<D,3>& global) 214 { 215 if (global[0] + global[1] + global[2] <= 0.5) 216 return 0; 217 else if (global[0] >= 0.5) 218 return 1; 219 else if (global[1] >= 0.5) 220 return 2; 221 else if (global[2] >= 0.5) 222 return 3; 223 else if ((global[0] + global[1] <= 0.5)and (global[1] + global[2] <= 0.5)) 224 return 4; 225 else if ((global[0] + global[1] >= 0.5)and (global[1] + global[2] <= 0.5)) 226 return 5; 227 else if ((global[0] + global[1] <= 0.5)and (global[1] + global[2] >= 0.5)) 228 return 6; 229 else if ((global[0] + global[1] >= 0.5)and (global[1] + global[2] >= 0.5)) 230 return 7; 231 232 DUNE_THROW(InvalidStateException, "no subelement defined"); 233 234 } 235 /** \brief Get local coordinates in the subsimplex 236 237 \param[in] global Coordinates in the reference simplex 238 \param[out] subElement Number of the subsimplex containing <tt>global</tt> 239 \param[out] local The local coordinates in the subsimplex 240 */ getSubElement(const FieldVector<D,3> & global,int & subElement,FieldVector<D,3> & local)241 static void getSubElement(const FieldVector<D,3>& global, 242 int& subElement, 243 FieldVector<D,3>& local) 244 { 245 if (global[0] + global[1] + global[2] <= 0.5) { 246 subElement = 0; 247 local = global; 248 local *= 2.0; 249 return; 250 } else if (global[0] >= 0.5) { 251 subElement = 1; 252 local = global; 253 local[0] -= 0.5; 254 local *= 2.0; 255 return; 256 } else if (global[1] >= 0.5) { 257 subElement = 2; 258 local = global; 259 local[1] -= 0.5; 260 local *= 2.0; 261 return; 262 } else if (global[2] >= 0.5) { 263 subElement = 3; 264 local = global; 265 local[2] -= 0.5; 266 local *= 2.0; 267 return; 268 } else if ((global[0] + global[1] <= 0.5)and (global[1] + global[2] <= 0.5)) { 269 subElement = 4; 270 local[0] = 2.0 * global[1]; 271 local[1] = 2.0 * (0.5 - global[0] - global[1]); 272 local[2] = 2.0 * (-0.5 + global[0] + global[1] + global[2]); 273 // Dune::FieldMatrix<double,3,3> A(0.0); 274 // A[0][1] = 2.0; 275 // A[1][0] = -2.0; 276 // A[1][1] = -2.0; 277 // A[2][0] = 2.0; 278 // A[2][1] = 2.0; 279 // A[2][2] = 2.0; 280 // A.mv(global,local); 281 // local[1] += 1.0; 282 // local[2] -= 1.0; 283 return; 284 } else if ((global[0] + global[1] >= 0.5)and (global[1] + global[2] <= 0.5)) { 285 subElement = 5; 286 local[0] = 2.0 * (0.5 - global[0]); 287 local[1] = 2.0 * (0.5 - global[1] - global[2]); 288 local[2] = 2.0 * global[2]; 289 // Dune::FieldMatrix<double,3,3> A(0.0); 290 // A[0][0] = -2.0; 291 // A[1][1] = -2.0; 292 // A[1][2] = -2.0; 293 // A[2][2] = 2.0; 294 // A.mv(global,local); 295 // local[0] += 1.0; 296 // local[1] += 1.0; 297 return; 298 } else if ((global[0] + global[1] <= 0.5)and (global[1] + global[2] >= 0.5)) { 299 subElement = 6; 300 local[0] = 2.0 * (0.5 - global[0] - global[1]); 301 local[1] = 2.0 * global[0]; 302 local[2] = 2.0 * (-0.5 + global[1] + global[2]); 303 // Dune::FieldMatrix<double,3,3> A(0.0); 304 // A[0][0] = -2.0; 305 // A[0][1] = -2.0; 306 // A[1][0] = 2.0; 307 // A[2][1] = 2.0; 308 // A[2][2] = 2.0; 309 // A.mv(global,local); 310 // local[0] += 1.0; 311 // local[2] -= 1.0; 312 return; 313 } else if ((global[0] + global[1] >= 0.5)and (global[1] + global[2] >= 0.5)) { 314 subElement = 7; 315 local[0] = 2.0 * (-0.5 + global[1] + global[2]); 316 local[1] = 2.0 * (0.5 - global[1]); 317 local[2] = 2.0 * (-0.5 + global[0] + global[1]); 318 // Dune::FieldMatrix<double,3,3> A(0.0); 319 // A[0][1] = 2.0; 320 // A[0][2] = 2.0; 321 // A[1][1] = -2.0; 322 // A[2][0] = 2.0; 323 // A[2][1] = 2.0; 324 // A.mv(global,local); 325 // local[0] -= 1.0; 326 // local[1] += 1.0; 327 // local[2] -= 1.0; 328 return; 329 } 330 331 DUNE_THROW(InvalidStateException, "no subelement defined"); 332 333 } 334 335 }; 336 337 338 } 339 340 #endif 341