MFEM: fe.hpp Source File

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 #ifndef MFEM_FE
00013 #define MFEM_FE
00014 
00015 // Base and derived classes for finite elements
00016 
00018 class FunctionSpace {
00019 public:
00020    enum {
00021       Pk, // polynomials of order k
00022       Qk, // tensor products of polynomials of order k
00023       rQk // refined tensor products of polynomials of order k
00024    };
00025 };
00026 
00027 class ElementTransformation;
00028 class Coefficient;
00029 class VectorCoefficient;
00030 class KnotVector;
00031 
00033 class FiniteElement
00034 {
00035 protected:
00036    int Dim, GeomType, Dof, Order, FuncSpace, RangeType;
00037    IntegrationRule Nodes;
00038 
00039 public:
00040    enum { SCALAR, VECTOR };
00041 
00048    FiniteElement(int D, int G, int Do, int O, int F = FunctionSpace::Pk);
00049 
00051    int GetDim() const { return Dim; }
00052 
00054    int GetGeomType() const { return GeomType; }
00055 
00057    int GetDof() const { return Dof; }
00058 
00060    int GetOrder() const { return Order; }
00061 
00063    int Space() const { return FuncSpace; }
00064 
00065    int GetRangeType() const { return RangeType; }
00066 
00070    virtual void CalcShape(const IntegrationPoint &ip,
00071                           Vector &shape) const = 0;
00072 
00077    virtual void CalcDShape(const IntegrationPoint &ip,
00078                            DenseMatrix &dshape) const = 0;
00079    const IntegrationRule & GetNodes() const { return Nodes; }
00080 
00081    // virtual functions for finite elements on vector spaces
00082 
00087    virtual void CalcVShape(const IntegrationPoint &ip,
00088                            DenseMatrix &shape) const;
00089 
00090    virtual void CalcVShape(ElementTransformation &Trans,
00091                            DenseMatrix &shape) const;
00092 
00096    virtual void CalcDivShape(const IntegrationPoint &ip,
00097                              Vector &divshape) const;
00098 
00103    virtual void CalcCurlShape(const IntegrationPoint &ip,
00104                               DenseMatrix &curl_shape) const;
00105 
00106    virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const;
00107 
00110    virtual void CalcHessian (const IntegrationPoint &ip,
00111                              DenseMatrix &h) const;
00112 
00116    virtual void GetLocalInterpolation (ElementTransformation &Trans,
00117                                        DenseMatrix &I) const;
00118 
00122    virtual void Project (Coefficient &coeff,
00123                          ElementTransformation &Trans, Vector &dofs) const;
00124 
00128    virtual void Project (VectorCoefficient &vc,
00129                          ElementTransformation &Trans, Vector &dofs) const;
00130 
00133    virtual void ProjectDelta(int vertex, Vector &dofs) const;
00134 
00138    virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
00139                         DenseMatrix &I) const;
00140 
00144    virtual void ProjectGrad(const FiniteElement &fe,
00145                             ElementTransformation &Trans,
00146                             DenseMatrix &grad) const;
00147 
00151    virtual void ProjectCurl(const FiniteElement &fe,
00152                             ElementTransformation &Trans,
00153                             DenseMatrix &curl) const;
00154 
00158    virtual void ProjectDiv(const FiniteElement &fe,
00159                            ElementTransformation &Trans,
00160                            DenseMatrix &div) const;
00161 
00162    virtual ~FiniteElement () { }
00163 };
00164 
00165 class NodalFiniteElement : public FiniteElement
00166 {
00167 protected:
00168    void NodalLocalInterpolation (ElementTransformation &Trans,
00169                                  DenseMatrix &I,
00170                                  const NodalFiniteElement &fine_fe) const;
00171 
00172 #ifndef MFEM_USE_OPENMP
00173    mutable Vector c_shape;
00174 #endif
00175 
00176 public:
00177    NodalFiniteElement(int D, int G, int Do, int O,
00178                       int F = FunctionSpace::Pk) :
00179 #ifdef MFEM_USE_OPENMP
00180       FiniteElement(D, G, Do, O, F)
00181 #else
00182       FiniteElement(D, G, Do, O, F), c_shape(Do)
00183 #endif
00184    { }
00185 
00186    virtual void GetLocalInterpolation (ElementTransformation &Trans,
00187                                        DenseMatrix &I) const
00188    { NodalLocalInterpolation (Trans, I, *this); }
00189 
00190    virtual void Project (Coefficient &coeff,
00191                          ElementTransformation &Trans, Vector &dofs) const;
00192 
00193    virtual void Project (VectorCoefficient &vc,
00194                          ElementTransformation &Trans, Vector &dofs) const;
00195 
00196    virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
00197                         DenseMatrix &I) const;
00198 
00199    virtual void ProjectGrad(const FiniteElement &fe,
00200                             ElementTransformation &Trans,
00201                             DenseMatrix &grad) const;
00202 
00203    virtual void ProjectDiv(const FiniteElement &fe,
00204                            ElementTransformation &Trans,
00205                            DenseMatrix &div) const;
00206 };
00207 
00208 class VectorFiniteElement : public FiniteElement
00209 {
00210    // Hide the scalar functions CalcShape and CalcDShape.
00211 private:
00213    virtual void CalcShape(const IntegrationPoint &ip,
00214                           Vector &shape) const;
00215 
00217    virtual void CalcDShape(const IntegrationPoint &ip,
00218                            DenseMatrix &dshape) const;
00219 
00220 protected:
00221 #ifndef MFEM_USE_OPENMP
00222    mutable DenseMatrix Jinv;
00223    mutable DenseMatrix vshape;
00224 #endif
00225 
00226    void CalcVShape_RT(ElementTransformation &Trans,
00227                       DenseMatrix &shape) const;
00228 
00229    void CalcVShape_ND(ElementTransformation &Trans,
00230                       DenseMatrix &shape) const;
00231 
00232    void Project_RT(const double *nk, const Array<int> &d2n,
00233                    VectorCoefficient &vc, ElementTransformation &Trans,
00234                    Vector &dofs) const;
00235 
00236    void Project_RT(const double *nk, const Array<int> &d2n,
00237                    const FiniteElement &fe, ElementTransformation &Trans,
00238                    DenseMatrix &I) const;
00239 
00240    // rotated gradient in 2D
00241    void ProjectGrad_RT(const double *nk, const Array<int> &d2n,
00242                        const FiniteElement &fe, ElementTransformation &Trans,
00243                        DenseMatrix &grad) const;
00244 
00245    void ProjectCurl_RT(const double *nk, const Array<int> &d2n,
00246                        const FiniteElement &fe, ElementTransformation &Trans,
00247                        DenseMatrix &curl) const;
00248 
00249    void Project_ND(const double *tk, const Array<int> &d2t,
00250                    VectorCoefficient &vc, ElementTransformation &Trans,
00251                    Vector &dofs) const;
00252 
00253    void Project_ND(const double *tk, const Array<int> &d2t,
00254                    const FiniteElement &fe, ElementTransformation &Trans,
00255                    DenseMatrix &I) const;
00256 
00257    void ProjectGrad_ND(const double *tk, const Array<int> &d2t,
00258                        const FiniteElement &fe, ElementTransformation &Trans,
00259                        DenseMatrix &grad) const;
00260 
00261    void LocalInterpolation_RT(const double *nk, const Array<int> &d2n,
00262                               ElementTransformation &Trans,
00263                               DenseMatrix &I) const;
00264 
00265    void LocalInterpolation_ND(const double *tk, const Array<int> &d2t,
00266                               ElementTransformation &Trans,
00267                               DenseMatrix &I) const;
00268 
00269 public:
00270    VectorFiniteElement (int D, int G, int Do, int O,
00271                         int F = FunctionSpace::Pk) :
00272 #ifdef MFEM_USE_OPENMP
00273       FiniteElement(D, G, Do, O, F)
00274 #else
00275       FiniteElement(D, G, Do, O, F), Jinv(D), vshape(Do, D)
00276 #endif
00277    { RangeType = VECTOR; }
00278 };
00279 
00280 class PointFiniteElement : public NodalFiniteElement
00281 {
00282 public:
00283    PointFiniteElement();
00284 
00285    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00286 
00287    virtual void CalcDShape(const IntegrationPoint &ip,
00288                            DenseMatrix &dshape) const;
00289 };
00290 
00292 class Linear1DFiniteElement : public NodalFiniteElement
00293 {
00294 public:
00296    Linear1DFiniteElement();
00297 
00301    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00302 
00307    virtual void CalcDShape(const IntegrationPoint &ip,
00308                            DenseMatrix &dshape) const;
00309 };
00310 
00312 class Linear2DFiniteElement : public NodalFiniteElement
00313 {
00314 public:
00316    Linear2DFiniteElement();
00317 
00321    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00322 
00327    virtual void CalcDShape(const IntegrationPoint &ip,
00328                            DenseMatrix &dshape) const;
00329    virtual void ProjectDelta(int vertex, Vector &dofs) const
00330    { dofs = 0.0; dofs(vertex) = 1.0; }
00331 };
00332 
00334 class BiLinear2DFiniteElement : public NodalFiniteElement
00335 {
00336 public:
00338    BiLinear2DFiniteElement();
00339 
00343    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00344 
00349    virtual void CalcDShape(const IntegrationPoint &ip,
00350                            DenseMatrix &dshape) const;
00351    virtual void CalcHessian (const IntegrationPoint &ip,
00352                              DenseMatrix &h) const;
00353    virtual void ProjectDelta(int vertex, Vector &dofs) const
00354    { dofs = 0.0; dofs(vertex) = 1.0; }
00355 //      { dofs = 1.0; }
00356 };
00357 
00359 class GaussLinear2DFiniteElement : public NodalFiniteElement
00360 {
00361 public:
00362    GaussLinear2DFiniteElement();
00363    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00364    virtual void CalcDShape(const IntegrationPoint &ip,
00365                            DenseMatrix &dshape) const;
00366    virtual void ProjectDelta(int vertex, Vector &dofs) const;
00367 };
00368 
00370 class GaussBiLinear2DFiniteElement : public NodalFiniteElement
00371 {
00372 private:
00373    static const double p[2];
00374 
00375 public:
00376    GaussBiLinear2DFiniteElement();
00377    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00378    virtual void CalcDShape(const IntegrationPoint &ip,
00379                            DenseMatrix &dshape) const;
00380    virtual void ProjectDelta(int vertex, Vector &dofs) const;
00381 };
00382 
00383 class P1OnQuadFiniteElement : public NodalFiniteElement
00384 {
00385 public:
00386    P1OnQuadFiniteElement();
00387    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00388    virtual void CalcDShape(const IntegrationPoint &ip,
00389                            DenseMatrix &dshape) const;
00390    virtual void ProjectDelta(int vertex, Vector &dofs) const
00391    { dofs = 1.0; }
00392 };
00393 
00395 class Quad1DFiniteElement : public NodalFiniteElement
00396 {
00397 public:
00399    Quad1DFiniteElement();
00400 
00404    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00405 
00410    virtual void CalcDShape(const IntegrationPoint &ip,
00411                            DenseMatrix &dshape) const;
00412 };
00413 
00414 class QuadPos1DFiniteElement : public FiniteElement
00415 {
00416 public:
00417    QuadPos1DFiniteElement();
00418    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00419    virtual void CalcDShape(const IntegrationPoint &ip,
00420                            DenseMatrix &dshape) const;
00421 };
00422 
00424 class Quad2DFiniteElement : public NodalFiniteElement
00425 {
00426 public:
00428    Quad2DFiniteElement();
00429 
00433    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00434 
00439    virtual void CalcDShape(const IntegrationPoint &ip,
00440                            DenseMatrix &dshape) const;
00441 
00442    virtual void CalcHessian (const IntegrationPoint &ip,
00443                              DenseMatrix &h) const;
00444    virtual void ProjectDelta(int vertex, Vector &dofs) const;
00445 };
00446 
00448 class GaussQuad2DFiniteElement : public NodalFiniteElement
00449 {
00450 private:
00451    static const double p[2];
00452    DenseMatrix A;
00453    mutable DenseMatrix D;
00454    mutable Vector pol;
00455 public:
00456    GaussQuad2DFiniteElement();
00457    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00458    virtual void CalcDShape(const IntegrationPoint &ip,
00459                            DenseMatrix &dshape) const;
00460    // virtual void ProjectDelta(int vertex, Vector &dofs) const;
00461 };
00462 
00464 class BiQuad2DFiniteElement : public NodalFiniteElement
00465 {
00466 public:
00468    BiQuad2DFiniteElement();
00469 
00473    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00474 
00479    virtual void CalcDShape(const IntegrationPoint &ip,
00480                            DenseMatrix &dshape) const;
00481    virtual void ProjectDelta(int vertex, Vector &dofs) const;
00482 };
00483 
00484 class BiQuadPos2DFiniteElement : public FiniteElement
00485 {
00486 public:
00487    BiQuadPos2DFiniteElement();
00488    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00489    virtual void CalcDShape(const IntegrationPoint &ip,
00490                            DenseMatrix &dshape) const;
00491    using FiniteElement::Project;
00492    virtual void Project(Coefficient &coeff, ElementTransformation &Trans,
00493                         Vector &dofs) const;
00494    virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans,
00495                         Vector &dofs) const;
00496    virtual void ProjectDelta(int vertex, Vector &dofs) const
00497    { dofs = 0.; dofs(vertex) = 1.; }
00498 };
00499 
00501 class GaussBiQuad2DFiniteElement : public NodalFiniteElement
00502 {
00503 public:
00504    GaussBiQuad2DFiniteElement();
00505    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00506    virtual void CalcDShape(const IntegrationPoint &ip,
00507                            DenseMatrix &dshape) const;
00508 //   virtual void ProjectDelta(int vertex, Vector &dofs) const { dofs = 1.; }
00509 };
00510 
00511 class BiCubic2DFiniteElement : public NodalFiniteElement
00512 {
00513 public:
00514    BiCubic2DFiniteElement();
00515    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00516    virtual void CalcDShape(const IntegrationPoint &ip,
00517                            DenseMatrix &dshape) const;
00518    virtual void CalcHessian (const IntegrationPoint &ip,
00519                              DenseMatrix &h) const;
00520 };
00521 
00522 class Cubic1DFiniteElement : public NodalFiniteElement
00523 {
00524 public:
00525    Cubic1DFiniteElement();
00526 
00527    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00528 
00529    virtual void CalcDShape(const IntegrationPoint &ip,
00530                            DenseMatrix &dshape) const;
00531 };
00532 
00533 class Cubic2DFiniteElement : public NodalFiniteElement
00534 {
00535 public:
00536    Cubic2DFiniteElement();
00537 
00538    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00539 
00540    virtual void CalcDShape(const IntegrationPoint &ip,
00541                            DenseMatrix &dshape) const;
00542 
00543    virtual void CalcHessian (const IntegrationPoint &ip,
00544                              DenseMatrix &h) const;
00545 };
00546 
00548 class Cubic3DFiniteElement : public NodalFiniteElement
00549 {
00550 public:
00552    Cubic3DFiniteElement();
00553 
00554    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00555 
00556    virtual void CalcDShape(const IntegrationPoint &ip,
00557                            DenseMatrix &dshape) const;
00558 };
00559 
00561 class P0TriangleFiniteElement : public NodalFiniteElement
00562 {
00563 public:
00565    P0TriangleFiniteElement();
00566 
00568    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00569 
00571    virtual void CalcDShape(const IntegrationPoint &ip,
00572                            DenseMatrix &dshape) const;
00573    virtual void ProjectDelta(int vertex, Vector &dofs) const
00574    { dofs(0) = 1.0; }
00575 };
00576 
00577 
00578 class P0QuadFiniteElement : public NodalFiniteElement
00579 {
00580 public:
00581    P0QuadFiniteElement();
00582    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00583    virtual void CalcDShape(const IntegrationPoint &ip,
00584                            DenseMatrix &dshape) const;
00585    virtual void ProjectDelta(int vertex, Vector &dofs) const
00586    { dofs(0) = 1.0; }
00587 };
00588 
00589 
00591 class Linear3DFiniteElement : public NodalFiniteElement
00592 {
00593 public:
00595    Linear3DFiniteElement();
00596 
00600    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00601 
00606    virtual void CalcDShape(const IntegrationPoint &ip,
00607                            DenseMatrix &dshape) const;
00608 
00609    virtual void ProjectDelta(int vertex, Vector &dofs) const
00610    { dofs = 0.0; dofs(vertex) = 1.0; }
00611 
00612    virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const;
00613 };
00614 
00616 class Quadratic3DFiniteElement : public NodalFiniteElement
00617 {
00618 public:
00620    Quadratic3DFiniteElement();
00621 
00622    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00623 
00624    virtual void CalcDShape(const IntegrationPoint &ip,
00625                            DenseMatrix &dshape) const;
00626 };
00627 
00629 class TriLinear3DFiniteElement : public NodalFiniteElement
00630 {
00631 public:
00633    TriLinear3DFiniteElement();
00634 
00638    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00639 
00644    virtual void CalcDShape(const IntegrationPoint &ip,
00645                            DenseMatrix &dshape) const;
00646 
00647    virtual void ProjectDelta(int vertex, Vector &dofs) const
00648    { dofs = 0.0; dofs(vertex) = 1.0; }
00649 };
00650 
00652 class CrouzeixRaviartFiniteElement : public NodalFiniteElement
00653 {
00654 public:
00655    CrouzeixRaviartFiniteElement();
00656    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00657    virtual void CalcDShape(const IntegrationPoint &ip,
00658                            DenseMatrix &dshape) const;
00659    virtual void ProjectDelta(int vertex, Vector &dofs) const
00660    { dofs = 1.0; }
00661 };
00662 
00664 class CrouzeixRaviartQuadFiniteElement : public NodalFiniteElement
00665 {
00666 public:
00667    CrouzeixRaviartQuadFiniteElement();
00668    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00669    virtual void CalcDShape(const IntegrationPoint &ip,
00670                            DenseMatrix &dshape) const;
00671 };
00672 
00673 class P0SegmentFiniteElement : public NodalFiniteElement
00674 {
00675 public:
00676    P0SegmentFiniteElement(int Ord = 0);
00677    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00678    virtual void CalcDShape(const IntegrationPoint &ip,
00679                            DenseMatrix &dshape) const;
00680 };
00681 
00682 class RT0TriangleFiniteElement : public VectorFiniteElement
00683 {
00684 private:
00685    static const double nk[3][2];
00686 
00687 public:
00688    RT0TriangleFiniteElement();
00689 
00690    virtual void CalcVShape(const IntegrationPoint &ip,
00691                            DenseMatrix &shape) const;
00692 
00693    virtual void CalcVShape(ElementTransformation &Trans,
00694                            DenseMatrix &shape) const
00695    { CalcVShape_RT(Trans, shape); };
00696 
00697    virtual void CalcDivShape(const IntegrationPoint &ip,
00698                              Vector &divshape) const;
00699 
00700    virtual void GetLocalInterpolation (ElementTransformation &Trans,
00701                                        DenseMatrix &I) const;
00702 
00703    using FiniteElement::Project;
00704 
00705    virtual void Project (VectorCoefficient &vc,
00706                          ElementTransformation &Trans, Vector &dofs) const;
00707 };
00708 
00709 class RT0QuadFiniteElement : public VectorFiniteElement
00710 {
00711 private:
00712    static const double nk[4][2];
00713 
00714 public:
00715    RT0QuadFiniteElement();
00716 
00717    virtual void CalcVShape(const IntegrationPoint &ip,
00718                            DenseMatrix &shape) const;
00719 
00720    virtual void CalcVShape(ElementTransformation &Trans,
00721                            DenseMatrix &shape) const
00722    { CalcVShape_RT(Trans, shape); };
00723 
00724    virtual void CalcDivShape(const IntegrationPoint &ip,
00725                              Vector &divshape) const;
00726 
00727    virtual void GetLocalInterpolation (ElementTransformation &Trans,
00728                                        DenseMatrix &I) const;
00729 
00730    using FiniteElement::Project;
00731 
00732    virtual void Project (VectorCoefficient &vc,
00733                          ElementTransformation &Trans, Vector &dofs) const;
00734 };
00735 
00736 class RT1TriangleFiniteElement : public VectorFiniteElement
00737 {
00738 private:
00739    static const double nk[8][2];
00740 
00741 public:
00742    RT1TriangleFiniteElement();
00743 
00744    virtual void CalcVShape(const IntegrationPoint &ip,
00745                            DenseMatrix &shape) const;
00746 
00747    virtual void CalcVShape(ElementTransformation &Trans,
00748                            DenseMatrix &shape) const
00749    { CalcVShape_RT(Trans, shape); };
00750 
00751    virtual void CalcDivShape(const IntegrationPoint &ip,
00752                              Vector &divshape) const;
00753 
00754    virtual void GetLocalInterpolation (ElementTransformation &Trans,
00755                                        DenseMatrix &I) const;
00756 
00757    using FiniteElement::Project;
00758 
00759    virtual void Project (VectorCoefficient &vc,
00760                          ElementTransformation &Trans, Vector &dofs) const;
00761 };
00762 
00763 class RT1QuadFiniteElement : public VectorFiniteElement
00764 {
00765 private:
00766    static const double nk[12][2];
00767 
00768 public:
00769    RT1QuadFiniteElement();
00770 
00771    virtual void CalcVShape(const IntegrationPoint &ip,
00772                            DenseMatrix &shape) const;
00773 
00774    virtual void CalcVShape(ElementTransformation &Trans,
00775                            DenseMatrix &shape) const
00776    { CalcVShape_RT(Trans, shape); };
00777 
00778    virtual void CalcDivShape(const IntegrationPoint &ip,
00779                              Vector &divshape) const;
00780 
00781    virtual void GetLocalInterpolation (ElementTransformation &Trans,
00782                                        DenseMatrix &I) const;
00783 
00784    using FiniteElement::Project;
00785 
00786    virtual void Project (VectorCoefficient &vc,
00787                          ElementTransformation &Trans, Vector &dofs) const;
00788 };
00789 
00790 class RT2TriangleFiniteElement : public VectorFiniteElement
00791 {
00792 private:
00793    static const double M[15][15];
00794 public:
00795    RT2TriangleFiniteElement();
00796 
00797    virtual void CalcVShape(const IntegrationPoint &ip,
00798                            DenseMatrix &shape) const;
00799 
00800    virtual void CalcVShape(ElementTransformation &Trans,
00801                            DenseMatrix &shape) const
00802    { CalcVShape_RT(Trans, shape); };
00803 
00804    virtual void CalcDivShape(const IntegrationPoint &ip,
00805                              Vector &divshape) const;
00806 };
00807 
00808 class RT2QuadFiniteElement : public VectorFiniteElement
00809 {
00810 private:
00811    static const double nk[24][2];
00812    static const double pt[4];
00813    static const double dpt[3];
00814 
00815 public:
00816    RT2QuadFiniteElement();
00817 
00818    virtual void CalcVShape(const IntegrationPoint &ip,
00819                            DenseMatrix &shape) const;
00820 
00821    virtual void CalcVShape(ElementTransformation &Trans,
00822                            DenseMatrix &shape) const
00823    { CalcVShape_RT(Trans, shape); };
00824 
00825    virtual void CalcDivShape(const IntegrationPoint &ip,
00826                              Vector &divshape) const;
00827 
00828    virtual void GetLocalInterpolation (ElementTransformation &Trans,
00829                                        DenseMatrix &I) const;
00830 
00831    using FiniteElement::Project;
00832 
00833    virtual void Project (VectorCoefficient &vc,
00834                          ElementTransformation &Trans, Vector &dofs) const;
00835 };
00836 
00838 class P1SegmentFiniteElement : public NodalFiniteElement
00839 {
00840 public:
00841    P1SegmentFiniteElement();
00842    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00843    virtual void CalcDShape(const IntegrationPoint &ip,
00844                            DenseMatrix &dshape) const;
00845 };
00846 
00848 class P2SegmentFiniteElement : public NodalFiniteElement
00849 {
00850 public:
00851    P2SegmentFiniteElement();
00852    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00853    virtual void CalcDShape(const IntegrationPoint &ip,
00854                            DenseMatrix &dshape) const;
00855 };
00856 
00857 class Lagrange1DFiniteElement : public NodalFiniteElement
00858 {
00859 private:
00860    Vector rwk;
00861 #ifndef MFEM_USE_OPENMP
00862    mutable Vector rxxk;
00863 #endif
00864 public:
00865    Lagrange1DFiniteElement (int degree);
00866    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00867    virtual void CalcDShape(const IntegrationPoint &ip,
00868                            DenseMatrix &dshape) const;
00869 };
00870 
00871 class P1TetNonConfFiniteElement : public NodalFiniteElement
00872 {
00873 public:
00874    P1TetNonConfFiniteElement();
00875    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00876    virtual void CalcDShape(const IntegrationPoint &ip,
00877                            DenseMatrix &dshape) const;
00878 };
00879 
00880 class P0TetFiniteElement : public NodalFiniteElement
00881 {
00882 public:
00883    P0TetFiniteElement ();
00884    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00885    virtual void CalcDShape(const IntegrationPoint &ip,
00886                            DenseMatrix &dshape) const;
00887    virtual void ProjectDelta(int vertex, Vector &dofs) const
00888    { dofs(0) = 1.0; }
00889 };
00890 
00891 class P0HexFiniteElement : public NodalFiniteElement
00892 {
00893 public:
00894    P0HexFiniteElement ();
00895    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00896    virtual void CalcDShape(const IntegrationPoint &ip,
00897                            DenseMatrix &dshape) const;
00898    virtual void ProjectDelta(int vertex, Vector &dofs) const
00899    { dofs(0) = 1.0; }
00900 };
00901 
00903 class LagrangeHexFiniteElement : public NodalFiniteElement
00904 {
00905 private:
00906    Lagrange1DFiniteElement * fe1d;
00907    int dof1d;
00908    int *I, *J, *K;
00909 #ifndef MFEM_USE_OPENMP
00910    mutable Vector shape1dx, shape1dy, shape1dz;
00911    mutable DenseMatrix dshape1dx, dshape1dy, dshape1dz;
00912 #endif
00913 
00914 public:
00915    LagrangeHexFiniteElement (int degree);
00916    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00917    virtual void CalcDShape(const IntegrationPoint &ip,
00918                            DenseMatrix &dshape) const;
00919    ~LagrangeHexFiniteElement ();
00920 };
00921 
00922 
00924 class RefinedLinear1DFiniteElement : public NodalFiniteElement
00925 {
00926 public:
00928    RefinedLinear1DFiniteElement();
00929 
00933    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00934 
00939    virtual void CalcDShape(const IntegrationPoint &ip,
00940                            DenseMatrix &dshape) const;
00941 };
00942 
00944 class RefinedLinear2DFiniteElement : public NodalFiniteElement
00945 {
00946 public:
00948    RefinedLinear2DFiniteElement();
00949 
00953    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00954 
00959    virtual void CalcDShape(const IntegrationPoint &ip,
00960                            DenseMatrix &dshape) const;
00961 };
00962 
00964 class RefinedLinear3DFiniteElement : public NodalFiniteElement
00965 {
00966 public:
00968    RefinedLinear3DFiniteElement();
00969 
00970    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00971 
00972    virtual void CalcDShape(const IntegrationPoint &ip,
00973                            DenseMatrix &dshape) const;
00974 };
00975 
00977 class RefinedBiLinear2DFiniteElement : public NodalFiniteElement
00978 {
00979 public:
00981    RefinedBiLinear2DFiniteElement();
00982 
00986    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
00987 
00992    virtual void CalcDShape(const IntegrationPoint &ip,
00993                            DenseMatrix &dshape) const;
00994 };
00995 
00997 class RefinedTriLinear3DFiniteElement : public NodalFiniteElement
00998 {
00999 public:
01001    RefinedTriLinear3DFiniteElement();
01002 
01006    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01007 
01012    virtual void CalcDShape(const IntegrationPoint &ip,
01013                            DenseMatrix &dshape) const;
01014 };
01015 
01016 
01017 class Nedelec1HexFiniteElement : public VectorFiniteElement
01018 {
01019 private:
01020    static const double tk[12][3];
01021 
01022 public:
01023    Nedelec1HexFiniteElement();
01024    virtual void CalcVShape(const IntegrationPoint &ip,
01025                            DenseMatrix &shape) const;
01026    virtual void CalcVShape(ElementTransformation &Trans,
01027                            DenseMatrix &shape) const
01028    { CalcVShape_ND(Trans, shape); };
01029    virtual void CalcCurlShape(const IntegrationPoint &ip,
01030                               DenseMatrix &curl_shape) const;
01031    virtual void GetLocalInterpolation (ElementTransformation &Trans,
01032                                        DenseMatrix &I) const;
01033    using FiniteElement::Project;
01034    virtual void Project (VectorCoefficient &vc,
01035                          ElementTransformation &Trans, Vector &dofs) const;
01036 };
01037 
01038 
01039 class Nedelec1TetFiniteElement : public VectorFiniteElement
01040 {
01041 private:
01042    static const double tk[6][3];
01043 
01044 public:
01045    Nedelec1TetFiniteElement();
01046    virtual void CalcVShape(const IntegrationPoint &ip,
01047                            DenseMatrix &shape) const;
01048    virtual void CalcVShape(ElementTransformation &Trans,
01049                            DenseMatrix &shape) const
01050    { CalcVShape_ND(Trans, shape); };
01051    virtual void CalcCurlShape(const IntegrationPoint &ip,
01052                               DenseMatrix &curl_shape) const;
01053    virtual void GetLocalInterpolation (ElementTransformation &Trans,
01054                                        DenseMatrix &I) const;
01055    using FiniteElement::Project;
01056    virtual void Project (VectorCoefficient &vc,
01057                          ElementTransformation &Trans, Vector &dofs) const;
01058 };
01059 
01060 
01061 class RT0HexFiniteElement : public VectorFiniteElement
01062 {
01063 private:
01064    static const double nk[6][3];
01065 
01066 public:
01067    RT0HexFiniteElement();
01068 
01069    virtual void CalcVShape(const IntegrationPoint &ip,
01070                            DenseMatrix &shape) const;
01071 
01072    virtual void CalcVShape(ElementTransformation &Trans,
01073                            DenseMatrix &shape) const
01074    { CalcVShape_RT(Trans, shape); };
01075 
01076    virtual void CalcDivShape(const IntegrationPoint &ip,
01077                              Vector &divshape) const;
01078 
01079    virtual void GetLocalInterpolation (ElementTransformation &Trans,
01080                                        DenseMatrix &I) const;
01081 
01082    using FiniteElement::Project;
01083 
01084    virtual void Project (VectorCoefficient &vc,
01085                          ElementTransformation &Trans, Vector &dofs) const;
01086 };
01087 
01088 
01089 class RT1HexFiniteElement : public VectorFiniteElement
01090 {
01091 private:
01092    static const double nk[36][3];
01093 
01094 public:
01095    RT1HexFiniteElement();
01096 
01097    virtual void CalcVShape(const IntegrationPoint &ip,
01098                            DenseMatrix &shape) const;
01099 
01100    virtual void CalcVShape(ElementTransformation &Trans,
01101                            DenseMatrix &shape) const
01102    { CalcVShape_RT(Trans, shape); };
01103 
01104    virtual void CalcDivShape(const IntegrationPoint &ip,
01105                              Vector &divshape) const;
01106 
01107    virtual void GetLocalInterpolation (ElementTransformation &Trans,
01108                                        DenseMatrix &I) const;
01109 
01110    using FiniteElement::Project;
01111 
01112    virtual void Project (VectorCoefficient &vc,
01113                          ElementTransformation &Trans, Vector &dofs) const;
01114 };
01115 
01116 
01117 class RT0TetFiniteElement : public VectorFiniteElement
01118 {
01119 private:
01120    static const double nk[4][3];
01121 
01122 public:
01123    RT0TetFiniteElement();
01124 
01125    virtual void CalcVShape(const IntegrationPoint &ip,
01126                            DenseMatrix &shape) const;
01127 
01128    virtual void CalcVShape(ElementTransformation &Trans,
01129                            DenseMatrix &shape) const
01130    { CalcVShape_RT(Trans, shape); };
01131 
01132    virtual void CalcDivShape(const IntegrationPoint &ip,
01133                              Vector &divshape) const;
01134 
01135    virtual void GetLocalInterpolation (ElementTransformation &Trans,
01136                                        DenseMatrix &I) const;
01137 
01138    using FiniteElement::Project;
01139 
01140    virtual void Project (VectorCoefficient &vc,
01141                          ElementTransformation &Trans, Vector &dofs) const;
01142 };
01143 
01144 
01145 class RotTriLinearHexFiniteElement : public NodalFiniteElement
01146 {
01147 public:
01148    RotTriLinearHexFiniteElement();
01149    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01150    virtual void CalcDShape(const IntegrationPoint &ip,
01151                            DenseMatrix &dshape) const;
01152 };
01153 
01154 
01155 class Poly_1D
01156 {
01157 public:
01158    class Basis
01159    {
01160    private:
01161       DenseMatrix A;
01162 #ifndef MFEM_USE_OPENMP
01163       mutable Vector a, b;
01164 #endif
01165 
01166    public:
01167       Basis(const int p, const double *nodes);
01168       void Eval(const double x, Vector &u) const;
01169       void Eval(const double x, Vector &u, Vector &d) const;
01170    };
01171 
01172 private:
01173    Array<double *> open_pts, closed_pts;
01174    Array<Basis *> open_basis, closed_basis;
01175 
01176    static void UniformPoints(const int p, double *x);
01177    static void GaussPoints(const int p, double *x);
01178    static void GaussLobattoPoints(const int p, double *x);
01179    static void ChebyshevPoints(const int p, double *x);
01180 
01181    static void CalcMono(const int p, const double x, double *u);
01182    static void CalcMono(const int p, const double x, double *u, double *d);
01183 
01184    static void CalcBernstein(const int p, const double x, double *u);
01185    static void CalcBernstein(const int p, const double x, double *u, double *d);
01186 
01187    static void CalcLegendre(const int p, const double x, double *u);
01188    static void CalcLegendre(const int p, const double x, double *u, double *d);
01189 
01190    static void CalcChebyshev(const int p, const double x, double *u);
01191    static void CalcChebyshev(const int p, const double x, double *u, double *d);
01192 
01193 public:
01194    Poly_1D() { }
01195 
01196    const double *OpenPoints(const int p);
01197    const double *ClosedPoints(const int p);
01198 
01199    Basis &OpenBasis(const int p);
01200    Basis &ClosedBasis(const int p);
01201 
01202    // Evaluate the values of a hierarchical 1D basis at point x
01203    // hierarchical = k-th basis function is degree k polynomial
01204    static void CalcBasis(const int p, const double x, double *u)
01205    // { CalcMono(p, x, u); }
01206    // Bernstein basis is not hierarchical --> does not work for triangles
01207    //  and tetrahedra
01208    // { CalcBernstein(p, x, u); }
01209    // { CalcLegendre(p, x, u); }
01210    { CalcChebyshev(p, x, u); }
01211 
01212    // Evaluate the values and derivatives of a hierarchical 1D basis at point x
01213    static void CalcBasis(const int p, const double x, double *u, double *d)
01214    // { CalcMono(p, x, u, d); }
01215    // { CalcBernstein(p, x, u, d); }
01216    // { CalcLegendre(p, x, u, d); }
01217    { CalcChebyshev(p, x, u, d); }
01218 
01219    // Evaluate a representation of a Delta function at point x
01220    static double CalcDelta(const int p, const double x){ return pow(x, (double) p); }
01221 
01222    ~Poly_1D();
01223 };
01224 
01225 extern Poly_1D poly1d;
01226 
01227 
01228 class H1_SegmentElement : public NodalFiniteElement
01229 {
01230 private:
01231    Poly_1D::Basis &basis1d;
01232 #ifndef MFEM_USE_OPENMP
01233    mutable Vector shape_x, dshape_x;
01234 #endif
01235 
01236 public:
01237    H1_SegmentElement(const int p);
01238    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01239    virtual void CalcDShape(const IntegrationPoint &ip,
01240                            DenseMatrix &dshape) const;
01241 };
01242 
01243 
01244 class H1_QuadrilateralElement : public NodalFiniteElement
01245 {
01246 private:
01247    Poly_1D::Basis &basis1d;
01248 #ifndef MFEM_USE_OPENMP
01249    mutable Vector shape_x, shape_y, dshape_x, dshape_y;
01250 #endif
01251    Array<int> dof_map;
01252 
01253 public:
01254    H1_QuadrilateralElement(const int p);
01255    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01256    virtual void CalcDShape(const IntegrationPoint &ip,
01257                            DenseMatrix &dshape) const;
01258 };
01259 
01260 
01261 class H1_HexahedronElement : public NodalFiniteElement
01262 {
01263 private:
01264    Poly_1D::Basis &basis1d;
01265 #ifndef MFEM_USE_OPENMP
01266    mutable Vector shape_x, shape_y, shape_z, dshape_x, dshape_y, dshape_z;
01267 #endif
01268    Array<int> dof_map;
01269 
01270 public:
01271    H1_HexahedronElement(const int p);
01272    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01273    virtual void CalcDShape(const IntegrationPoint &ip,
01274                            DenseMatrix &dshape) const;
01275 };
01276 
01277 
01278 class H1_TriangleElement : public NodalFiniteElement
01279 {
01280 private:
01281 #ifndef MFEM_USE_OPENMP
01282    mutable Vector shape_x, shape_y, shape_l, dshape_x, dshape_y, dshape_l, u;
01283    mutable DenseMatrix du;
01284 #endif
01285    DenseMatrix T;
01286 
01287 public:
01288    H1_TriangleElement(const int p);
01289    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01290    virtual void CalcDShape(const IntegrationPoint &ip,
01291                            DenseMatrix &dshape) const;
01292 };
01293 
01294 
01295 class H1_TetrahedronElement : public NodalFiniteElement
01296 {
01297 private:
01298 #ifndef MFEM_USE_OPENMP
01299    mutable Vector shape_x, shape_y, shape_z, shape_l;
01300    mutable Vector dshape_x, dshape_y, dshape_z, dshape_l, u;
01301    mutable DenseMatrix du;
01302 #endif
01303    DenseMatrix T;
01304 
01305 public:
01306    H1_TetrahedronElement(const int p);
01307    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01308    virtual void CalcDShape(const IntegrationPoint &ip,
01309                            DenseMatrix &dshape) const;
01310 };
01311 
01312 
01313 class L2_SegmentElement : public NodalFiniteElement
01314 {
01315 private:
01316    Poly_1D::Basis &basis1d;
01317 #ifndef MFEM_USE_OPENMP
01318    mutable Vector shape_x, dshape_x;
01319 #endif
01320 
01321 public:
01322    L2_SegmentElement(const int p);
01323    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01324    virtual void CalcDShape(const IntegrationPoint &ip,
01325                            DenseMatrix &dshape) const;
01326    virtual void ProjectDelta(int vertex, Vector &dofs) const;
01327 };
01328 
01329 
01330 class L2_QuadrilateralElement : public NodalFiniteElement
01331 {
01332 private:
01333    Poly_1D::Basis &basis1d;
01334 #ifndef MFEM_USE_OPENMP
01335    mutable Vector shape_x, shape_y, dshape_x, dshape_y;
01336 #endif
01337 
01338 public:
01339    L2_QuadrilateralElement(const int p);
01340    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01341    virtual void CalcDShape(const IntegrationPoint &ip,
01342                            DenseMatrix &dshape) const;
01343    virtual void ProjectDelta(int vertex, Vector &dofs) const;
01344 };
01345 
01346 
01347 class L2_HexahedronElement : public NodalFiniteElement
01348 {
01349 private:
01350    Poly_1D::Basis &basis1d;
01351 #ifndef MFEM_USE_OPENMP
01352    mutable Vector shape_x, shape_y, shape_z, dshape_x, dshape_y, dshape_z;
01353 #endif
01354 
01355 public:
01356    L2_HexahedronElement(const int p);
01357    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01358    virtual void CalcDShape(const IntegrationPoint &ip,
01359                            DenseMatrix &dshape) const;
01360    virtual void ProjectDelta(int vertex, Vector &dofs) const;
01361 };
01362 
01363 
01364 class L2_TriangleElement : public NodalFiniteElement
01365 {
01366 private:
01367 #ifndef MFEM_USE_OPENMP
01368    mutable Vector shape_x, shape_y, shape_l, dshape_x, dshape_y, dshape_l, u;
01369    mutable DenseMatrix du;
01370 #endif
01371    DenseMatrix T;
01372 
01373 public:
01374    L2_TriangleElement(const int p);
01375    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01376    virtual void CalcDShape(const IntegrationPoint &ip,
01377                            DenseMatrix &dshape) const;
01378 };
01379 
01380 
01381 class L2_TetrahedronElement : public NodalFiniteElement
01382 {
01383 private:
01384 #ifndef MFEM_USE_OPENMP
01385    mutable Vector shape_x, shape_y, shape_z, shape_l;
01386    mutable Vector dshape_x, dshape_y, dshape_z, dshape_l, u;
01387    mutable DenseMatrix du;
01388 #endif
01389    DenseMatrix T;
01390 
01391 public:
01392    L2_TetrahedronElement(const int p);
01393    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01394    virtual void CalcDShape(const IntegrationPoint &ip,
01395                            DenseMatrix &dshape) const;
01396 };
01397 
01398 
01399 class RT_QuadrilateralElement : public VectorFiniteElement
01400 {
01401 private:
01402    static const double nk[8];
01403 
01404    Poly_1D::Basis &cbasis1d, &obasis1d;
01405 #ifndef MFEM_USE_OPENMP
01406    mutable Vector shape_cx, shape_ox, shape_cy, shape_oy;
01407    mutable Vector dshape_cx, dshape_cy;
01408 #endif
01409    Array<int> dof_map, dof2nk;
01410 
01411 public:
01412    RT_QuadrilateralElement(const int p);
01413    virtual void CalcVShape(const IntegrationPoint &ip,
01414                            DenseMatrix &shape) const;
01415    virtual void CalcVShape(ElementTransformation &Trans,
01416                            DenseMatrix &shape) const
01417    { CalcVShape_RT(Trans, shape); }
01418    virtual void CalcDivShape(const IntegrationPoint &ip,
01419                              Vector &divshape) const;
01420    virtual void GetLocalInterpolation(ElementTransformation &Trans,
01421                                       DenseMatrix &I) const
01422    { LocalInterpolation_RT(nk, dof2nk, Trans, I); }
01423    using FiniteElement::Project;
01424    virtual void Project(VectorCoefficient &vc,
01425                         ElementTransformation &Trans, Vector &dofs) const
01426    { Project_RT(nk, dof2nk, vc, Trans, dofs); }
01427    virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
01428                         DenseMatrix &I) const
01429    { Project_RT(nk, dof2nk, fe, Trans, I); }
01430    virtual void ProjectGrad(const FiniteElement &fe,
01431                             ElementTransformation &Trans,
01432                             DenseMatrix &grad) const
01433    { ProjectGrad_RT(nk, dof2nk, fe, Trans, grad); }
01434 };
01435 
01436 
01437 class RT_HexahedronElement : public VectorFiniteElement
01438 {
01439    static const double nk[18];
01440 
01441    Poly_1D::Basis &cbasis1d, &obasis1d;
01442 #ifndef MFEM_USE_OPENMP
01443    mutable Vector shape_cx, shape_ox, shape_cy, shape_oy, shape_cz, shape_oz;
01444    mutable Vector dshape_cx, dshape_cy, dshape_cz;
01445 #endif
01446    Array<int> dof_map, dof2nk;
01447 
01448 public:
01449    RT_HexahedronElement(const int p);
01450    virtual void CalcVShape(const IntegrationPoint &ip,
01451                            DenseMatrix &shape) const;
01452    virtual void CalcVShape(ElementTransformation &Trans,
01453                            DenseMatrix &shape) const
01454    { CalcVShape_RT(Trans, shape); }
01455    virtual void CalcDivShape(const IntegrationPoint &ip,
01456                              Vector &divshape) const;
01457    virtual void GetLocalInterpolation(ElementTransformation &Trans,
01458                                       DenseMatrix &I) const
01459    { LocalInterpolation_RT(nk, dof2nk, Trans, I); }
01460    using FiniteElement::Project;
01461    virtual void Project(VectorCoefficient &vc,
01462                         ElementTransformation &Trans, Vector &dofs) const
01463    { Project_RT(nk, dof2nk, vc, Trans, dofs); }
01464    virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
01465                         DenseMatrix &I) const
01466    { Project_RT(nk, dof2nk, fe, Trans, I); }
01467    virtual void ProjectCurl(const FiniteElement &fe,
01468                             ElementTransformation &Trans,
01469                             DenseMatrix &curl) const
01470    { ProjectCurl_RT(nk, dof2nk, fe, Trans, curl); }
01471 };
01472 
01473 
01474 class RT_TriangleElement : public VectorFiniteElement
01475 {
01476    static const double nk[6], c;
01477 
01478 #ifndef MFEM_USE_OPENMP
01479    mutable Vector shape_x, shape_y, shape_l;
01480    mutable Vector dshape_x, dshape_y, dshape_l;
01481    mutable DenseMatrix u;
01482    mutable Vector divu;
01483 #endif
01484    Array<int> dof2nk;
01485    DenseMatrix T;
01486 
01487 public:
01488    RT_TriangleElement(const int p);
01489    virtual void CalcVShape(const IntegrationPoint &ip,
01490                            DenseMatrix &shape) const;
01491    virtual void CalcVShape(ElementTransformation &Trans,
01492                            DenseMatrix &shape) const
01493    { CalcVShape_RT(Trans, shape); }
01494    virtual void CalcDivShape(const IntegrationPoint &ip,
01495                              Vector &divshape) const;
01496    virtual void GetLocalInterpolation(ElementTransformation &Trans,
01497                                       DenseMatrix &I) const
01498    { LocalInterpolation_RT(nk, dof2nk, Trans, I); }
01499    using FiniteElement::Project;
01500    virtual void Project(VectorCoefficient &vc,
01501                         ElementTransformation &Trans, Vector &dofs) const
01502    { Project_RT(nk, dof2nk, vc, Trans, dofs); }
01503    virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
01504                         DenseMatrix &I) const
01505    { Project_RT(nk, dof2nk, fe, Trans, I); }
01506    virtual void ProjectGrad(const FiniteElement &fe,
01507                             ElementTransformation &Trans,
01508                             DenseMatrix &grad) const
01509    { ProjectGrad_RT(nk, dof2nk, fe, Trans, grad); }
01510 };
01511 
01512 
01513 class RT_TetrahedronElement : public VectorFiniteElement
01514 {
01515    static const double nk[12], c;
01516 
01517 #ifndef MFEM_USE_OPENMP
01518    mutable Vector shape_x, shape_y, shape_z, shape_l;
01519    mutable Vector dshape_x, dshape_y, dshape_z, dshape_l;
01520    mutable DenseMatrix u;
01521    mutable Vector divu;
01522 #endif
01523    Array<int> dof2nk;
01524    DenseMatrix T;
01525 
01526 public:
01527    RT_TetrahedronElement(const int p);
01528    virtual void CalcVShape(const IntegrationPoint &ip,
01529                            DenseMatrix &shape) const;
01530    virtual void CalcVShape(ElementTransformation &Trans,
01531                            DenseMatrix &shape) const
01532    { CalcVShape_RT(Trans, shape); }
01533    virtual void CalcDivShape(const IntegrationPoint &ip,
01534                              Vector &divshape) const;
01535    virtual void GetLocalInterpolation(ElementTransformation &Trans,
01536                                       DenseMatrix &I) const
01537    { LocalInterpolation_RT(nk, dof2nk, Trans, I); }
01538    using FiniteElement::Project;
01539    virtual void Project(VectorCoefficient &vc,
01540                         ElementTransformation &Trans, Vector &dofs) const
01541    { Project_RT(nk, dof2nk, vc, Trans, dofs); }
01542    virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
01543                         DenseMatrix &I) const
01544    { Project_RT(nk, dof2nk, fe, Trans, I); }
01545    virtual void ProjectCurl(const FiniteElement &fe,
01546                             ElementTransformation &Trans,
01547                             DenseMatrix &curl) const
01548    { ProjectCurl_RT(nk, dof2nk, fe, Trans, curl); }
01549 };
01550 
01551 
01552 class ND_HexahedronElement : public VectorFiniteElement
01553 {
01554    static const double tk[18];
01555 
01556    Poly_1D::Basis &cbasis1d, &obasis1d;
01557 #ifndef MFEM_USE_OPENMP
01558    mutable Vector shape_cx, shape_ox, shape_cy, shape_oy, shape_cz, shape_oz;
01559    mutable Vector dshape_cx, dshape_cy, dshape_cz;
01560 #endif
01561    Array<int> dof_map, dof2tk;
01562 
01563 public:
01564    ND_HexahedronElement(const int p);
01565    virtual void CalcVShape(const IntegrationPoint &ip,
01566                            DenseMatrix &shape) const;
01567    virtual void CalcVShape(ElementTransformation &Trans,
01568                            DenseMatrix &shape) const
01569    { CalcVShape_ND(Trans, shape); }
01570    virtual void CalcCurlShape(const IntegrationPoint &ip,
01571                               DenseMatrix &curl_shape) const;
01572    virtual void GetLocalInterpolation(ElementTransformation &Trans,
01573                                       DenseMatrix &I) const
01574    { LocalInterpolation_ND(tk, dof2tk, Trans, I); }
01575    using FiniteElement::Project;
01576    virtual void Project(VectorCoefficient &vc,
01577                         ElementTransformation &Trans, Vector &dofs) const
01578    { Project_ND(tk, dof2tk, vc, Trans, dofs); }
01579    virtual void Project(const FiniteElement &fe,
01580                         ElementTransformation &Trans,
01581                         DenseMatrix &I) const
01582    { Project_ND(tk, dof2tk, fe, Trans, I); }
01583    virtual void ProjectGrad(const FiniteElement &fe,
01584                             ElementTransformation &Trans,
01585                             DenseMatrix &grad) const
01586    { ProjectGrad_ND(tk, dof2tk, fe, Trans, grad); }
01587 };
01588 
01589 
01590 class ND_QuadrilateralElement : public VectorFiniteElement
01591 {
01592    static const double tk[8];
01593 
01594    Poly_1D::Basis &cbasis1d, &obasis1d;
01595 #ifndef MFEM_USE_OPENMP
01596    mutable Vector shape_cx, shape_ox, shape_cy, shape_oy;
01597    mutable Vector dshape_cx, dshape_cy;
01598 #endif
01599    Array<int> dof_map, dof2tk;
01600 
01601 public:
01602    ND_QuadrilateralElement(const int p);
01603    virtual void CalcVShape(const IntegrationPoint &ip,
01604                            DenseMatrix &shape) const;
01605    virtual void CalcVShape(ElementTransformation &Trans,
01606                            DenseMatrix &shape) const
01607    { CalcVShape_ND(Trans, shape); }
01608    virtual void CalcCurlShape(const IntegrationPoint &ip,
01609                               DenseMatrix &curl_shape) const;
01610    virtual void GetLocalInterpolation(ElementTransformation &Trans,
01611                                       DenseMatrix &I) const
01612    { LocalInterpolation_ND(tk, dof2tk, Trans, I); }
01613    using FiniteElement::Project;
01614    virtual void Project(VectorCoefficient &vc,
01615                         ElementTransformation &Trans, Vector &dofs) const
01616    { Project_ND(tk, dof2tk, vc, Trans, dofs); }
01617    virtual void Project(const FiniteElement &fe,
01618                         ElementTransformation &Trans,
01619                         DenseMatrix &I) const
01620    { Project_ND(tk, dof2tk, fe, Trans, I); }
01621    virtual void ProjectGrad(const FiniteElement &fe,
01622                             ElementTransformation &Trans,
01623                             DenseMatrix &grad) const
01624    { ProjectGrad_ND(tk, dof2tk, fe, Trans, grad); }
01625 };
01626 
01627 
01628 class ND_TetrahedronElement : public VectorFiniteElement
01629 {
01630    static const double tk[18], c;
01631 
01632 #ifndef MFEM_USE_OPENMP
01633    mutable Vector shape_x, shape_y, shape_z, shape_l;
01634    mutable Vector dshape_x, dshape_y, dshape_z, dshape_l;
01635    mutable DenseMatrix u;
01636 #endif
01637    Array<int> dof2tk;
01638    DenseMatrix T;
01639 
01640 public:
01641    ND_TetrahedronElement(const int p);
01642    virtual void CalcVShape(const IntegrationPoint &ip,
01643                            DenseMatrix &shape) const;
01644    virtual void CalcVShape(ElementTransformation &Trans,
01645                            DenseMatrix &shape) const
01646    { CalcVShape_ND(Trans, shape); }
01647    virtual void CalcCurlShape(const IntegrationPoint &ip,
01648                               DenseMatrix &curl_shape) const;
01649    virtual void GetLocalInterpolation(ElementTransformation &Trans,
01650                                       DenseMatrix &I) const
01651    { LocalInterpolation_ND(tk, dof2tk, Trans, I); }
01652    using FiniteElement::Project;
01653    virtual void Project(VectorCoefficient &vc,
01654                         ElementTransformation &Trans, Vector &dofs) const
01655    { Project_ND(tk, dof2tk, vc, Trans, dofs); }
01656    virtual void Project(const FiniteElement &fe,
01657                         ElementTransformation &Trans,
01658                         DenseMatrix &I) const
01659    { Project_ND(tk, dof2tk, fe, Trans, I); }
01660    virtual void ProjectGrad(const FiniteElement &fe,
01661                             ElementTransformation &Trans,
01662                             DenseMatrix &grad) const
01663    { ProjectGrad_ND(tk, dof2tk, fe, Trans, grad); }
01664 };
01665 
01666 class ND_TriangleElement : public VectorFiniteElement
01667 {
01668    static const double tk[8], c;
01669 
01670 #ifndef MFEM_USE_OPENMP
01671    mutable Vector shape_x, shape_y, shape_l;
01672    mutable Vector dshape_x, dshape_y, dshape_l;
01673    mutable DenseMatrix u;
01674 #endif
01675    Array<int> dof2tk;
01676    DenseMatrix T;
01677 
01678 public:
01679    ND_TriangleElement(const int p);
01680    virtual void CalcVShape(const IntegrationPoint &ip,
01681                            DenseMatrix &shape) const;
01682    virtual void CalcVShape(ElementTransformation &Trans,
01683                            DenseMatrix &shape) const
01684    { CalcVShape_ND(Trans, shape); }
01685    virtual void CalcCurlShape(const IntegrationPoint &ip,
01686                               DenseMatrix &curl_shape) const;
01687    virtual void GetLocalInterpolation(ElementTransformation &Trans,
01688                                       DenseMatrix &I) const
01689    { LocalInterpolation_ND(tk, dof2tk, Trans, I); }
01690    using FiniteElement::Project;
01691    virtual void Project(VectorCoefficient &vc,
01692                         ElementTransformation &Trans, Vector &dofs) const
01693    { Project_ND(tk, dof2tk, vc, Trans, dofs); }
01694    virtual void Project(const FiniteElement &fe,
01695                         ElementTransformation &Trans,
01696                         DenseMatrix &I) const
01697    { Project_ND(tk, dof2tk, fe, Trans, I); }
01698    virtual void ProjectGrad(const FiniteElement &fe,
01699                             ElementTransformation &Trans,
01700                             DenseMatrix &grad) const
01701    { ProjectGrad_ND(tk, dof2tk, fe, Trans, grad); }
01702 };
01703 
01704 
01705 class NURBSFiniteElement : public FiniteElement
01706 {
01707 protected:
01708    mutable Array <KnotVector*> kv;
01709    mutable int *ijk, patch, elem;
01710    mutable Vector weights;
01711 
01712 public:
01713    NURBSFiniteElement(int D, int G, int Do, int O, int F)
01714       : FiniteElement(D, G, Do, O, F)
01715    {
01716       ijk = NULL;
01717       patch = elem = -1;
01718       kv.SetSize(Dim);
01719       weights.SetSize(Dof);
01720       weights = 1.0;
01721    }
01722 
01723    void                 Reset      ()         const { patch = elem = -1; }
01724    void                 SetIJK     (int *IJK) const { ijk = IJK; }
01725    int                  GetPatch   ()         const { return patch; }
01726    void                 SetPatch   (int p)    const { patch = p; }
01727    int                  GetElement ()         const { return elem; }
01728    void                 SetElement (int e)    const { elem = e; }
01729    Array <KnotVector*> &KnotVectors()         const { return kv; }
01730    Vector              &Weights    ()         const { return weights; }
01731 };
01732 
01733 class NURBS1DFiniteElement : public NURBSFiniteElement
01734 {
01735 protected:
01736    mutable Vector shape_x;
01737 
01738 public:
01739    NURBS1DFiniteElement(int p)
01740       : NURBSFiniteElement(1, Geometry::SEGMENT, p + 1, p, FunctionSpace::Qk),
01741         shape_x(p + 1) { }
01742 
01743    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01744    virtual void CalcDShape(const IntegrationPoint &ip,
01745                            DenseMatrix &dshape) const;
01746 };
01747 
01748 class NURBS2DFiniteElement : public NURBSFiniteElement
01749 {
01750 protected:
01751    mutable Vector u, shape_x, shape_y, dshape_x, dshape_y;
01752 
01753 public:
01754    NURBS2DFiniteElement(int p)
01755       : NURBSFiniteElement(2, Geometry::SQUARE, (p + 1)*(p + 1), p,
01756                            FunctionSpace::Qk), u(Dof),
01757         shape_x(p + 1), shape_y(p + 1), dshape_x(p + 1), dshape_y(p + 1) { }
01758 
01759    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01760    virtual void CalcDShape(const IntegrationPoint &ip,
01761                            DenseMatrix &dshape) const;
01762 };
01763 
01764 class NURBS3DFiniteElement : public NURBSFiniteElement
01765 {
01766 protected:
01767    mutable Vector u, shape_x, shape_y, shape_z, dshape_x, dshape_y, dshape_z;
01768 
01769 public:
01770    NURBS3DFiniteElement(int p)
01771       : NURBSFiniteElement(3, Geometry::CUBE, (p + 1)*(p + 1)*(p + 1), p,
01772                            FunctionSpace::Qk), u(Dof),
01773         shape_x(p + 1), shape_y(p + 1), shape_z(p + 1),
01774         dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1) { }
01775 
01776    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
01777    virtual void CalcDShape(const IntegrationPoint &ip,
01778                            DenseMatrix &dshape) const;
01779 };
01780 
01781 #endif