fe_subdivision_2D.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 00019 // Local includes 00020 #include "libmesh/fe.h" 00021 #include "libmesh/libmesh_logging.h" 00022 #include "libmesh/fe_type.h" 00023 #include "libmesh/quadrature.h" 00024 #include "libmesh/face_tri3_subdivision.h" 00025 #include "libmesh/fe_macro.h" 00026 #include "libmesh/dense_matrix.h" 00027 #include "libmesh/utility.h" 00028 00029 00030 namespace libMesh 00031 { 00032 00033 FESubdivision::FESubdivision(const FEType& fet) : 00034 FE<2,SUBDIVISION>(fet) 00035 { 00036 // Only 2D meshes in 3D space are supported 00037 libmesh_assert_equal_to(LIBMESH_DIM, 3); 00038 } 00039 00040 00041 00042 void FESubdivision::init_subdivision_matrix(DenseMatrix<Real> &A, 00043 unsigned int valence) 00044 { 00045 A.resize(valence + 12, valence + 12); 00046 00047 // A = (S11 0; S21 S22), see Cirak et al., 00048 // Int. J. Numer. Meth. Engng. 2000; 47:2039-2072, Appendix A.2. 00049 00050 // First, set the static S21 part 00051 A(valence+ 1,0 ) = 0.125; 00052 A(valence+ 1,1 ) = 0.375; 00053 A(valence+ 1,valence ) = 0.375; 00054 A(valence+ 2,0 ) = 0.0625; 00055 A(valence+ 2,1 ) = 0.625; 00056 A(valence+ 2,2 ) = 0.0625; 00057 A(valence+ 2,valence ) = 0.0625; 00058 A(valence+ 3,0 ) = 0.125; 00059 A(valence+ 3,1 ) = 0.375; 00060 A(valence+ 3,2 ) = 0.375; 00061 A(valence+ 4,0 ) = 0.0625; 00062 A(valence+ 4,1 ) = 0.0625; 00063 A(valence+ 4,valence-1) = 0.0625; 00064 A(valence+ 4,valence ) = 0.625; 00065 A(valence+ 5,0 ) = 0.125; 00066 A(valence+ 5,valence-1) = 0.375; 00067 A(valence+ 5,valence ) = 0.375; 00068 A(valence+ 6,1 ) = 0.375; 00069 A(valence+ 6,valence ) = 0.125; 00070 A(valence+ 7,1 ) = 0.375; 00071 A(valence+ 8,1 ) = 0.375; 00072 A(valence+ 8,2 ) = 0.125; 00073 A(valence+ 9,1 ) = 0.125; 00074 A(valence+ 9,valence ) = 0.375; 00075 A(valence+10,valence ) = 0.375; 00076 A(valence+11,valence-1) = 0.125; 00077 A(valence+11,valence ) = 0.375; 00078 00079 // Next, set the static S22 part 00080 A(valence+ 1,valence+1) = 0.125; 00081 A(valence+ 2,valence+1) = 0.0625; 00082 A(valence+ 2,valence+2) = 0.0625; 00083 A(valence+ 2,valence+3) = 0.0625; 00084 A(valence+ 3,valence+3) = 0.125; 00085 A(valence+ 4,valence+1) = 0.0625; 00086 A(valence+ 4,valence+4) = 0.0625; 00087 A(valence+ 4,valence+5) = 0.0625; 00088 A(valence+ 5,valence+5) = 0.125; 00089 A(valence+ 6,valence+1) = 0.375; 00090 A(valence+ 6,valence+2) = 0.125; 00091 A(valence+ 7,valence+1) = 0.125; 00092 A(valence+ 7,valence+2) = 0.375; 00093 A(valence+ 7,valence+3) = 0.125; 00094 A(valence+ 8,valence+2) = 0.125; 00095 A(valence+ 8,valence+3) = 0.375; 00096 A(valence+ 9,valence+1) = 0.375; 00097 A(valence+ 9,valence+4) = 0.125; 00098 A(valence+10,valence+1) = 0.125; 00099 A(valence+10,valence+4) = 0.375; 00100 A(valence+10,valence+5) = 0.125; 00101 A(valence+11,valence+4) = 0.125; 00102 A(valence+11,valence+5) = 0.375; 00103 00104 // Last, set the S11 part: first row 00105 std::vector<Real> weights; 00106 loop_subdivision_mask(weights, valence); 00107 for (unsigned int i = 0; i <= valence; ++i) 00108 A(0,i) = weights[i]; 00109 00110 // second row 00111 A(1,0) = 0.375; 00112 A(1,1) = 0.375; 00113 A(1,2) = 0.125; 00114 A(1,valence) = 0.125; 00115 00116 // third to second-to-last rows 00117 for (unsigned int i = 2; i < valence; ++i) 00118 { 00119 A(i,0 ) = 0.375; 00120 A(i,i-1) = 0.125; 00121 A(i,i ) = 0.375; 00122 A(i,i+1) = 0.125; 00123 } 00124 00125 // last row 00126 A(valence,0) = 0.375; 00127 A(valence,1) = 0.125; 00128 A(valence,valence-1) = 0.125; 00129 A(valence,valence ) = 0.375; 00130 } 00131 00132 00133 00134 Real FESubdivision::regular_shape(const unsigned int i, 00135 const Real v, 00136 const Real w) 00137 { 00138 // These are the 12 quartic box splines, see Cirak et al., 00139 // Int. J. Numer. Meth. Engng. 2000; 47:2039-2072, Appendix A.1. 00140 00141 const Real u = 1 - v - w; 00142 libmesh_assert_less_equal(0, v); 00143 libmesh_assert_less_equal(0, w); 00144 libmesh_assert_less_equal(0, u); 00145 00146 using Utility::pow; 00147 const Real factor = 1. / 12; 00148 00149 switch (i) 00150 { 00151 case 0: 00152 return factor*(pow<4>(u) + 2*u*u*u*v); 00153 case 1: 00154 return factor*(pow<4>(u) + 2*u*u*u*w); 00155 case 2: 00156 return factor*(pow<4>(u) + 2*u*u*u*w + 6*u*u*u*v + 6*u*u*v*w + 12*u*u*v*v + 6*u*v*v*w + 6*u*v*v*v + 00157 2*v*v*v*w + pow<4>(v)); 00158 case 3: 00159 return factor*(6*pow<4>(u) + 24*u*u*u*w + 24*u*u*w*w + 8*u*w*w*w + pow<4>(w) + 24*u*u*u*v + 00160 60*u*u*v*w + 36*u*v*w*w + 6*v*w*w*w + 24*u*u*v*v + 36*u*v*v*w + 12*v*v*w*w + 8*u*v*v*v + 00161 6*v*v*v*w + pow<4>(v)); 00162 case 4: 00163 return factor*(pow<4>(u) + 6*u*u*u*w + 12*u*u*w*w + 6*u*w*w*w + pow<4>(w) + 2*u*u*u*v + 6*u*u*v*w + 00164 6*u*v*w*w + 2*v*w*w*w); 00165 case 5: 00166 return factor*(2*u*v*v*v + pow<4>(v)); 00167 case 6: 00168 return factor*(pow<4>(u) + 6*u*u*u*w + 12*u*u*w*w + 6*u*w*w*w + pow<4>(w) + 8*u*u*u*v + 36*u*u*v*w + 00169 36*u*v*w*w + 8*v*w*w*w + 24*u*u*v*v + 60*u*v*v*w + 24*v*v*w*w + 24*u*v*v*v + 24*v*v*v*w + 6*pow<4>(v)); 00170 case 7: 00171 return factor*(pow<4>(u) + 8*u*u*u*w + 24*u*u*w*w + 24*u*w*w*w + 6*pow<4>(w) + 6*u*u*u*v + 36*u*u*v*w + 00172 60*u*v*w*w + 24*v*w*w*w + 12*u*u*v*v + 36*u*v*v*w + 24*v*v*w*w + 6*u*v*v*v + 8*v*v*v*w + pow<4>(v)); 00173 case 8: 00174 return factor*(2*u*w*w*w + pow<4>(w)); 00175 case 9: 00176 return factor*(2*v*v*v*w + pow<4>(v)); 00177 case 10: 00178 return factor*(2*u*w*w*w + pow<4>(w) + 6*u*v*w*w + 6*v*w*w*w + 6*u*v*v*w + 12*v*v*w*w + 2*u*v*v*v + 00179 6*v*v*v*w + pow<4>(v)); 00180 case 11: 00181 return factor*(pow<4>(w) + 2*v*w*w*w); 00182 00183 default: 00184 libmesh_error_msg("Invalid i = " << i); 00185 } 00186 00187 libmesh_error_msg("We'll never get here!"); 00188 return 0.; 00189 } 00190 00191 00192 00193 Real FESubdivision::regular_shape_deriv(const unsigned int i, 00194 const unsigned int j, 00195 const Real v, 00196 const Real w) 00197 { 00198 const Real u = 1 - v - w; 00199 const Real factor = 1. / 12; 00200 00201 switch (j) // j=0: xi-directional derivative, j=1: eta-directional derivative 00202 { 00203 case 0: // xi derivatives 00204 { 00205 switch (i) // shape function number 00206 { 00207 case 0: 00208 return factor*(-6*v*u*u - 2*u*u*u); 00209 case 1: 00210 return factor*(-4*u*u*u - 6*u*u*w); 00211 case 2: 00212 return factor*(-2*v*v*v - 6*v*v*u + 6*v*u*u + 2*u*u*u); 00213 case 3: 00214 return factor*(-4*v*v*v - 24*v*v*u - 24*v*u*u - 18*v*v*w - 48*v*u*w - 12*u*u*w - 00215 12*v*w*w - 12*u*w*w - 2*w*w*w); 00216 case 4: 00217 return factor*(-6*v*u*u - 2*u*u*u - 12*v*u*w-12*u*u*w - 6*v*w*w - 18*u*w*w - 4*w*w*w); 00218 case 5: 00219 return factor*(2*v*v*v + 6*v*v*u); 00220 case 6: 00221 return factor*(24*v*v*u + 24*v*u*u + 4*u*u*u + 12*v*v*w + 48*v*u*w + 18*u*u*w + 00222 12*v*w*w + 12*u*w*w + 2*w*w*w); 00223 case 7: 00224 return factor*(-2*v*v*v - 6*v*v*u + 6*v*u*u + 2*u*u*u - 12*v*v*w + 12*u*u*w - 00225 12*v*w*w + 12*u*w*w); 00226 case 8: 00227 return -w*w*w/6; 00228 case 9: 00229 return factor*(4*v*v*v + 6*v*v*w); 00230 case 10: 00231 return factor*(2*v*v*v + 6*v*v*u + 12*v*v*w + 12*v*u*w + 18*v*w*w + 6*u*w*w + 4*w*w*w); 00232 case 11: 00233 return w*w*w/6; 00234 default: 00235 libmesh_error_msg("Invalid i = " << i); 00236 } 00237 } 00238 case 1: // eta derivatives 00239 { 00240 switch (i) // shape function number 00241 { 00242 case 0: 00243 return factor*(-6*v*u*u - 4*u*u*u); 00244 case 1: 00245 return factor*(-2*u*u*u - 6*u*u*w); 00246 case 2: 00247 return factor*(-4*v*v*v - 18*v*v*u - 12*v*u*u - 2*u*u*u - 6*v*v*w - 12*v*u*w - 00248 6*u*u*w); 00249 case 3: 00250 return factor*(-2*v*v*v-12*v*v*u - 12*v*u*u - 12*v*v*w - 48*v*u*w - 24*u*u*w - 00251 18*v*w*w - 24*u*w*w - 4*w*w*w); 00252 case 4: 00253 return factor*(2*u*u*u + 6*u*u*w - 6*u*w*w - 2*w*w*w); 00254 case 5: 00255 return -v*v*v/6; 00256 case 6: 00257 return factor*(12*v*v*u + 12*v*u*u + 2*u*u*u - 12*v*v*w + 6*u*u*w - 12*v*w*w - 00258 6*u*w*w - 2*w*w*w); 00259 case 7: 00260 return factor*(2*v*v*v + 12*v*v*u + 18*v*u*u + 4*u*u*u + 12*v*v*w + 48*v*u*w + 00261 24*u*u*w + 12*v*w*w + 24*u*w*w); 00262 case 8: 00263 return factor*(6*u*w*w + 2*w*w*w); 00264 case 9: 00265 return v*v*v/6; 00266 case 10: 00267 return factor*(4*v*v*v + 6*v*v*u + 18*v*v*w + 12*v*u*w + 12*v*w*w + 6*u*w*w + 00268 2*w*w*w); 00269 case 11: 00270 return factor*(6*v*w*w + 4*w*w*w); 00271 default: 00272 libmesh_error_msg("Invalid i = " << i); 00273 } 00274 } 00275 default: 00276 libmesh_error_msg("Invalid j = " << j); 00277 } 00278 00279 libmesh_error_msg("We'll never get here!"); 00280 return 0.; 00281 } 00282 00283 00284 00285 Real FESubdivision::regular_shape_second_deriv(const unsigned int i, 00286 const unsigned int j, 00287 const Real v, 00288 const Real w) 00289 { 00290 const Real u = 1 - v - w; 00291 const Real factor = 1. / 12; 00292 00293 switch (j) 00294 { 00295 case 0: // xi-xi derivative 00296 { 00297 switch (i) // shape function number 00298 { 00299 case 0: 00300 return v*u; 00301 case 1: 00302 return u*u + u*w; 00303 case 2: 00304 return -2*v*u; 00305 case 3: 00306 return v*v - 2*u*u + v*w - 2*u*w; 00307 case 4: 00308 return v*u + v*w + u*w + w*w; 00309 case 5: 00310 return v*u; 00311 case 6: 00312 return factor*(-24*v*v + 12*u*u - 24*v*w + 12*u*w); 00313 case 7: 00314 return -2*v*u - 2*v*w - 2*u*w - 2*w*w; 00315 case 8: 00316 return 0.; 00317 case 9: 00318 return v*v + v*w; 00319 case 10: 00320 return v*u + v*w + u*w + w*w; 00321 case 11: 00322 return 0.; 00323 default: 00324 libmesh_error_msg("Invalid i = " << i); 00325 } 00326 } 00327 case 1: //eta-xi derivative 00328 { 00329 switch (i) 00330 { 00331 case 0: 00332 return factor*(12*v*u + 6*u*u); 00333 case 1: 00334 return factor*(6*u*u + 12*u*w); 00335 case 2: 00336 return factor*(6*v*v - 12*v*u - 6*u*u); 00337 case 3: 00338 return factor*(6*v*v - 12*u*u + 24*v*w + 6*w*w); 00339 case 4: 00340 return factor*(-6*u*u - 12*u*w + 6*w*w); 00341 case 5: 00342 return -v*v/2.; 00343 case 6: 00344 return factor*(-12*v*v + 6*u*u - 24*v*w - 12*u*w - 6*w*w); 00345 case 7: 00346 return factor*(-6*v*v - 12*v*u + 6*u*u - 24*v*w - 12*w*w); 00347 case 8: 00348 return -w*w/2.; 00349 case 9: 00350 return v*v/2.; 00351 case 10: 00352 return factor*(6*v*v + 12*v*u + 24*v*w + 12*u*w + 6*w*w); 00353 case 11: 00354 return w*w/2.; 00355 default: 00356 libmesh_error_msg("Invalid i = " << i); 00357 } 00358 } 00359 case 2: // eta-eta derivative 00360 { 00361 switch (i) 00362 { 00363 case 0: 00364 return v*u + u*u; 00365 case 1: 00366 return u*w; 00367 case 2: 00368 return v*v + v*u + v*w + u*w; 00369 case 3: 00370 return -2*v*u - 2*u*u + v*w + w*w; 00371 case 4: 00372 return -2*u*w; 00373 case 5: 00374 return 0.; 00375 case 6: 00376 return -2*v*v - 2*v*u - 2*v*w - 2*u*w; 00377 case 7: 00378 return v*u + u*u - 2*v*w - 2*w*w; 00379 case 8: 00380 return u*w; 00381 case 9: 00382 return 0.; 00383 case 10: 00384 return v*v + v*u + v*w + u*w; 00385 case 11: 00386 return v*w + w*w; 00387 default: 00388 libmesh_error_msg("Invalid i = " << i); 00389 } 00390 } 00391 default: 00392 libmesh_error_msg("Invalid j = " << j); 00393 } 00394 00395 libmesh_error_msg("We'll never get here!"); 00396 return 0.; 00397 } 00398 00399 00400 00401 void FESubdivision::loop_subdivision_mask(std::vector<Real> & weights, 00402 const unsigned int valence) 00403 { 00404 libmesh_assert_greater(valence, 0); 00405 const Real cs = std::cos(2 * libMesh::pi / valence); 00406 const Real nb_weight = (0.625 - Utility::pow<2>(0.375 + 0.25 * cs)) / valence; 00407 weights.resize(1 + valence, nb_weight); 00408 weights[0] = 1.0 - valence * nb_weight; 00409 } 00410 00411 00412 00413 void FESubdivision::init_shape_functions(const std::vector<Point> &qp, 00414 const Elem *elem) 00415 { 00416 libmesh_assert(elem); 00417 libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION); 00418 const Tri3Subdivision* sd_elem = static_cast<const Tri3Subdivision*>(elem); 00419 00420 START_LOG("init_shape_functions()", "FESubdivision"); 00421 00422 calculations_started = true; 00423 00424 // If the user forgot to request anything, we'll be safe and calculate everything: 00425 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00426 if (!calculate_phi && !calculate_dphi && !calculate_d2phi) 00427 calculate_phi = calculate_dphi = calculate_d2phi = true; 00428 #else 00429 if (!calculate_phi && !calculate_dphi) 00430 calculate_phi = calculate_dphi = true; 00431 #endif 00432 00433 const unsigned int valence = sd_elem->get_ordered_valence(0); 00434 const unsigned int n_qp = cast_int<unsigned int>(qp.size()); 00435 const unsigned int n_approx_shape_functions = valence + 6; 00436 00437 // resize the vectors to hold current data 00438 phi.resize (n_approx_shape_functions); 00439 dphi.resize (n_approx_shape_functions); 00440 dphidxi.resize (n_approx_shape_functions); 00441 dphideta.resize (n_approx_shape_functions); 00442 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00443 d2phi.resize (n_approx_shape_functions); 00444 d2phidxi2.resize (n_approx_shape_functions); 00445 d2phidxideta.resize(n_approx_shape_functions); 00446 d2phideta2.resize (n_approx_shape_functions); 00447 #endif 00448 00449 for (unsigned int i = 0; i < n_approx_shape_functions; ++i) 00450 { 00451 phi[i].resize (n_qp); 00452 dphi[i].resize (n_qp); 00453 dphidxi[i].resize (n_qp); 00454 dphideta[i].resize (n_qp); 00455 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00456 d2phi[i].resize (n_qp); 00457 d2phidxi2[i].resize (n_qp); 00458 d2phidxideta[i].resize(n_qp); 00459 d2phideta2[i].resize (n_qp); 00460 #endif 00461 } 00462 00463 // Renumbering of the shape functions 00464 static const unsigned int cvi[12] = {3,6,2,0,1,4,7,10,9,5,11,8}; 00465 00466 if (valence == 6) // This means that all vertices are regular, i.e. we have 12 shape functions 00467 { 00468 for (unsigned int i = 0; i < n_approx_shape_functions; ++i) 00469 { 00470 for (unsigned int p = 0; p < n_qp; ++p) 00471 { 00472 phi[i][p] = FE<2,SUBDIVISION>::shape (elem, fe_type.order, cvi[i], qp[p]); 00473 dphidxi[i][p] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, cvi[i], 0, qp[p]); 00474 dphideta[i][p] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, cvi[i], 1, qp[p]); 00475 dphi[i][p](0) = dphidxi[i][p]; 00476 dphi[i][p](1) = dphideta[i][p]; 00477 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00478 d2phidxi2[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 0, qp[p]); 00479 d2phidxideta[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 1, qp[p]); 00480 d2phideta2[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 2, qp[p]); 00481 d2phi[i][p](0,0) = d2phidxi2[i][p]; 00482 d2phi[i][p](0,1) = d2phi[i][p](1,0) = d2phidxideta[i][p]; 00483 d2phi[i][p](1,1) = d2phideta2[i][p]; 00484 #endif 00485 } 00486 } 00487 } 00488 else // vertex 0 is irregular by construction of the mesh 00489 { 00490 static const Real eps = 1e-10; 00491 00492 // temporary values 00493 std::vector<Real> tphi(12); 00494 std::vector<Real> tdphidxi(12); 00495 std::vector<Real> tdphideta(12); 00496 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00497 std::vector<Real> td2phidxi2(12); 00498 std::vector<Real> td2phidxideta(12); 00499 std::vector<Real> td2phideta2(12); 00500 #endif 00501 00502 for (unsigned int p = 0; p < n_qp; ++p) 00503 { 00504 // evaluate the number of the required subdivisions 00505 Real v = qp[p](0); 00506 Real w = qp[p](1); 00507 Real u = 1 - v - w; 00508 Real min = 0, max = 0.5; 00509 int n = 0; 00510 while (!(u > min-eps && u < max+eps)) 00511 { 00512 ++n; 00513 min = max; 00514 max += std::pow((Real)(2), -n-1); 00515 } 00516 00517 // transform u, v and w according to the number of subdivisions required. 00518 const Real pow2 = std::pow((Real)(2), n); 00519 v *= pow2; 00520 w *= pow2; 00521 u = 1 - v - w; 00522 libmesh_assert_less(u, 0.5 + eps); 00523 libmesh_assert_greater(u, -eps); 00524 00525 // find out in which subdivided patch we are and setup the "selection matrix" P and the transformation Jacobian 00526 // (see Int. J. Numer. Meth. Engng. 2000; 47:2039-2072, Appendix A.2.) 00527 const int k = n+1; 00528 Real jfac; // the additional factor per derivative order 00529 DenseMatrix<Real> P(12, valence+12); 00530 if (v > 0.5 - eps) 00531 { 00532 v = 2*v - 1; 00533 w = 2*w; 00534 jfac = std::pow((Real)(2), k); 00535 P( 0,2 ) = 1; 00536 P( 1,0 ) = 1; 00537 P( 2,valence+3) = 1; 00538 P( 3,1 ) = 1; 00539 P( 4,valence ) = 1; 00540 P( 5,valence+8) = 1; 00541 P( 6,valence+2) = 1; 00542 P( 7,valence+1) = 1; 00543 P( 8,valence+4) = 1; 00544 P( 9,valence+7) = 1; 00545 P(10,valence+6) = 1; 00546 P(11,valence+9) = 1; 00547 } 00548 else if (w > 0.5 - eps) 00549 { 00550 v = 2*v; 00551 w = 2*w - 1; 00552 jfac = std::pow((Real)(2), k); 00553 P( 0,0 ) = 1; 00554 P( 1,valence- 1) = 1; 00555 P( 2,1 ) = 1; 00556 P( 3,valence ) = 1; 00557 P( 4,valence+ 5) = 1; 00558 P( 5,valence+ 2) = 1; 00559 P( 6,valence+ 1) = 1; 00560 P( 7,valence+ 4) = 1; 00561 P( 8,valence+11) = 1; 00562 P( 9,valence+ 6) = 1; 00563 P(10,valence+ 9) = 1; 00564 P(11,valence+10) = 1; 00565 } 00566 else 00567 { 00568 v = 1 - 2*v; 00569 w = 1 - 2*w; 00570 jfac = std::pow((Real)(-2), k); 00571 P( 0,valence+9) = 1; 00572 P( 1,valence+6) = 1; 00573 P( 2,valence+4) = 1; 00574 P( 3,valence+1) = 1; 00575 P( 4,valence+2) = 1; 00576 P( 5,valence+5) = 1; 00577 P( 6,valence ) = 1; 00578 P( 7,1 ) = 1; 00579 P( 8,valence+3) = 1; 00580 P( 9,valence-1) = 1; 00581 P(10,0 ) = 1; 00582 P(11,2 ) = 1; 00583 } 00584 00585 u = 1 - v - w; 00586 if ((u > 1 + eps) || (u < -eps)) 00587 libmesh_error_msg("SUBDIVISION irregular patch: u is outside valid range!"); 00588 00589 DenseMatrix<Real> A; 00590 init_subdivision_matrix(A, valence); 00591 00592 // compute P*A^k 00593 if (k > 1) 00594 { 00595 DenseMatrix<Real> Acopy(A); 00596 for (int e = 1; e < k; ++e) 00597 A.right_multiply(Acopy); 00598 } 00599 P.right_multiply(A); 00600 00601 const Point transformed_p(v,w); 00602 00603 for (unsigned int i = 0; i < 12; ++i) 00604 { 00605 tphi[i] = FE<2,SUBDIVISION>::shape (elem, fe_type.order, i, transformed_p); 00606 tdphidxi[i] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, i, 0, transformed_p); 00607 tdphideta[i] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, i, 1, transformed_p); 00608 00609 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00610 td2phidxi2[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 0, transformed_p); 00611 td2phidxideta[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 1, transformed_p); 00612 td2phideta2[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 2, transformed_p); 00613 #endif 00614 } 00615 00616 // Finally, we can compute the irregular shape functions as the product of P 00617 // and the regular shape functions: 00618 Real sum1, sum2, sum3, sum4, sum5, sum6; 00619 for (unsigned int j = 0; j < n_approx_shape_functions; ++j) 00620 { 00621 sum1 = sum2 = sum3 = sum4 = sum5 = sum6 = 0; 00622 for (unsigned int i = 0; i < 12; ++i) 00623 { 00624 sum1 += P(i,j) * tphi[i]; 00625 sum2 += P(i,j) * tdphidxi[i]; 00626 sum3 += P(i,j) * tdphideta[i]; 00627 00628 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00629 sum4 += P(i,j) * td2phidxi2[i]; 00630 sum5 += P(i,j) * td2phidxideta[i]; 00631 sum6 += P(i,j) * td2phideta2[i]; 00632 #endif 00633 } 00634 phi[j][p] = sum1; 00635 dphidxi[j][p] = sum2 * jfac; 00636 dphideta[j][p] = sum3 * jfac; 00637 dphi[j][p](0) = dphidxi[j][p]; 00638 dphi[j][p](1) = dphideta[j][p]; 00639 00640 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00641 d2phidxi2[j][p] = sum4 * jfac * jfac; 00642 d2phidxideta[j][p] = sum5 * jfac * jfac; 00643 d2phideta2[j][p] = sum6 * jfac * jfac; 00644 d2phi[j][p](0,0) = d2phidxi2[j][p]; 00645 d2phi[j][p](0,1) = d2phi[j][p](1,0) = d2phidxideta[j][p]; 00646 d2phi[j][p](1,1) = d2phideta2[j][p]; 00647 #endif 00648 } 00649 } // end quadrature loop 00650 } // end irregular vertex 00651 00652 // Let the FEMap use the same initialized shape functions 00653 this->_fe_map->get_phi_map() = phi; 00654 this->_fe_map->get_dphidxi_map() = dphidxi; 00655 this->_fe_map->get_dphideta_map() = dphideta; 00656 00657 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES 00658 this->_fe_map->get_d2phidxi2_map() = d2phidxi2; 00659 this->_fe_map->get_d2phideta2_map() = d2phideta2; 00660 this->_fe_map->get_d2phidxideta_map() = d2phidxideta; 00661 #endif 00662 00663 STOP_LOG("init_shape_functions()", "FESubdivision"); 00664 } 00665 00666 00667 00668 void FESubdivision::attach_quadrature_rule(QBase *q) 00669 { 00670 libmesh_assert(q); 00671 00672 qrule = q; 00673 // make sure we don't cache results from a previous quadrature rule 00674 elem_type = INVALID_ELEM; 00675 return; 00676 } 00677 00678 00679 00680 void FESubdivision::reinit(const Elem* elem, 00681 const std::vector<Point>* const libmesh_dbg_var(pts), 00682 const std::vector<Real>* const) 00683 { 00684 libmesh_assert(elem); 00685 libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION); 00686 #ifndef NDEBUG 00687 const Tri3Subdivision* sd_elem = static_cast<const Tri3Subdivision*>(elem); 00688 #endif 00689 00690 START_LOG("reinit()", "FESubdivision"); 00691 00692 libmesh_assert(!sd_elem->is_ghost()); 00693 libmesh_assert(sd_elem->is_subdivision_updated()); 00694 00695 // check if vertices 1 and 2 are regular 00696 libmesh_assert_equal_to(sd_elem->get_ordered_valence(1), 6); 00697 libmesh_assert_equal_to(sd_elem->get_ordered_valence(2), 6); 00698 00699 // no custom quadrature support 00700 libmesh_assert(pts == NULL); 00701 libmesh_assert(qrule); 00702 qrule->init(elem->type()); 00703 00704 // Initialize the shape functions 00705 this->init_shape_functions(this->qrule->get_points(), elem); 00706 00707 // The shape functions correspond to the qrule 00708 shapes_on_quadrature = true; 00709 00710 // Compute the map for this element. 00711 this->_fe_map->compute_map (this->dim, this->qrule->get_weights(), elem); 00712 00713 STOP_LOG("reinit()", "FESubdivision"); 00714 } 00715 00716 00717 00718 template <> 00719 Real FE<2,SUBDIVISION>::shape(const ElemType type, 00720 const Order order, 00721 const unsigned int i, 00722 const Point& p) 00723 { 00724 switch (order) 00725 { 00726 case FOURTH: 00727 { 00728 switch (type) 00729 { 00730 case TRI3SUBDIVISION: 00731 libmesh_assert_less(i, 12); 00732 return FESubdivision::regular_shape(i,p(0),p(1)); 00733 default: 00734 libmesh_error_msg("ERROR: Unsupported element type!"); 00735 } 00736 } 00737 default: 00738 libmesh_error_msg("ERROR: Unsupported polynomial order!"); 00739 } 00740 00741 libmesh_error_msg("We'll never get here!"); 00742 return 0.; 00743 } 00744 00745 00746 00747 template <> 00748 Real FE<2,SUBDIVISION>::shape(const Elem* elem, 00749 const Order order, 00750 const unsigned int i, 00751 const Point& p) 00752 { 00753 libmesh_assert(elem); 00754 return FE<2,SUBDIVISION>::shape(elem->type(), order, i, p); 00755 } 00756 00757 00758 00759 template <> 00760 Real FE<2,SUBDIVISION>::shape_deriv(const ElemType type, 00761 const Order order, 00762 const unsigned int i, 00763 const unsigned int j, 00764 const Point& p) 00765 { 00766 switch (order) 00767 { 00768 case FOURTH: 00769 { 00770 switch (type) 00771 { 00772 case TRI3SUBDIVISION: 00773 libmesh_assert_less(i, 12); 00774 return FESubdivision::regular_shape_deriv(i,j,p(0),p(1)); 00775 default: 00776 libmesh_error_msg("ERROR: Unsupported element type!"); 00777 } 00778 } 00779 default: 00780 libmesh_error_msg("ERROR: Unsupported polynomial order!"); 00781 } 00782 00783 libmesh_error_msg("We'll never get here!"); 00784 return 0.; 00785 } 00786 00787 00788 00789 template <> 00790 Real FE<2,SUBDIVISION>::shape_deriv(const Elem* elem, 00791 const Order order, 00792 const unsigned int i, 00793 const unsigned int j, 00794 const Point& p) 00795 { 00796 libmesh_assert(elem); 00797 return FE<2,SUBDIVISION>::shape_deriv(elem->type(), order, i, j, p); 00798 } 00799 00800 00801 00802 template <> 00803 Real FE<2,SUBDIVISION>::shape_second_deriv(const ElemType type, 00804 const Order order, 00805 const unsigned int i, 00806 const unsigned int j, 00807 const Point& p) 00808 { 00809 switch (order) 00810 { 00811 case FOURTH: 00812 { 00813 switch (type) 00814 { 00815 case TRI3SUBDIVISION: 00816 libmesh_assert_less(i, 12); 00817 return FESubdivision::regular_shape_second_deriv(i,j,p(0),p(1)); 00818 default: 00819 libmesh_error_msg("ERROR: Unsupported element type!"); 00820 } 00821 } 00822 default: 00823 libmesh_error_msg("ERROR: Unsupported polynomial order!"); 00824 } 00825 00826 libmesh_error_msg("We'll never get here!"); 00827 return 0.; 00828 } 00829 00830 00831 00832 template <> 00833 Real FE<2,SUBDIVISION>::shape_second_deriv(const Elem* elem, 00834 const Order order, 00835 const unsigned int i, 00836 const unsigned int j, 00837 const Point& p) 00838 { 00839 libmesh_assert(elem); 00840 return FE<2,SUBDIVISION>::shape_second_deriv(elem->type(), order, i, j, p); 00841 } 00842 00843 00844 00845 template <> 00846 void FE<2,SUBDIVISION>::nodal_soln(const Elem* elem, 00847 const Order, 00848 const std::vector<Number>& elem_soln, 00849 std::vector<Number>& nodal_soln) 00850 { 00851 libmesh_assert(elem); 00852 libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION); 00853 const Tri3Subdivision* sd_elem = static_cast<const Tri3Subdivision*>(elem); 00854 00855 nodal_soln.resize(3); // three nodes per element 00856 00857 // Ghost nodes are auxiliary. 00858 if (sd_elem->is_ghost()) 00859 { 00860 nodal_soln[0] = 0; 00861 nodal_soln[1] = 0; 00862 nodal_soln[2] = 0; 00863 return; 00864 } 00865 00866 // First node (node 0 in the element patch): 00867 unsigned int j = sd_elem->local_node_number(sd_elem->get_ordered_node(0)->id()); 00868 nodal_soln[j] = elem_soln[0]; 00869 00870 // Second node (node 1 in the element patch): 00871 j = sd_elem->local_node_number(sd_elem->get_ordered_node(1)->id()); 00872 nodal_soln[j] = elem_soln[1]; 00873 00874 // Third node (node 'valence' in the element patch): 00875 j = sd_elem->local_node_number(sd_elem->get_ordered_node(2)->id()); 00876 nodal_soln[j] = elem_soln[sd_elem->get_ordered_valence(0)]; 00877 } 00878 00879 00880 00881 // the empty template specializations below are needed to avoid 00882 // linker reference errors, but should never get called 00883 template <> 00884 void FE<2,SUBDIVISION>::side_map(const Elem*, 00885 const Elem*, 00886 const unsigned int, 00887 const std::vector<Point>&, 00888 std::vector<Point>&) 00889 { 00890 libmesh_not_implemented(); 00891 } 00892 00893 template <> 00894 void FE<2,SUBDIVISION>::edge_reinit(Elem const*, 00895 unsigned int, 00896 Real, 00897 const std::vector<Point>* const, 00898 const std::vector<Real>* const) 00899 { 00900 libmesh_not_implemented(); 00901 } 00902 00903 template <> 00904 Point FE<2,SUBDIVISION>::inverse_map(const Elem*, 00905 const Point&, 00906 const Real, 00907 const bool) 00908 { 00909 libmesh_not_implemented(); 00910 } 00911 00912 template <> 00913 void FE<2,SUBDIVISION>::inverse_map(const Elem*, 00914 const std::vector<Point>&, 00915 std::vector<Point>&, 00916 Real, 00917 bool) 00918 { 00919 libmesh_not_implemented(); 00920 } 00921 00922 00923 00924 // The number of dofs per subdivision element depends on the valence of its nodes, so it is non-static 00925 template <> unsigned int FE<2,SUBDIVISION>::n_dofs(const ElemType, const Order) { libmesh_not_implemented(); return 0; } 00926 00927 // Loop subdivision elements have only a single dof per node 00928 template <> unsigned int FE<2,SUBDIVISION>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 1; } 00929 00930 // Subdivision FEMs have dofs only at the nodes 00931 template <> unsigned int FE<2,SUBDIVISION>::n_dofs_per_elem(const ElemType, const Order) { return 0; } 00932 00933 // Subdivision FEMs have dofs only at the nodes 00934 template <> void FE<2,SUBDIVISION>::dofs_on_side(const Elem *const, const Order, unsigned int, std::vector<unsigned int> &di) { di.resize(0); } 00935 template <> void FE<2,SUBDIVISION>::dofs_on_edge(const Elem *const, const Order, unsigned int, std::vector<unsigned int> &di) { di.resize(0); } 00936 00937 // Subdivision FEMs are C^1 continuous 00938 template <> FEContinuity FE<2,SUBDIVISION>::get_continuity() const { return C_ONE; } 00939 00940 // Subdivision FEMs are not hierarchic 00941 template <> bool FE<2,SUBDIVISION>::is_hierarchic() const { return false; } 00942 00943 // Subdivision FEM shapes need reinit 00944 template <> bool FE<2,SUBDIVISION>::shapes_need_reinit() const { return true; } 00945 00946 } // namespace libMesh
generated by
