/usr/include/dune/localfunctions/dualmortarbasis/dualq1.hh is in libdune-localfunctions-dev 2.4.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 | // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
#define DUNE_DUAL_Q1_LOCALFINITEELEMENT_HH
#include <array>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include <dune/localfunctions/lagrange/q1/q1localbasis.hh>
#include "dualq1/dualq1localbasis.hh"
#include "dualq1/dualq1localcoefficients.hh"
#include "dualq1/dualq1localinterpolation.hh"
namespace Dune
{
/**
* \brief The local dual Q1 finite element on cubes
*
* Note that if the dual functions are chosen to be dual on the faces,
* the integrated product of a Lagrange \f$\lambda_p\f$ and dual
* function \f$\theta_q\f$ over faces not containing \f$q\f$ does in
* general not vanish.
*
* \ingroup DualMortar
*
* \tparam D Domain data type
* \tparam R Range data type
* \tparam dim Dimension of the hypercube
* \tparam faceDual If set, the basis functions are bi-orthogonal only on faces containing the corresponding vertex.
*/
template<class D, class R, int dim, bool faceDual=false>
class DualQ1LocalFiniteElement
{
public:
/** \todo Please doc me !
*/
typedef LocalFiniteElementTraits<DualQ1LocalBasis<D,R,dim>,DualQ1LocalCoefficients<dim>,
DualQ1LocalInterpolation<dim,DualQ1LocalBasis<D,R,dim> > > Traits;
/** \todo Please doc me !
*/
DualQ1LocalFiniteElement ()
{
gt.makeCube(dim);
if (faceDual)
setupFaceDualCoefficients();
else
setupDualCoefficients();
}
/** \todo Please doc me !
*/
const typename Traits::LocalBasisType& localBasis () const
{
return basis;
}
/** \todo Please doc me !
*/
const typename Traits::LocalCoefficientsType& localCoefficients () const
{
return coefficients;
}
/** \todo Please doc me !
*/
const typename Traits::LocalInterpolationType& localInterpolation () const
{
return interpolation;
}
/** \brief Number of shape functions in this finite element */
unsigned int size () const
{
return basis.size();
}
/** \todo Please doc me !
*/
GeometryType type () const
{
return gt;
}
DualQ1LocalFiniteElement* clone () const
{
return new DualQ1LocalFiniteElement(*this);
}
private:
/** \brief Setup the coefficients for the face dual basis functions. */
void setupFaceDualCoefficients();
/** \brief Setup the coefficients for the dual basis functions. */
void setupDualCoefficients();
DualQ1LocalBasis<D,R,dim> basis;
DualQ1LocalCoefficients<dim> coefficients;
DualQ1LocalInterpolation<dim,DualQ1LocalBasis<D,R,dim> > interpolation;
GeometryType gt;
};
template<class D, class R, int dim, bool faceDual>
void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupDualCoefficients()
{
const int size = 1 <<dim;
std::array<Dune::FieldVector<R, size>, size> coeffs;
// dual basis functions are linear combinations of Lagrange elements
// compute these coefficients here because the basis and the local interpolation needs them
const Dune::QuadratureRule<D,dim>& quad = Dune::QuadratureRules<D,dim>::rule(gt, 2*dim);
// assemble mass matrix on the reference element
Dune::FieldMatrix<R, size, size> massMat;
massMat = 0;
// and the integrals of the lagrange shape functions
std::vector<Dune::FieldVector<R,1> > integral(size);
for (int i=0; i<size; i++)
integral[i] = 0;
Dune::Q1LocalBasis<D,R,dim> q1Basis;
for(size_t pt=0; pt<quad.size(); pt++) {
const Dune::FieldVector<D ,dim>& pos = quad[pt].position();
std::vector<Dune::FieldVector<R,1> > q1Values(size);
q1Basis.evaluateFunction(pos,q1Values);
D weight = quad[pt].weight();
for (int k=0; k<size; k++) {
integral[k] += q1Values[k]*weight;
for (int l=0; l<=k; l++)
massMat[k][l] += weight*(q1Values[k]*q1Values[l]);
}
}
// make matrix symmetric
for (int i=0; i<size-1; i++)
for (int j=i+1; j<size; j++)
massMat[i][j] = massMat[j][i];
//solve for the coefficients
for (int i=0; i<size; i++) {
Dune::FieldVector<R, size> rhs(0);
rhs[i] = integral[i];
coeffs[i] = 0;
massMat.solve(coeffs[i] ,rhs);
}
basis.setCoefficients(coeffs);
interpolation.setCoefficients(coeffs);
}
template<class D, class R, int dim, bool faceDual>
void DualQ1LocalFiniteElement<D,R,dim,faceDual>::setupFaceDualCoefficients()
{
const int size = 1 <<dim;
std::array<Dune::FieldVector<R, size>, size> coeffs;
// dual basis functions are linear combinations of Lagrange elements
Dune::Q1LocalBasis<D,R,dim> q1Basis;
typedef Dune::ReferenceElement<D,dim> RefElement;
const RefElement& refElement = Dune::ReferenceElements<D,dim>::general(gt);
// loop over faces
for (size_t i=0; i<refElement.size(1);i++) {
const Dune::QuadratureRule<D,dim-1>& quad = Dune::QuadratureRules<D,dim-1>::rule(refElement.type(i,1),2*dim);
// for each face assemble the mass matrix over the face of all
// non-vanishing basis functions,
// for cubes refElement.size(i,1,dim)=size/2
Dune::FieldMatrix<R, size/2, size/2> massMat;
massMat = 0;
// get geometry
const typename RefElement::template Codim<1>::Geometry& geometry = refElement.template geometry<1>(i);
// and the integrals of the lagrange shape functions
std::vector<Dune::FieldVector<R,1> > integral(size/2);
for (int k=0; k<size/2; k++)
integral[k] = 0;
for(size_t pt=0; pt<quad.size(); pt++) {
const Dune::FieldVector<D,dim-1>& pos = quad[pt].position();
const Dune::FieldVector<D,dim> elementPos = geometry.global(pos);
std::vector<Dune::FieldVector<R,1> > q1Values(size);
q1Basis.evaluateFunction(elementPos,q1Values);
D weight = quad[pt].weight();
for (int k=0; k<refElement.size(i,1,dim); k++) {
int row = refElement.subEntity(i,1,k,dim);
integral[k] += q1Values[row]*weight;
for (int l=0; l<refElement.size(i,1,dim); l++) {
int col = refElement.subEntity(i,1,l,dim);
massMat[k][l] += weight*(q1Values[row]*q1Values[col]);
}
}
}
// solve for the coefficients
// note that we possibly overwrite coefficients for neighbouring faces
// this is okay since the coefficients are symmetric
for (int l=0; l<refElement.size(i,1,dim); l++) {
int row = refElement.subEntity(i,1,l,dim);
Dune::FieldVector<R, size/2> rhs(0);
rhs[l] = integral[l];
Dune::FieldVector<R, size/2> x(0);
massMat.solve(x ,rhs);
for (int k=0; k<refElement.size(i,1,dim); k++) {
int col = refElement.subEntity(i,1,k,dim);
coeffs[row][col]=x[k];
}
}
}
basis.setCoefficients(coeffs);
interpolation.setCoefficients(coeffs);
}
}
#endif
|