MFEM v2.0
|
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at 00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights 00003 // reserved. See file COPYRIGHT for details. 00004 // 00005 // This file is part of the MFEM library. For more information and source code 00006 // availability see http://mfem.googlecode.com. 00007 // 00008 // MFEM is free software; you can redistribute it and/or modify it under the 00009 // terms of the GNU Lesser General Public License (as published by the Free 00010 // Software Foundation) version 2.1 dated February 1999. 00011 00012 #include "fem.hpp" 00013 #include <cstdlib> 00014 #include <cstring> 00015 #include <cstdio> 00016 00017 int FiniteElementCollection::HasFaceDofs(int GeomType) const 00018 { 00019 switch (GeomType) 00020 { 00021 case Geometry::TETRAHEDRON: return DofForGeometry (Geometry::TRIANGLE); 00022 case Geometry::CUBE: return DofForGeometry (Geometry::SQUARE); 00023 default: 00024 mfem_error ("FiniteElementCollection::HasFaceDofs:" 00025 " unknown geometry type."); 00026 } 00027 return 0; 00028 } 00029 00030 FiniteElementCollection *FiniteElementCollection::New(const char *name) 00031 { 00032 FiniteElementCollection *fec = NULL; 00033 00034 if (!strcmp(name, "Linear")) 00035 fec = new LinearFECollection; 00036 else if (!strcmp(name, "Quadratic")) 00037 fec = new QuadraticFECollection; 00038 else if (!strcmp(name, "QuadraticPos")) 00039 fec = new QuadraticPosFECollection; 00040 else if (!strcmp(name, "Cubic")) 00041 fec = new CubicFECollection; 00042 else if (!strcmp(name, "Const3D")) 00043 fec = new Const3DFECollection; 00044 else if (!strcmp(name, "Const2D")) 00045 fec = new Const2DFECollection; 00046 else if (!strcmp(name, "LinearDiscont2D")) 00047 fec = new LinearDiscont2DFECollection; 00048 else if (!strcmp(name, "GaussLinearDiscont2D")) 00049 fec = new GaussLinearDiscont2DFECollection; 00050 else if (!strcmp(name, "P1OnQuad")) 00051 fec = new P1OnQuadFECollection; 00052 else if (!strcmp(name, "QuadraticDiscont2D")) 00053 fec = new QuadraticDiscont2DFECollection; 00054 else if (!strcmp(name, "QuadraticPosDiscont2D")) 00055 fec = new QuadraticPosDiscont2DFECollection; 00056 else if (!strcmp(name, "GaussQuadraticDiscont2D")) 00057 fec = new GaussQuadraticDiscont2DFECollection; 00058 else if (!strcmp(name, "CubicDiscont2D")) 00059 fec = new CubicDiscont2DFECollection; 00060 else if (!strcmp(name, "LinearDiscont3D")) 00061 fec = new LinearDiscont3DFECollection; 00062 else if (!strcmp(name, "QuadraticDiscont3D")) 00063 fec = new QuadraticDiscont3DFECollection; 00064 else if (!strcmp(name, "LinearNonConf3D")) 00065 fec = new LinearNonConf3DFECollection; 00066 else if (!strcmp(name, "CrouzeixRaviart")) 00067 fec = new CrouzeixRaviartFECollection; 00068 else if (!strcmp(name, "ND1_3D")) 00069 fec = new ND1_3DFECollection; 00070 else if (!strcmp(name, "RT0_2D")) 00071 fec = new RT0_2DFECollection; 00072 else if (!strcmp(name, "RT1_2D")) 00073 fec = new RT1_2DFECollection; 00074 else if (!strcmp(name, "RT2_2D")) 00075 fec = new RT2_2DFECollection; 00076 else if (!strcmp(name, "RT0_3D")) 00077 fec = new RT0_3DFECollection; 00078 else if (!strcmp(name, "RT1_3D")) 00079 fec = new RT1_3DFECollection; 00080 else if (!strncmp(name, "H1_", 3)) 00081 fec = new H1_FECollection(atoi(name + 7), atoi(name + 3)); 00082 else if (!strncmp(name, "L2_", 3)) 00083 fec = new L2_FECollection(atoi(name + 7), atoi(name + 3)); 00084 else if (!strncmp(name, "RT_", 3)) 00085 fec = new RT_FECollection(atoi(name + 7), atoi(name + 3)); 00086 else if (!strncmp(name, "ND_", 3)) 00087 fec = new ND_FECollection(atoi(name + 7), atoi(name + 3)); 00088 else if (!strncmp(name, "Local_", 6)) 00089 fec = new Local_FECollection(name + 6); 00090 else if (!strncmp(name, "NURBS", 5)) 00091 fec = new NURBSFECollection(atoi(name + 5)); 00092 else 00093 mfem_error ("FiniteElementCollection::New : " 00094 "Unknown FiniteElementCollection!"); 00095 00096 return fec; 00097 } 00098 00099 const FiniteElement * 00100 LinearFECollection::FiniteElementForGeometry(int GeomType) const 00101 { 00102 switch (GeomType) 00103 { 00104 case Geometry::POINT: return &PointFE; 00105 case Geometry::SEGMENT: return &SegmentFE; 00106 case Geometry::TRIANGLE: return &TriangleFE; 00107 case Geometry::SQUARE: return &QuadrilateralFE; 00108 case Geometry::TETRAHEDRON: return &TetrahedronFE; 00109 case Geometry::CUBE: return &ParallelepipedFE; 00110 default: 00111 mfem_error ("LinearFECollection: unknown geometry type."); 00112 } 00113 return &SegmentFE; // Make some compilers happy 00114 } 00115 00116 int LinearFECollection::DofForGeometry(int GeomType) const 00117 { 00118 switch (GeomType) 00119 { 00120 case Geometry::POINT: return 1; 00121 case Geometry::SEGMENT: return 0; 00122 case Geometry::TRIANGLE: return 0; 00123 case Geometry::SQUARE: return 0; 00124 case Geometry::TETRAHEDRON: return 0; 00125 case Geometry::CUBE: return 0; 00126 default: 00127 mfem_error ("LinearFECollection: unknown geometry type."); 00128 } 00129 return 0; // Make some compilers happy 00130 } 00131 00132 int * LinearFECollection::DofOrderForOrientation(int GeomType, int Or) const 00133 { 00134 return NULL; 00135 } 00136 00137 00138 const FiniteElement * 00139 QuadraticFECollection::FiniteElementForGeometry(int GeomType) const 00140 { 00141 switch (GeomType) 00142 { 00143 case Geometry::POINT: return &PointFE; 00144 case Geometry::SEGMENT: return &SegmentFE; 00145 case Geometry::TRIANGLE: return &TriangleFE; 00146 case Geometry::SQUARE: return &QuadrilateralFE; 00147 case Geometry::TETRAHEDRON: return &TetrahedronFE; 00148 case Geometry::CUBE: return &ParallelepipedFE; 00149 default: 00150 mfem_error ("QuadraticFECollection: unknown geometry type."); 00151 } 00152 return &SegmentFE; // Make some compilers happy 00153 } 00154 00155 int QuadraticFECollection::DofForGeometry(int GeomType) const 00156 { 00157 switch (GeomType) 00158 { 00159 case Geometry::POINT: return 1; 00160 case Geometry::SEGMENT: return 1; 00161 case Geometry::TRIANGLE: return 0; 00162 case Geometry::SQUARE: return 1; 00163 case Geometry::TETRAHEDRON: return 0; 00164 case Geometry::CUBE: return 1; 00165 default: 00166 mfem_error ("QuadraticFECollection: unknown geometry type."); 00167 } 00168 return 0; // Make some compilers happy 00169 } 00170 00171 int * QuadraticFECollection::DofOrderForOrientation(int GeomType, int Or) const 00172 { 00173 static int indexes[] = { 0 }; 00174 00175 return indexes; 00176 } 00177 00178 00179 const FiniteElement * 00180 QuadraticPosFECollection::FiniteElementForGeometry(int GeomType) const 00181 { 00182 switch (GeomType) 00183 { 00184 case Geometry::SEGMENT: return &SegmentFE; 00185 case Geometry::SQUARE: return &QuadrilateralFE; 00186 default: 00187 mfem_error ("QuadraticPosFECollection: unknown geometry type."); 00188 } 00189 return NULL; // Make some compilers happy 00190 } 00191 00192 int QuadraticPosFECollection::DofForGeometry(int GeomType) const 00193 { 00194 switch (GeomType) 00195 { 00196 case Geometry::POINT: return 1; 00197 case Geometry::SEGMENT: return 1; 00198 case Geometry::SQUARE: return 1; 00199 default: 00200 mfem_error ("QuadraticPosFECollection: unknown geometry type."); 00201 } 00202 return 0; // Make some compilers happy 00203 } 00204 00205 int * QuadraticPosFECollection::DofOrderForOrientation(int GeomType, int Or) 00206 const 00207 { 00208 static int indexes[] = { 0 }; 00209 00210 return indexes; 00211 } 00212 00213 00214 const FiniteElement * 00215 CubicFECollection::FiniteElementForGeometry(int GeomType) const 00216 { 00217 switch (GeomType) 00218 { 00219 case Geometry::POINT: return &PointFE; 00220 case Geometry::SEGMENT: return &SegmentFE; 00221 case Geometry::TRIANGLE: return &TriangleFE; 00222 case Geometry::SQUARE: return &QuadrilateralFE; 00223 case Geometry::TETRAHEDRON: return &TetrahedronFE; 00224 case Geometry::CUBE: return &ParallelepipedFE; 00225 default: 00226 mfem_error ("CubicFECollection: unknown geometry type."); 00227 } 00228 return &SegmentFE; // Make some compilers happy 00229 } 00230 00231 int CubicFECollection::DofForGeometry(int GeomType) const 00232 { 00233 switch (GeomType) 00234 { 00235 case Geometry::POINT: return 1; 00236 case Geometry::SEGMENT: return 2; 00237 case Geometry::TRIANGLE: return 1; 00238 case Geometry::SQUARE: return 4; 00239 case Geometry::TETRAHEDRON: return 0; 00240 case Geometry::CUBE: return 8; 00241 default: 00242 mfem_error ("CubicFECollection: unknown geometry type."); 00243 } 00244 return 0; // Make some compilers happy 00245 } 00246 00247 int * CubicFECollection::DofOrderForOrientation(int GeomType, int Or) const 00248 { 00249 if (GeomType == Geometry::SEGMENT) 00250 { 00251 static int ind_pos[] = { 0, 1 }; 00252 static int ind_neg[] = { 1, 0 }; 00253 00254 if (Or < 0) 00255 return ind_neg; 00256 return ind_pos; 00257 } 00258 else if (GeomType == Geometry::TRIANGLE) 00259 { 00260 static int indexes[] = { 0 }; 00261 00262 return indexes; 00263 } 00264 else if (GeomType == Geometry::SQUARE) 00265 { 00266 static int sq_ind[8][4] = {{0, 1, 2, 3}, {0, 2, 1, 3}, 00267 {2, 0, 3, 1}, {1, 0, 3, 2}, 00268 {3, 2, 1, 0}, {3, 1, 2, 0}, 00269 {1, 3, 0, 2}, {2, 3, 0, 1}}; 00270 return sq_ind[Or]; 00271 } 00272 00273 return NULL; 00274 } 00275 00276 00277 const FiniteElement * 00278 CrouzeixRaviartFECollection::FiniteElementForGeometry(int GeomType) const 00279 { 00280 switch (GeomType) 00281 { 00282 case Geometry::SEGMENT: return &SegmentFE; 00283 case Geometry::TRIANGLE: return &TriangleFE; 00284 case Geometry::SQUARE: return &QuadrilateralFE; 00285 default: 00286 mfem_error ("CrouzeixRaviartFECollection: unknown geometry type."); 00287 } 00288 return &SegmentFE; // Make some compilers happy 00289 } 00290 00291 int CrouzeixRaviartFECollection::DofForGeometry(int GeomType) const 00292 { 00293 switch (GeomType) 00294 { 00295 case Geometry::POINT: return 0; 00296 case Geometry::SEGMENT: return 1; 00297 case Geometry::TRIANGLE: return 0; 00298 case Geometry::SQUARE: return 0; 00299 default: 00300 mfem_error ("CrouzeixRaviartFECollection: unknown geometry type."); 00301 } 00302 return 0; // Make some compilers happy 00303 } 00304 00305 int * CrouzeixRaviartFECollection::DofOrderForOrientation(int GeomType, int Or) 00306 const 00307 { 00308 static int indexes[] = { 0 }; 00309 00310 return indexes; 00311 } 00312 00313 00314 const FiniteElement * 00315 RT0_2DFECollection::FiniteElementForGeometry(int GeomType) const 00316 { 00317 switch (GeomType) 00318 { 00319 case Geometry::SEGMENT: return &SegmentFE; 00320 case Geometry::TRIANGLE: return &TriangleFE; 00321 case Geometry::SQUARE: return &QuadrilateralFE; 00322 default: 00323 mfem_error ("RT0_2DFECollection: unknown geometry type."); 00324 } 00325 return &SegmentFE; // Make some compilers happy 00326 } 00327 00328 int RT0_2DFECollection::DofForGeometry(int GeomType) const 00329 { 00330 switch (GeomType) 00331 { 00332 case Geometry::POINT: return 0; 00333 case Geometry::SEGMENT: return 1; 00334 case Geometry::TRIANGLE: return 0; 00335 case Geometry::SQUARE: return 0; 00336 default: 00337 mfem_error ("RT0_2DFECollection: unknown geometry type."); 00338 } 00339 return 0; // Make some compilers happy 00340 } 00341 00342 int * RT0_2DFECollection::DofOrderForOrientation(int GeomType, int Or) 00343 const 00344 { 00345 static int ind_pos[] = { 0 }; 00346 static int ind_neg[] = { -1 }; 00347 00348 if (Or > 0) 00349 return ind_pos; 00350 return ind_neg; 00351 } 00352 00353 00354 const FiniteElement * 00355 RT1_2DFECollection::FiniteElementForGeometry(int GeomType) const 00356 { 00357 switch (GeomType) 00358 { 00359 case Geometry::SEGMENT: return &SegmentFE; 00360 case Geometry::TRIANGLE: return &TriangleFE; 00361 case Geometry::SQUARE: return &QuadrilateralFE; 00362 default: 00363 mfem_error ("RT1_2DFECollection: unknown geometry type."); 00364 } 00365 return &SegmentFE; // Make some compilers happy 00366 } 00367 00368 int RT1_2DFECollection::DofForGeometry(int GeomType) const 00369 { 00370 switch (GeomType) 00371 { 00372 case Geometry::POINT: return 0; 00373 case Geometry::SEGMENT: return 2; 00374 case Geometry::TRIANGLE: return 2; 00375 case Geometry::SQUARE: return 4; 00376 default: 00377 mfem_error ("RT1_2DFECollection: unknown geometry type."); 00378 } 00379 return 0; // Make some compilers happy 00380 } 00381 00382 int * RT1_2DFECollection::DofOrderForOrientation(int GeomType, int Or) 00383 const 00384 { 00385 static int ind_pos[] = { 0, 1 }; 00386 static int ind_neg[] = { -2, -1 }; 00387 00388 if (Or > 0) 00389 return ind_pos; 00390 return ind_neg; 00391 } 00392 00393 const FiniteElement * 00394 RT2_2DFECollection::FiniteElementForGeometry(int GeomType) const 00395 { 00396 switch (GeomType) 00397 { 00398 case Geometry::SEGMENT: return &SegmentFE; 00399 case Geometry::TRIANGLE: return &TriangleFE; 00400 case Geometry::SQUARE: return &QuadrilateralFE; 00401 default: 00402 mfem_error ("RT2_2DFECollection: unknown geometry type."); 00403 } 00404 return &SegmentFE; // Make some compilers happy 00405 } 00406 00407 int RT2_2DFECollection::DofForGeometry(int GeomType) const 00408 { 00409 switch (GeomType) 00410 { 00411 case Geometry::POINT: return 0; 00412 case Geometry::SEGMENT: return 3; 00413 case Geometry::TRIANGLE: return 6; 00414 case Geometry::SQUARE: return 12; 00415 default: 00416 mfem_error ("RT2_2DFECollection: unknown geometry type."); 00417 } 00418 return 0; // Make some compilers happy 00419 } 00420 00421 int * RT2_2DFECollection::DofOrderForOrientation(int GeomType, int Or) 00422 const 00423 { 00424 static int ind_pos[] = { 0, 1, 2 }; 00425 static int ind_neg[] = { -3, -2, -1 }; 00426 00427 if (Or > 0) 00428 return ind_pos; 00429 return ind_neg; 00430 } 00431 00432 00433 const FiniteElement * 00434 Const2DFECollection::FiniteElementForGeometry(int GeomType) const 00435 { 00436 switch (GeomType) 00437 { 00438 case Geometry::TRIANGLE: return &TriangleFE; 00439 case Geometry::SQUARE: return &QuadrilateralFE; 00440 default: 00441 mfem_error ("Const2DFECollection: unknown geometry type."); 00442 } 00443 return &TriangleFE; // Make some compilers happy 00444 } 00445 00446 int Const2DFECollection::DofForGeometry(int GeomType) const 00447 { 00448 switch (GeomType) 00449 { 00450 case Geometry::POINT: return 0; 00451 case Geometry::SEGMENT: return 0; 00452 case Geometry::TRIANGLE: return 1; 00453 case Geometry::SQUARE: return 1; 00454 default: 00455 mfem_error ("Const2DFECollection: unknown geometry type."); 00456 } 00457 return 0; // Make some compilers happy 00458 } 00459 00460 int * Const2DFECollection::DofOrderForOrientation(int GeomType, int Or) 00461 const 00462 { 00463 return NULL; 00464 } 00465 00466 00467 const FiniteElement * 00468 LinearDiscont2DFECollection::FiniteElementForGeometry(int GeomType) const 00469 { 00470 switch (GeomType) 00471 { 00472 case Geometry::TRIANGLE: return &TriangleFE; 00473 case Geometry::SQUARE: return &QuadrilateralFE; 00474 default: 00475 mfem_error ("LinearDiscont2DFECollection: unknown geometry type."); 00476 } 00477 return &TriangleFE; // Make some compilers happy 00478 } 00479 00480 int LinearDiscont2DFECollection::DofForGeometry(int GeomType) const 00481 { 00482 switch (GeomType) 00483 { 00484 case Geometry::POINT: return 0; 00485 case Geometry::SEGMENT: return 0; 00486 case Geometry::TRIANGLE: return 3; 00487 case Geometry::SQUARE: return 4; 00488 default: 00489 mfem_error ("LinearDiscont2DFECollection: unknown geometry type."); 00490 } 00491 return 0; // Make some compilers happy 00492 } 00493 00494 int * LinearDiscont2DFECollection::DofOrderForOrientation(int GeomType, int Or) 00495 const 00496 { 00497 return NULL; 00498 } 00499 00500 00501 const FiniteElement * 00502 GaussLinearDiscont2DFECollection::FiniteElementForGeometry(int GeomType) const 00503 { 00504 switch (GeomType) 00505 { 00506 case Geometry::TRIANGLE: return &TriangleFE; 00507 case Geometry::SQUARE: return &QuadrilateralFE; 00508 default: 00509 mfem_error ("GaussLinearDiscont2DFECollection:" 00510 " unknown geometry type."); 00511 } 00512 return &TriangleFE; // Make some compilers happy 00513 } 00514 00515 int GaussLinearDiscont2DFECollection::DofForGeometry(int GeomType) const 00516 { 00517 switch (GeomType) 00518 { 00519 case Geometry::POINT: return 0; 00520 case Geometry::SEGMENT: return 0; 00521 case Geometry::TRIANGLE: return 3; 00522 case Geometry::SQUARE: return 4; 00523 default: 00524 mfem_error ("GaussLinearDiscont2DFECollection:" 00525 " unknown geometry type."); 00526 } 00527 return 0; // Make some compilers happy 00528 } 00529 00530 int * GaussLinearDiscont2DFECollection::DofOrderForOrientation( 00531 int GeomType, int Or) const 00532 { 00533 return NULL; 00534 } 00535 00536 00537 const FiniteElement * 00538 P1OnQuadFECollection::FiniteElementForGeometry(int GeomType) const 00539 { 00540 if (GeomType != Geometry::SQUARE) 00541 { 00542 mfem_error ("P1OnQuadFECollection: unknown geometry type."); 00543 } 00544 return &QuadrilateralFE; 00545 } 00546 00547 int P1OnQuadFECollection::DofForGeometry(int GeomType) const 00548 { 00549 switch (GeomType) 00550 { 00551 case Geometry::POINT: return 0; 00552 case Geometry::SEGMENT: return 0; 00553 case Geometry::SQUARE: return 3; 00554 default: 00555 mfem_error ("P1OnQuadFECollection: unknown geometry type."); 00556 } 00557 return 0; // Make some compilers happy 00558 } 00559 00560 int * P1OnQuadFECollection::DofOrderForOrientation( 00561 int GeomType, int Or) const 00562 { 00563 return NULL; 00564 } 00565 00566 00567 const FiniteElement * 00568 QuadraticDiscont2DFECollection::FiniteElementForGeometry(int GeomType) const 00569 { 00570 switch (GeomType) 00571 { 00572 case Geometry::TRIANGLE: return &TriangleFE; 00573 case Geometry::SQUARE: return &QuadrilateralFE; 00574 default: 00575 mfem_error ("QuadraticDiscont2DFECollection: unknown geometry type."); 00576 } 00577 return &TriangleFE; // Make some compilers happy 00578 } 00579 00580 int QuadraticDiscont2DFECollection::DofForGeometry(int GeomType) const 00581 { 00582 switch (GeomType) 00583 { 00584 case Geometry::POINT: return 0; 00585 case Geometry::SEGMENT: return 0; 00586 case Geometry::TRIANGLE: return 6; 00587 case Geometry::SQUARE: return 9; 00588 default: 00589 mfem_error ("QuadraticDiscont2DFECollection: unknown geometry type."); 00590 } 00591 return 0; // Make some compilers happy 00592 } 00593 00594 int * QuadraticDiscont2DFECollection::DofOrderForOrientation( 00595 int GeomType, int Or) const 00596 { 00597 return NULL; 00598 } 00599 00600 00601 const FiniteElement * 00602 QuadraticPosDiscont2DFECollection::FiniteElementForGeometry(int GeomType) const 00603 { 00604 switch (GeomType) 00605 { 00606 case Geometry::SQUARE: return &QuadrilateralFE; 00607 default: 00608 mfem_error ("QuadraticPosDiscont2DFECollection: unknown geometry type."); 00609 } 00610 return NULL; // Make some compilers happy 00611 } 00612 00613 int QuadraticPosDiscont2DFECollection::DofForGeometry(int GeomType) const 00614 { 00615 switch (GeomType) 00616 { 00617 case Geometry::POINT: return 0; 00618 case Geometry::SEGMENT: return 0; 00619 case Geometry::SQUARE: return 9; 00620 default: 00621 mfem_error ("QuadraticPosDiscont2DFECollection: unknown geometry type."); 00622 } 00623 return 0; // Make some compilers happy 00624 } 00625 00626 00627 const FiniteElement * 00628 GaussQuadraticDiscont2DFECollection::FiniteElementForGeometry(int GeomType) 00629 const 00630 { 00631 switch (GeomType) 00632 { 00633 case Geometry::TRIANGLE: return &TriangleFE; 00634 case Geometry::SQUARE: return &QuadrilateralFE; 00635 default: 00636 mfem_error ("GaussQuadraticDiscont2DFECollection:" 00637 " unknown geometry type."); 00638 } 00639 return &QuadrilateralFE; // Make some compilers happy 00640 } 00641 00642 int GaussQuadraticDiscont2DFECollection::DofForGeometry(int GeomType) const 00643 { 00644 switch (GeomType) 00645 { 00646 case Geometry::POINT: return 0; 00647 case Geometry::SEGMENT: return 0; 00648 case Geometry::TRIANGLE: return 6; 00649 case Geometry::SQUARE: return 9; 00650 default: 00651 mfem_error ("GaussQuadraticDiscont2DFECollection:" 00652 " unknown geometry type."); 00653 } 00654 return 0; // Make some compilers happy 00655 } 00656 00657 int * GaussQuadraticDiscont2DFECollection::DofOrderForOrientation( 00658 int GeomType, int Or) const 00659 { 00660 return NULL; 00661 } 00662 00663 00664 const FiniteElement * 00665 CubicDiscont2DFECollection::FiniteElementForGeometry(int GeomType) const 00666 { 00667 switch (GeomType) 00668 { 00669 case Geometry::TRIANGLE: return &TriangleFE; 00670 case Geometry::SQUARE: return &QuadrilateralFE; 00671 default: 00672 mfem_error ("CubicDiscont2DFECollection: unknown geometry type."); 00673 } 00674 return &TriangleFE; // Make some compilers happy 00675 } 00676 00677 int CubicDiscont2DFECollection::DofForGeometry(int GeomType) const 00678 { 00679 switch (GeomType) 00680 { 00681 case Geometry::POINT: return 0; 00682 case Geometry::SEGMENT: return 0; 00683 case Geometry::TRIANGLE: return 10; 00684 case Geometry::SQUARE: return 16; 00685 default: 00686 mfem_error ("CubicDiscont2DFECollection: unknown geometry type."); 00687 } 00688 return 0; // Make some compilers happy 00689 } 00690 00691 int * CubicDiscont2DFECollection::DofOrderForOrientation(int GeomType, int Or) 00692 const 00693 { 00694 return NULL; 00695 } 00696 00697 00698 const FiniteElement * 00699 LinearNonConf3DFECollection::FiniteElementForGeometry(int GeomType) const 00700 { 00701 switch (GeomType) 00702 { 00703 case Geometry::TRIANGLE: return &TriangleFE; 00704 case Geometry::SQUARE: return &QuadrilateralFE; 00705 case Geometry::TETRAHEDRON: return &TetrahedronFE; 00706 case Geometry::CUBE: return &ParallelepipedFE; 00707 default: 00708 mfem_error ("LinearNonConf3DFECollection: unknown geometry type."); 00709 } 00710 return &TriangleFE; // Make some compilers happy 00711 } 00712 00713 int LinearNonConf3DFECollection::DofForGeometry(int GeomType) const 00714 { 00715 switch (GeomType) 00716 { 00717 case Geometry::POINT: return 0; 00718 case Geometry::SEGMENT: return 0; 00719 case Geometry::TRIANGLE: return 1; 00720 case Geometry::SQUARE: return 1; 00721 case Geometry::TETRAHEDRON: return 0; 00722 case Geometry::CUBE: return 0; 00723 default: 00724 mfem_error ("LinearNonConf3DFECollection: unknown geometry type."); 00725 } 00726 return 0; // Make some compilers happy 00727 } 00728 00729 int * LinearNonConf3DFECollection::DofOrderForOrientation(int GeomType, int Or) 00730 const 00731 { 00732 static int indexes[] = { 0 }; 00733 00734 return indexes; 00735 } 00736 00737 00738 const FiniteElement * 00739 Const3DFECollection::FiniteElementForGeometry(int GeomType) const 00740 { 00741 switch (GeomType) 00742 { 00743 case Geometry::TETRAHEDRON: return &TetrahedronFE; 00744 case Geometry::CUBE: return &ParallelepipedFE; 00745 default: 00746 mfem_error ("Const3DFECollection: unknown geometry type."); 00747 } 00748 return &TetrahedronFE; // Make some compilers happy 00749 } 00750 00751 int Const3DFECollection::DofForGeometry(int GeomType) const 00752 { 00753 switch (GeomType) 00754 { 00755 case Geometry::POINT: return 0; 00756 case Geometry::SEGMENT: return 0; 00757 case Geometry::TRIANGLE: return 0; 00758 case Geometry::TETRAHEDRON: return 1; 00759 case Geometry::SQUARE: return 0; 00760 case Geometry::CUBE: return 1; 00761 default: 00762 mfem_error ("Const3DFECollection: unknown geometry type."); 00763 } 00764 return 0; // Make some compilers happy 00765 } 00766 00767 int * Const3DFECollection::DofOrderForOrientation(int GeomType, int Or) 00768 const 00769 { 00770 return NULL; 00771 } 00772 00773 00774 const FiniteElement * 00775 LinearDiscont3DFECollection::FiniteElementForGeometry(int GeomType) const 00776 { 00777 switch (GeomType) 00778 { 00779 case Geometry::TETRAHEDRON: return &TetrahedronFE; 00780 case Geometry::CUBE: return &ParallelepipedFE; 00781 default: 00782 mfem_error ("LinearDiscont3DFECollection: unknown geometry type."); 00783 } 00784 return &TetrahedronFE; // Make some compilers happy 00785 } 00786 00787 int LinearDiscont3DFECollection::DofForGeometry(int GeomType) const 00788 { 00789 switch (GeomType) 00790 { 00791 case Geometry::POINT: return 0; 00792 case Geometry::SEGMENT: return 0; 00793 case Geometry::TRIANGLE: return 0; 00794 case Geometry::SQUARE: return 0; 00795 case Geometry::TETRAHEDRON: return 4; 00796 case Geometry::CUBE: return 8; 00797 default: 00798 mfem_error ("LinearDiscont3DFECollection: unknown geometry type."); 00799 } 00800 return 0; // Make some compilers happy 00801 } 00802 00803 int * LinearDiscont3DFECollection::DofOrderForOrientation(int GeomType, int Or) 00804 const 00805 { 00806 return NULL; 00807 } 00808 00809 00810 const FiniteElement * 00811 QuadraticDiscont3DFECollection::FiniteElementForGeometry(int GeomType) const 00812 { 00813 switch (GeomType) 00814 { 00815 case Geometry::TETRAHEDRON: return &TetrahedronFE; 00816 case Geometry::CUBE: return &ParallelepipedFE; 00817 default: 00818 mfem_error ("QuadraticDiscont3DFECollection: unknown geometry type."); 00819 } 00820 return &TetrahedronFE; // Make some compilers happy 00821 } 00822 00823 int QuadraticDiscont3DFECollection::DofForGeometry(int GeomType) const 00824 { 00825 switch (GeomType) 00826 { 00827 case Geometry::POINT: return 0; 00828 case Geometry::SEGMENT: return 0; 00829 case Geometry::TRIANGLE: return 0; 00830 case Geometry::SQUARE: return 0; 00831 case Geometry::TETRAHEDRON: return 10; 00832 case Geometry::CUBE: return 27; 00833 default: 00834 mfem_error ("QuadraticDiscont3DFECollection: unknown geometry type."); 00835 } 00836 return 0; // Make some compilers happy 00837 } 00838 00839 int * QuadraticDiscont3DFECollection::DofOrderForOrientation( 00840 int GeomType, int Or) const 00841 { 00842 return NULL; 00843 } 00844 00845 const FiniteElement * 00846 RefinedLinearFECollection::FiniteElementForGeometry(int GeomType) const 00847 { 00848 switch (GeomType) 00849 { 00850 case Geometry::POINT: return &PointFE; 00851 case Geometry::SEGMENT: return &SegmentFE; 00852 case Geometry::TRIANGLE: return &TriangleFE; 00853 case Geometry::SQUARE: return &QuadrilateralFE; 00854 case Geometry::TETRAHEDRON: return &TetrahedronFE; 00855 case Geometry::CUBE: return &ParallelepipedFE; 00856 default: 00857 mfem_error ("RefinedLinearFECollection: unknown geometry type."); 00858 } 00859 return &SegmentFE; // Make some compilers happy 00860 } 00861 00862 int RefinedLinearFECollection::DofForGeometry(int GeomType) const 00863 { 00864 switch (GeomType) 00865 { 00866 case Geometry::POINT: return 1; 00867 case Geometry::SEGMENT: return 1; 00868 case Geometry::TRIANGLE: return 0; 00869 case Geometry::SQUARE: return 1; 00870 case Geometry::TETRAHEDRON: return 0; 00871 case Geometry::CUBE: return 1; 00872 default: 00873 mfem_error ("RefinedLinearFECollection: unknown geometry type."); 00874 } 00875 return 0; // Make some compilers happy 00876 } 00877 00878 int * RefinedLinearFECollection::DofOrderForOrientation(int GeomType, int Or) const 00879 { 00880 static int indexes[] = { 0 }; 00881 00882 return indexes; 00883 } 00884 00885 00886 const FiniteElement * 00887 ND1_3DFECollection::FiniteElementForGeometry(int GeomType) const 00888 { 00889 switch (GeomType) 00890 { 00891 case Geometry::CUBE: return &HexahedronFE; 00892 case Geometry::TETRAHEDRON: return &TetrahedronFE; 00893 default: 00894 mfem_error ("ND1_3DFECollection: unknown geometry type."); 00895 } 00896 return &HexahedronFE; // Make some compilers happy 00897 } 00898 00899 int ND1_3DFECollection::DofForGeometry(int GeomType) const 00900 { 00901 switch (GeomType) 00902 { 00903 case Geometry::POINT: return 0; 00904 case Geometry::SEGMENT: return 1; 00905 case Geometry::TRIANGLE: return 0; 00906 case Geometry::SQUARE: return 0; 00907 case Geometry::TETRAHEDRON: return 0; 00908 case Geometry::CUBE: return 0; 00909 default: 00910 mfem_error ("ND1_3DFECollection: unknown geometry type."); 00911 } 00912 return 0; // Make some compilers happy 00913 } 00914 00915 int * ND1_3DFECollection::DofOrderForOrientation(int GeomType, int Or) 00916 const 00917 { 00918 static int ind_pos[] = { 0 }; 00919 static int ind_neg[] = { -1 }; 00920 00921 if (Or > 0) 00922 return ind_pos; 00923 return ind_neg; 00924 } 00925 00926 00927 const FiniteElement * 00928 RT0_3DFECollection::FiniteElementForGeometry(int GeomType) const 00929 { 00930 switch (GeomType) 00931 { 00932 case Geometry::TRIANGLE: return &TriangleFE; 00933 case Geometry::SQUARE: return &QuadrilateralFE; 00934 case Geometry::CUBE: return &HexahedronFE; 00935 case Geometry::TETRAHEDRON: return &TetrahedronFE; 00936 default: 00937 mfem_error ("RT0_3DFECollection: unknown geometry type."); 00938 } 00939 return &HexahedronFE; // Make some compilers happy 00940 } 00941 00942 int RT0_3DFECollection::DofForGeometry(int GeomType) const 00943 { 00944 switch (GeomType) 00945 { 00946 case Geometry::POINT: return 0; 00947 case Geometry::SEGMENT: return 0; 00948 case Geometry::TRIANGLE: return 1; 00949 case Geometry::SQUARE: return 1; 00950 case Geometry::TETRAHEDRON: return 0; 00951 case Geometry::CUBE: return 0; 00952 default: 00953 mfem_error ("RT0_3DFECollection: unknown geometry type."); 00954 } 00955 return 0; // Make some compilers happy 00956 } 00957 00958 int * RT0_3DFECollection::DofOrderForOrientation(int GeomType, int Or) 00959 const 00960 { 00961 static int ind_pos[] = { 0 }; 00962 static int ind_neg[] = { -1 }; 00963 00964 if ((GeomType == Geometry::TRIANGLE) || (GeomType == Geometry::SQUARE)) 00965 { 00966 if (Or % 2 == 0) 00967 return ind_pos; 00968 return ind_neg; 00969 } 00970 return NULL; 00971 } 00972 00973 const FiniteElement * 00974 RT1_3DFECollection::FiniteElementForGeometry(int GeomType) const 00975 { 00976 switch (GeomType) 00977 { 00978 case Geometry::TRIANGLE: return &TriangleFE; 00979 case Geometry::SQUARE: return &QuadrilateralFE; 00980 case Geometry::CUBE: return &HexahedronFE; 00981 default: 00982 mfem_error ("RT1_3DFECollection: unknown geometry type."); 00983 } 00984 return &HexahedronFE; // Make some compilers happy 00985 } 00986 00987 int RT1_3DFECollection::DofForGeometry(int GeomType) const 00988 { 00989 switch (GeomType) 00990 { 00991 case Geometry::POINT: return 0; 00992 case Geometry::SEGMENT: return 0; 00993 case Geometry::TRIANGLE: return 2; 00994 case Geometry::SQUARE: return 4; 00995 case Geometry::CUBE: return 12; 00996 default: 00997 mfem_error ("RT1_3DFECollection: unknown geometry type."); 00998 } 00999 return 0; // Make some compilers happy 01000 } 01001 01002 int * RT1_3DFECollection::DofOrderForOrientation(int GeomType, int Or) 01003 const 01004 { 01005 if (GeomType == Geometry::SQUARE) 01006 { 01007 static int sq_ind[8][4] = { 01008 {0, 1, 2, 3}, {-1, -3, -2, -4}, 01009 {2, 0, 3, 1}, {-2, -1, -4, -3}, 01010 {3, 2, 1, 0}, {-4, -2, -3, -1}, 01011 {1, 3, 0, 2}, {-3, -4, -1, -2} 01012 }; 01013 01014 return sq_ind[Or]; 01015 } 01016 else 01017 return NULL; 01018 } 01019 01020 01021 H1_FECollection::H1_FECollection(const int p, const int dim) 01022 { 01023 const int pm1 = p - 1, pm2 = pm1 - 1, pm3 = pm2 - 1; 01024 snprintf(h1_name, 32, "H1_%dD_P%d", dim, p); 01025 01026 for (int g = 0; g < Geometry::NumGeom; g++) 01027 { 01028 H1_dof[g] = 0; 01029 H1_Elements[g] = NULL; 01030 } 01031 for (int i = 0; i < 2; i++) 01032 SegDofOrd[i] = NULL; 01033 for (int i = 0; i < 6; i++) 01034 TriDofOrd[i] = NULL; 01035 for (int i = 0; i < 8; i++) 01036 QuadDofOrd[i] = NULL; 01037 01038 H1_dof[Geometry::POINT] = 1; 01039 H1_dof[Geometry::SEGMENT] = pm1; 01040 H1_Elements[Geometry::SEGMENT] = new H1_SegmentElement(p); 01041 01042 if (dim >= 2) 01043 { 01044 H1_dof[Geometry::TRIANGLE] = (pm1*pm2)/2; 01045 H1_dof[Geometry::SQUARE] = pm1*pm1; 01046 H1_Elements[Geometry::TRIANGLE] = new H1_TriangleElement(p); 01047 H1_Elements[Geometry::SQUARE] = new H1_QuadrilateralElement(p); 01048 01049 SegDofOrd[0] = new int[2*pm1]; 01050 SegDofOrd[1] = SegDofOrd[0] + pm1; 01051 for (int i = 0; i < pm1; i++) 01052 { 01053 SegDofOrd[0][i] = i; 01054 SegDofOrd[1][i] = pm2 - i; 01055 } 01056 01057 if (dim >= 3) 01058 { 01059 const int &TriDof = H1_dof[Geometry::TRIANGLE]; 01060 const int &QuadDof = H1_dof[Geometry::SQUARE]; 01061 H1_dof[Geometry::TETRAHEDRON] = (TriDof*pm3)/3; 01062 H1_dof[Geometry::CUBE] = QuadDof*pm1; 01063 H1_Elements[Geometry::TETRAHEDRON] = new H1_TetrahedronElement(p); 01064 H1_Elements[Geometry::CUBE] = new H1_HexahedronElement(p); 01065 01066 TriDofOrd[0] = new int[6*TriDof]; 01067 for (int i = 1; i < 6; i++) 01068 TriDofOrd[i] = TriDofOrd[i-1] + TriDof; 01069 // see Mesh::GetTriOrientation in mesh/mesh.cpp 01070 for (int j = 0; j < pm2; j++) 01071 for (int i = 0; i + j < pm2; i++) 01072 { 01073 int o = TriDof - ((pm1 - j)*(pm2 - j))/2 + i; 01074 int k = pm3 - j - i; 01075 TriDofOrd[0][o] = o; // (0,1,2) 01076 TriDofOrd[1][o] = TriDof - ((pm1-j)*(pm2-j))/2 + k; // (1,0,2) 01077 TriDofOrd[2][o] = TriDof - ((pm1-i)*(pm2-i))/2 + k; // (2,0,1) 01078 TriDofOrd[3][o] = TriDof - ((pm1-k)*(pm2-k))/2 + i; // (2,1,0) 01079 TriDofOrd[4][o] = TriDof - ((pm1-k)*(pm2-k))/2 + j; // (1,2,0) 01080 TriDofOrd[5][o] = TriDof - ((pm1-i)*(pm2-i))/2 + j; // (0,2,1) 01081 } 01082 01083 QuadDofOrd[0] = new int[8*QuadDof]; 01084 for (int i = 1; i < 8; i++) 01085 QuadDofOrd[i] = QuadDofOrd[i-1] + QuadDof; 01086 // see Mesh::GetQuadOrientation in mesh/mesh.cpp 01087 for (int j = 0; j < pm1; j++) 01088 for (int i = 0; i < pm1; i++) 01089 { 01090 int o = i + j*pm1; 01091 QuadDofOrd[0][o] = i + j*pm1; // (0,1,2,3) 01092 QuadDofOrd[1][o] = j + i*pm1; // (0,3,2,1) 01093 QuadDofOrd[2][o] = j + (pm2 - i)*pm1; // (1,2,3,0) 01094 QuadDofOrd[3][o] = (pm2 - i) + j*pm1; // (1,0,3,2) 01095 QuadDofOrd[4][o] = (pm2 - i) + (pm2 - j)*pm1; // (2,3,0,1) 01096 QuadDofOrd[5][o] = (pm2 - j) + (pm2 - i)*pm1; // (2,1,0,3) 01097 QuadDofOrd[6][o] = (pm2 - j) + i*pm1; // (3,0,1,2) 01098 QuadDofOrd[7][o] = i + (pm2 - j)*pm1; // (3,2,1,0) 01099 } 01100 } 01101 } 01102 } 01103 01104 int *H1_FECollection::DofOrderForOrientation(int GeomType, int Or) const 01105 { 01106 if (GeomType == Geometry::SEGMENT) 01107 { 01108 if (Or > 0) 01109 return SegDofOrd[0]; 01110 return SegDofOrd[1]; 01111 } 01112 else if (GeomType == Geometry::TRIANGLE) 01113 { 01114 return TriDofOrd[Or%6]; 01115 } 01116 else if (GeomType == Geometry::SQUARE) 01117 { 01118 return QuadDofOrd[Or%8]; 01119 } 01120 return NULL; 01121 } 01122 01123 H1_FECollection::~H1_FECollection() 01124 { 01125 delete [] SegDofOrd[0]; 01126 delete [] TriDofOrd[0]; 01127 delete [] QuadDofOrd[0]; 01128 for (int g = 0; g < Geometry::NumGeom; g++) 01129 delete H1_Elements[g]; 01130 } 01131 01132 01133 L2_FECollection::L2_FECollection(const int p, const int dim) 01134 { 01135 snprintf(d_name, 32, "L2_%dD_P%d", dim, p); 01136 01137 for (int g = 0; g < Geometry::NumGeom; g++) 01138 L2_Elements[g] = NULL; 01139 01140 if (dim == 1) 01141 { 01142 L2_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p); 01143 } 01144 else if (dim == 2) 01145 { 01146 L2_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p); 01147 L2_Elements[Geometry::SQUARE] = new L2_QuadrilateralElement(p); 01148 } 01149 else if (dim == 3) 01150 { 01151 L2_Elements[Geometry::TETRAHEDRON] = new L2_TetrahedronElement(p); 01152 L2_Elements[Geometry::CUBE] = new L2_HexahedronElement(p); 01153 } 01154 else 01155 { 01156 cerr << "L2_FECollection::L2_FECollection : dim = " 01157 << dim << endl; 01158 mfem_error(); 01159 } 01160 } 01161 01162 L2_FECollection::~L2_FECollection() 01163 { 01164 for (int i = 0; i < Geometry::NumGeom; i++) 01165 delete L2_Elements[i]; 01166 } 01167 01168 01169 RT_FECollection::RT_FECollection(const int p, const int dim) 01170 { 01171 const int pp1 = p + 1, pp2 = p + 2; 01172 01173 snprintf(rt_name, 32, "RT_%dD_P%d", dim, p); 01174 01175 for (int g = 0; g < Geometry::NumGeom; g++) 01176 { 01177 RT_Elements[g] = NULL; 01178 RT_dof[g] = 0; 01179 } 01180 for (int i = 0; i < 2; i++) 01181 SegDofOrd[i] = NULL; 01182 for (int i = 0; i < 6; i++) 01183 TriDofOrd[i] = NULL; 01184 for (int i = 0; i < 8; i++) 01185 QuadDofOrd[i] = NULL; 01186 01187 if (dim == 2) 01188 { 01189 RT_Elements[Geometry::TRIANGLE] = new RT_TriangleElement(p); 01190 RT_dof[Geometry::TRIANGLE] = p*pp1; 01191 01192 RT_Elements[Geometry::SQUARE] = new RT_QuadrilateralElement(p); 01193 RT_dof[Geometry::SQUARE] = 2*p*pp1; 01194 01195 RT_Elements[Geometry::SEGMENT] = new L2_SegmentElement(p); 01196 RT_dof[Geometry::SEGMENT] = pp1; 01197 01198 SegDofOrd[0] = new int[2*pp1]; 01199 SegDofOrd[1] = SegDofOrd[0] + pp1; 01200 for (int i = 0; i <= p; i++) 01201 { 01202 SegDofOrd[0][i] = i; 01203 SegDofOrd[1][i] = -1 - (p - i); 01204 } 01205 } 01206 else if (dim == 3) 01207 { 01208 RT_Elements[Geometry::TETRAHEDRON] = new RT_TetrahedronElement(p); 01209 RT_dof[Geometry::TETRAHEDRON] = p*pp1*pp2/2; 01210 01211 RT_Elements[Geometry::CUBE] = new RT_HexahedronElement(p); 01212 RT_dof[Geometry::CUBE] = 3*p*pp1*pp1; 01213 01214 RT_Elements[Geometry::TRIANGLE] = new L2_TriangleElement(p); 01215 RT_dof[Geometry::TRIANGLE] = pp1*pp2/2; 01216 01217 RT_Elements[Geometry::SQUARE] = new L2_QuadrilateralElement(p); 01218 RT_dof[Geometry::SQUARE] = pp1*pp1; 01219 01220 int TriDof = RT_dof[Geometry::TRIANGLE]; 01221 TriDofOrd[0] = new int[6*TriDof]; 01222 for (int i = 1; i < 6; i++) 01223 TriDofOrd[i] = TriDofOrd[i-1] + TriDof; 01224 // see Mesh::GetTriOrientation in mesh/mesh.cpp, 01225 // the constructor of H1_FECollection 01226 for (int j = 0; j <= p; j++) 01227 for (int i = 0; i + j <= p; i++) 01228 { 01229 int o = TriDof - ((pp2 - j)*(pp1 - j))/2 + i; 01230 int k = p - j - i; 01231 TriDofOrd[0][o] = o; // (0,1,2) 01232 TriDofOrd[1][o] = -1-(TriDof-((pp2-j)*(pp1-j))/2+k); // (1,0,2) 01233 TriDofOrd[2][o] = TriDof-((pp2-i)*(pp1-i))/2+k; // (2,0,1) 01234 TriDofOrd[3][o] = -1-(TriDof-((pp2-k)*(pp1-k))/2+i); // (2,1,0) 01235 TriDofOrd[4][o] = TriDof-((pp2-k)*(pp1-k))/2+j; // (1,2,0) 01236 TriDofOrd[5][o] = -1-(TriDof-((pp2-i)*(pp1-i))/2+j); // (0,2,1) 01237 } 01238 01239 int QuadDof = RT_dof[Geometry::SQUARE]; 01240 QuadDofOrd[0] = new int[8*QuadDof]; 01241 for (int i = 1; i < 8; i++) 01242 QuadDofOrd[i] = QuadDofOrd[i-1] + QuadDof; 01243 // see Mesh::GetQuadOrientation in mesh/mesh.cpp 01244 for (int j = 0; j <= p; j++) 01245 for (int i = 0; i <= p; i++) 01246 { 01247 int o = i + j*pp1; 01248 QuadDofOrd[0][o] = i + j*pp1; // (0,1,2,3) 01249 QuadDofOrd[1][o] = -1 - (j + i*pp1); // (0,3,2,1) 01250 QuadDofOrd[2][o] = j + (p - i)*pp1; // (1,2,3,0) 01251 QuadDofOrd[3][o] = -1 - ((p - i) + j*pp1); // (1,0,3,2) 01252 QuadDofOrd[4][o] = (p - i) + (p - j)*pp1; // (2,3,0,1) 01253 QuadDofOrd[5][o] = -1 - ((p - j) + (p - i)*pp1); // (2,1,0,3) 01254 QuadDofOrd[6][o] = (p - j) + i*pp1; // (3,0,1,2) 01255 QuadDofOrd[7][o] = -1 - (i + (p - j)*pp1); // (3,2,1,0) 01256 } 01257 } 01258 else 01259 { 01260 cerr << "RT_FECollection::RT_FECollection : dim = " << dim << endl; 01261 mfem_error(); 01262 } 01263 } 01264 01265 int *RT_FECollection::DofOrderForOrientation(int GeomType, int Or) const 01266 { 01267 if (GeomType == Geometry::SEGMENT) 01268 { 01269 if (Or > 0) 01270 return SegDofOrd[0]; 01271 return SegDofOrd[1]; 01272 } 01273 else if (GeomType == Geometry::TRIANGLE) 01274 { 01275 return TriDofOrd[Or%6]; 01276 } 01277 else if (GeomType == Geometry::SQUARE) 01278 { 01279 return QuadDofOrd[Or%8]; 01280 } 01281 return NULL; 01282 } 01283 01284 RT_FECollection::~RT_FECollection() 01285 { 01286 delete [] SegDofOrd[0]; 01287 delete [] TriDofOrd[0]; 01288 delete [] QuadDofOrd[0]; 01289 for (int g = 0; g < Geometry::NumGeom; g++) 01290 delete RT_Elements[g]; 01291 } 01292 01293 ND_FECollection::ND_FECollection(const int p, const int dim) 01294 { 01295 const int pm1 = p - 1, pm2 = p - 2; 01296 01297 snprintf(nd_name, 32, "ND_%dD_P%d", dim, p); 01298 01299 for (int g = 0; g < Geometry::NumGeom; g++) 01300 { 01301 ND_Elements[g] = NULL; 01302 ND_dof[g] = 0; 01303 } 01304 for (int i = 0; i < 2; i++) 01305 SegDofOrd[i] = NULL; 01306 for (int i = 0; i < 6; i++) 01307 TriDofOrd[i] = NULL; 01308 for (int i = 0; i < 8; i++) 01309 QuadDofOrd[i] = NULL; 01310 01311 if (dim == 2 || dim == 3) 01312 { 01313 ND_Elements[Geometry::SQUARE] = new ND_QuadrilateralElement(p); 01314 ND_dof[Geometry::SQUARE] = 2*p*pm1; 01315 01316 ND_Elements[Geometry::TRIANGLE] = new ND_TriangleElement(p); 01317 ND_dof[Geometry::TRIANGLE] = p*pm1; 01318 01319 // ND_Elements[Geometry::SEGMENT] = NULL; 01320 ND_dof[Geometry::SEGMENT] = p; 01321 01322 SegDofOrd[0] = new int[2*p]; 01323 SegDofOrd[1] = SegDofOrd[0] + p; 01324 for (int i = 0; i < p; i++) 01325 { 01326 SegDofOrd[0][i] = i; 01327 SegDofOrd[1][i] = -1 - (pm1 - i); 01328 } 01329 } 01330 else 01331 { 01332 mfem_error("ND_FECollection::ND_FECollection : dim != 2 or 3"); 01333 } 01334 01335 if (dim == 3) 01336 { 01337 ND_Elements[Geometry::CUBE] = new ND_HexahedronElement(p); 01338 ND_dof[Geometry::CUBE] = 3*p*pm1*pm1; 01339 01340 ND_Elements[Geometry::TETRAHEDRON] = new ND_TetrahedronElement(p); 01341 ND_dof[Geometry::TETRAHEDRON] = p*pm1*pm2/2; 01342 01343 int QuadDof = ND_dof[Geometry::SQUARE]; 01344 QuadDofOrd[0] = new int[8*QuadDof]; 01345 for (int i = 1; i < 8; i++) 01346 QuadDofOrd[i] = QuadDofOrd[i-1] + QuadDof; 01347 // see Mesh::GetQuadOrientation in mesh/mesh.cpp 01348 for (int j = 0; j < pm1; j++) 01349 for (int i = 0; i < p; i++) 01350 { 01351 int d1 = i + j*p; // x-component 01352 int d2 = p*pm1 + j + i*pm1; // y-component 01353 // (0,1,2,3) 01354 QuadDofOrd[0][d1] = d1; 01355 QuadDofOrd[0][d2] = d2; 01356 // (0,3,2,1) 01357 QuadDofOrd[1][d1] = d2; 01358 QuadDofOrd[1][d2] = d1; 01359 // (1,2,3,0) 01360 // QuadDofOrd[2][d1] = p*pm1 + (pm2 - j) + i*pm1; 01361 // QuadDofOrd[2][d2] = -1 - ((pm1 - i) + j*p); 01362 QuadDofOrd[2][d1] = -1 - (p*pm1 + j + (pm1 - i)*pm1); 01363 QuadDofOrd[2][d2] = i + (pm2 - j)*p; 01364 // (1,0,3,2) 01365 QuadDofOrd[3][d1] = -1 - ((pm1 - i) + j*p); 01366 QuadDofOrd[3][d2] = p*pm1 + (pm2 - j) + i*pm1; 01367 // (2,3,0,1) 01368 QuadDofOrd[4][d1] = -1 - ((pm1 - i) + (pm2 - j)*p); 01369 QuadDofOrd[4][d2] = -1 - (p*pm1 + (pm2 - j) + (pm1 - i)*pm1); 01370 // (2,1,0,3) 01371 QuadDofOrd[5][d1] = -1 - (p*pm1 + (pm2 - j) + (pm1 - i)*pm1); 01372 QuadDofOrd[5][d2] = -1 - ((pm1 - i) + (pm2 - j)*p); 01373 // (3,0,1,2) 01374 // QuadDofOrd[6][d1] = -1 - (p*pm1 + j + (pm1 - i)*pm1); 01375 // QuadDofOrd[6][d2] = i + (pm2 - j)*p; 01376 QuadDofOrd[6][d1] = p*pm1 + (pm2 - j) + i*pm1; 01377 QuadDofOrd[6][d2] = -1 - ((pm1 - i) + j*p); 01378 // (3,2,1,0) 01379 QuadDofOrd[7][d1] = i + (pm2 - j)*p; 01380 QuadDofOrd[7][d2] = -1 - (p*pm1 + j + (pm1 - i)*pm1); 01381 } 01382 01383 int TriDof = ND_dof[Geometry::TRIANGLE]; 01384 TriDofOrd[0] = new int[6*TriDof]; 01385 for (int i = 1; i < 6; i++) 01386 TriDofOrd[i] = TriDofOrd[i-1] + TriDof; 01387 // see Mesh::GetTriOrientation in mesh/mesh.cpp, 01388 // the constructor of H1_FECollection 01389 for (int j = 0; j <= pm2; j++) 01390 for (int i = 0; i + j <= pm2; i++) 01391 { 01392 int k1 = p*pm1 - (p - j)*(pm1 - j) + 2*i; 01393 int k2 = p*pm1 - (p - i)*(pm1 - i) + 2*j; 01394 // (0,1,2) 01395 TriDofOrd[0][k1 ] = k1; 01396 TriDofOrd[0][k1+1] = k1 + 1; 01397 // (0,2,1) 01398 TriDofOrd[5][k1 ] = k2 + 1; 01399 TriDofOrd[5][k1+1] = k2; 01400 01401 // The other orientations can not be supported with the current 01402 // interface. The method Mesh::ReorientTetMesh will ensure that 01403 // only orientations 0 and 5 are generated. 01404 } 01405 } 01406 } 01407 01408 int *ND_FECollection::DofOrderForOrientation(int GeomType, int Or) const 01409 { 01410 if (GeomType == Geometry::SEGMENT) 01411 { 01412 if (Or > 0) 01413 return SegDofOrd[0]; 01414 return SegDofOrd[1]; 01415 } 01416 else if (GeomType == Geometry::TRIANGLE) 01417 { 01418 if (Or != 0 && Or != 5) 01419 { 01420 cerr << 01421 "ND_FECollection::DofOrderForOrientation :\n" 01422 " triangle face orientation " << Or << " is not supported!\n" 01423 " Use Mesh::ReorientTetMesh to fix it." << endl; 01424 mfem_error(); 01425 } 01426 return TriDofOrd[Or%6]; 01427 } 01428 else if (GeomType == Geometry::SQUARE) 01429 { 01430 return QuadDofOrd[Or%8]; 01431 } 01432 return NULL; 01433 } 01434 01435 ND_FECollection::~ND_FECollection() 01436 { 01437 delete [] SegDofOrd[0]; 01438 delete [] TriDofOrd[0]; 01439 delete [] QuadDofOrd[0]; 01440 for (int g = 0; g < Geometry::NumGeom; g++) 01441 delete ND_Elements[g]; 01442 } 01443 01444 01445 Local_FECollection::Local_FECollection(const char *fe_name) 01446 { 01447 snprintf(d_name, 32, "Local_%s", fe_name); 01448 01449 Local_Element = NULL; 01450 01451 if (!strcmp(fe_name, "BiCubic2DFiniteElement") || 01452 !strcmp(fe_name, "Quad_Q3")) 01453 { 01454 GeomType = Geometry::SQUARE; 01455 Local_Element = new BiCubic2DFiniteElement; 01456 } 01457 else if (!strcmp(fe_name, "Nedelec1HexFiniteElement") || 01458 !strcmp(fe_name, "Hex_ND1")) 01459 { 01460 GeomType = Geometry::CUBE; 01461 Local_Element = new Nedelec1HexFiniteElement; 01462 } 01463 else 01464 { 01465 cerr << "Local_FECollection::Local_FECollection : fe_name = " 01466 << fe_name << endl; 01467 mfem_error(); 01468 } 01469 } 01470 01471 01472 void NURBSFECollection::Allocate(int Order) 01473 { 01474 SegmentFE = new NURBS1DFiniteElement(Order); 01475 QuadrilateralFE = new NURBS2DFiniteElement(Order); 01476 ParallelepipedFE = new NURBS3DFiniteElement(Order); 01477 01478 sprintf(name, "NURBS%i", Order); 01479 } 01480 01481 void NURBSFECollection::Deallocate() 01482 { 01483 delete ParallelepipedFE; 01484 delete QuadrilateralFE; 01485 delete SegmentFE; 01486 } 01487 01488 const FiniteElement * 01489 NURBSFECollection::FiniteElementForGeometry(int GeomType) const 01490 { 01491 switch (GeomType) 01492 { 01493 case Geometry::SEGMENT: return SegmentFE; 01494 case Geometry::SQUARE: return QuadrilateralFE; 01495 case Geometry::CUBE: return ParallelepipedFE; 01496 default: 01497 mfem_error ("NURBSFECollection: unknown geometry type."); 01498 } 01499 return SegmentFE; // Make some compilers happy 01500 } 01501 01502 int NURBSFECollection::DofForGeometry(int GeomType) const 01503 { 01504 mfem_error("NURBSFECollection::DofForGeometry"); 01505 return 0; // Make some compilers happy 01506 } 01507 01508 int *NURBSFECollection::DofOrderForOrientation(int GeomType, int Or) const 01509 { 01510 mfem_error("NURBSFECollection::DofOrderForOrientation"); 01511 return NULL; 01512 }