MFEM v2.0
fe_coll.cpp
Go to the documentation of this file.
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 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines