dune-localfunctions  2.2.0
whitney/edges0.5/basis.hh
Go to the documentation of this file.
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set ts=8 sw=2 et sts=2:
00003 
00004 #ifndef DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH
00005 #define DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH
00006 
00007 #include <cstddef>
00008 #include <vector>
00009 
00010 #include <dune/common/fmatrix.hh>
00011 #include <dune/common/fvector.hh>
00012 // #include <dune/common/geometrytype.hh>
00013 
00014 // #include <dune/grid/common/referenceelements.hh>
00015 
00016 // #include <dune/localfunctions/common/localkey.hh>
00017 #include <dune/localfunctions/common/localtoglobaladaptors.hh>
00018 #include <dune/localfunctions/lagrange/p1/p1localbasis.hh>
00019 #include <dune/localfunctions/whitney/edges0.5/common.hh>
00020 
00021 namespace Dune {
00022 
00024   //
00025   //  Basis
00026   //
00027 
00029 
00037   template<class Geometry, class RF>
00038   class EdgeS0_5Basis :
00039     private EdgeS0_5Common<Geometry::mydimension, typename Geometry::ctype>
00040   {
00041   public:
00043     struct Traits {
00044       typedef typename Geometry::ctype DomainField;
00045       static const std::size_t dimDomainLocal = Geometry::mydimension;
00046       static const std::size_t dimDomainGlobal = Geometry::coorddimension;
00047       typedef FieldVector<DomainField, dimDomainLocal> DomainLocal;
00048       typedef FieldVector<DomainField, dimDomainGlobal> DomainGlobal;
00049 
00050       typedef RF RangeField;
00051       static const std::size_t dimRange = dimDomainLocal;
00052       typedef FieldVector<RangeField, dimRange> Range;
00053 
00054       typedef FieldMatrix<RangeField, dimRange, dimDomainGlobal> Jacobian;
00055 
00056       static const std::size_t diffOrder = 1;
00057     };
00058 
00059   private:
00060     typedef Dune::P1LocalBasis<typename Traits::DomainField,
00061                                typename Traits::RangeField,
00062                                Traits::dimDomainLocal
00063                                > P1LocalBasis;
00064     typedef ScalarLocalToGlobalBasisAdaptor<P1LocalBasis, Geometry> P1Basis;
00065 
00066     static const P1LocalBasis& p1LocalBasis;
00067     static const std::size_t dim = Traits::dimDomainLocal;
00068 
00069     typedef EdgeS0_5Common<dim, typename Geometry::ctype> Base;
00070     using Base::refelem;
00071     using Base::s;
00072 
00073     // global values of the Jacobians (gradients) of the p1 basis
00074     std::vector<typename P1Basis::Traits::Jacobian> p1j;
00075     // edge sizes and orientations
00076     std::vector<typename Traits::DomainField> edgel;
00077 
00078   public:
00080 
00086     template<typename VertexOrder>
00087     EdgeS0_5Basis(const Geometry& geo, const VertexOrder& vertexOrder) :
00088       p1j(s, typename P1Basis::Traits::Jacobian(0)), edgel(s)
00089     {
00090       // use some arbitrary position to evaluate jacobians, they are constant
00091       static const typename Traits::DomainLocal xl(0);
00092 
00093       // precompute Jacobian (gradients) of the p1 element
00094       P1Basis(p1LocalBasis, geo).evaluateJacobian(xl, p1j);
00095 
00096       // calculate edge sizes and orientations
00097       for(std::size_t i = 0; i < s; ++i) {
00098         edgel[i] = (geo.corner(refelem.subEntity(i,dim-1,0,dim))-
00099                     geo.corner(refelem.subEntity(i,dim-1,1,dim))
00100                     ).two_norm();
00101         const typename VertexOrder::iterator& edgeVertexOrder =
00102           vertexOrder.begin(dim-1, i);
00103         if(edgeVertexOrder[0] > edgeVertexOrder[1])
00104           edgel[i] *= -1;
00105       }
00106     }
00107 
00109     std::size_t size () const { return s; }
00110 
00112     void evaluateFunction(const typename Traits::DomainLocal& xl,
00113                           std::vector<typename Traits::Range>& out) const
00114     {
00115       out.assign(s, typename Traits::Range(0));
00116 
00117       // compute p1 values -- use the local basis directly for that, local and
00118       // global values are identical for scalars
00119       std::vector<typename P1LocalBasis::Traits::RangeType> p1v;
00120       p1LocalBasis.evaluateFunction(xl, p1v);
00121 
00122       for(std::size_t i = 0; i < s; i++) {
00123         const std::size_t i0 = refelem.subEntity(i,dim-1,0,dim);
00124         const std::size_t i1 = refelem.subEntity(i,dim-1,1,dim);
00125         out[i].axpy( p1v[i0], p1j[i1][0]);
00126         out[i].axpy(-p1v[i1], p1j[i0][0]);
00127         out[i] *= edgel[i];
00128       }
00129     }
00130 
00132     void evaluateJacobian(const typename Traits::DomainLocal&,
00133                           std::vector<typename Traits::Jacobian>& out) const
00134     {
00135       out.resize(s);
00136 
00137       for(std::size_t i = 0; i < s; i++) {
00138         const std::size_t i0 = refelem.subEntity(i,dim-1,0,dim);
00139         const std::size_t i1 = refelem.subEntity(i,dim-1,1,dim);
00140         for(std::size_t j = 0; j < dim; j++)
00141           for(std::size_t k = 0; k < dim; k++)
00142             out[i][j][k] = edgel[i] *
00143               (p1j[i0][0][k]*p1j[i1][0][j]-p1j[i1][0][k]*p1j[i0][0][j]);
00144       }
00145     }
00146 
00148     std::size_t order () const { return 1; }
00149   };
00150 
00151   template<class Geometry, class RF>
00152   const typename EdgeS0_5Basis<Geometry, RF>::P1LocalBasis&
00153   EdgeS0_5Basis<Geometry, RF>::p1LocalBasis = P1LocalBasis();
00154 
00155 } // namespace Dune
00156 
00157 #endif // DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH