bump local functions

parent 7ae09a4b
Pipeline #9307 failed with stage
in 0 seconds
add_subdirectory(brezzidouglasmarini)
add_subdirectory(common)
add_subdirectory(bump)
add_subdirectory(dualmortarbasis)
add_subdirectory(hierarchical)
add_subdirectory(lagrange)
......
add_subdirectory(test)
install(FILES
bumplocalbasis.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/localfunctions/bump)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_BUMPLOCALFINITEELEMENT_HH
#define DUNE_BUMPLOCALFINITEELEMENT_HH
#include <dune/geometry/type.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include "bumplocalbasis.hh"
#include "bumplocalcoefficients.hh"
#include "bumplocalinterpolation.hh"
namespace Dune
{
/** \brief The local p0 finite element on all types of reference elements
\tparam D Domain data type
\tparam R Range data type
\tparam d Dimension of the reference element
*/
template<class D, class R, int d>
class BumpLocalFiniteElement
{
public:
/** \todo Please doc me !
*/
typedef LocalFiniteElementTraits<BumpLocalBasis<D,R,d>, BumpLocalCoefficients<d>,
BumpLocalInterpolation<BumpLocalBasis<D,R,d> > > Traits;
/** \todo Please doc me !
*/
BumpLocalFiniteElement (const GeometryType& type)
: basis(type)
, coefficients(type)
, interpolation(type)
, gt(type)
{}
/** \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 The number of shape functions -- here: 1 */
unsigned int size () const
{
return basis.size();
}
/** \todo Please doc me !
*/
GeometryType type () const
{
return gt;
}
BumpLocalFiniteElement* clone () const
{
return new BumpLocalFiniteElement(*this);
}
private:
BumpLocalBasis<D,R,d> basis;
BumpLocalCoefficients<d> coefficients;
BumpLocalInterpolation<BumpLocalBasis<D,R,d> > interpolation;
GeometryType gt;
};
}
#endif
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_BUMP_LOCALBASIS_HH
#define DUNE_BUMP_LOCALBASIS_HH
#include <array>
#include <numeric>
#include <iostream>
#include <dune/common/deprecated.hh>
#include <dune/common/fmatrix.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/localfunctions/common/localbasis.hh>
namespace Dune
{
/**@ingroup LocalBasisImplementation
\brief Linear Lagrange shape functions on the simplex.
Defines the linear shape functions on the simplex.
\tparam D Type to represent the field in the domain.
\tparam R Type to represent the field in the range.
\tparam dim The dimension of the simplex
\nosubgrouping
*/
template<class D, class R, int dim>
class BumpLocalBasis;
template<class D, class R, int dim>
class BumpLocalBasis
{
public:
BumpLocalBasis(const GeometryType& type) : ref_geo(ReferenceElements<double, dim>::general(type))
{
assert(type.dim()==dim);
}
//! \brief export type traits for function signature
typedef LocalBasisTraits<D,dim,Dune::FieldVector<D,dim>,R,1,Dune::FieldVector<R,1>,
Dune::FieldMatrix<R,1,dim>, 0> Traits;
//! \brief number of shape functions
unsigned int size () const
{
unsigned int dofs = 0;
for (unsigned int codim_i=0; codim_i<(dim+1); codim_i++)
dofs += ref_geo.size(codim_i);
return dofs;
}
//! \brief Evaluate all shape functions
inline void evaluateFunction (const typename Traits::DomainType& in,
std::vector<typename Traits::RangeType>& out) const
{
out.resize(size());
typename Traits::RangeType sum(0.);
unsigned int i = 0; /* counter */
/* loop over all codimensions but 0 */
for (unsigned int codim_i=dim; codim_i>0; codim_i--)
{
/* loop over all entities in codimension_i */
for (unsigned int subentity_i = 0; subentity_i < ref_geo.size(codim_i); ++subentity_i)
{
auto center = ref_geo.position(subentity_i,codim_i);
// double dist = std::min(0.5,(center - ref_geo.position(0,0)).two_norm());
double dist = 0.5;
center -= in;
/* Bump function in a radious of "dist". */
typename Traits::RangeType value = std::exp(-1./(dist*dist-center.two_norm2()))*std::exp(1./(dist*dist));
out[i] = (center.two_norm()<=dist) ? value : typename Traits::RangeType(0.);
sum += out[i];
i++;
}
}
/* Balance basis function in codimension 0 */
out[size()-1] =(typename Traits::RangeType(1.) - sum);
}
//! \brief Evaluate Jacobian of all shape functions
inline void
evaluateJacobian (const typename Traits::DomainType& in, // position
std::vector<typename Traits::JacobianType>& out) const // return value
{
out.resize(size());
// typename Traits::RangeType sum(0.);
// typename Traits::JacobianType sum(0.);
// loop in vertices
for (size_t i=0; i<ref_geo.size(dim); i++) {
auto center = ref_geo.position(i,dim);
center -= in;
double dist = 0.5; /* distance from the basis point to the oposite face in direction "center-in",
That, would do an affine transformation for these 1d functions */
// typename Traits::RangeType value = std::exp(-1./(dist*dist-center.two_norm2()))*std::exp(1./(dist*dist));
// value = (center.two_norm()<=dist) ? value : typename Traits::RangeType(0.);
// double value_jac = 0.;
// // if ((dist*dist-center.two_norm2())>1e-15)
// value_jac = 2.* center * value/((dist*dist-center.two_norm2())*(dist*dist-center.two_norm2()));
// out[i][0][0] = (center.two_norm()<=dist) ? value_jac : 0.;
// out[ref_geo.size(dim)][0][0] -= out[i][0][0];
}
}
/** \brief Evaluate partial derivatives of any order of all shape functions
* \param order Order of the partial derivatives, in the classic multi-index notation
* \param in Position where to evaluate the derivatives
* \param[out] out Return value: the desired partial derivatives
*/
inline void partial(const std::array<unsigned int,dim>& order,
const typename Traits::DomainType& in,
std::vector<typename Traits::RangeType>& out) const
{
// auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
out.resize(size());
// typename Traits::RangeType sum(0.);
for (int o = 0; o < order.size(); ++o)
{
if (order[o]==0)
{
evaluateFunction(in,out);
}
else
{
// typename Traits::JacobianType sum(0.);
// loop in vertices
for (size_t i=0; i<ref_geo.size(dim); i++) {
auto center = ref_geo.position(i,dim);
center -= in;
double dist = 0.5; /* distance from the basis point to the oposite face in direction "center-in",
That, would do an affine transformation for these 1d functions */
// typename Traits::RangeType value = std::exp(-1./(dist*dist-center.two_norm2()))*std::exp(1./(dist*dist));
// value = (center.two_norm()<=dist) ? value : typename Traits::RangeType(0.);
// double value_jac = 0.;
// // if ((dist*dist-center.two_norm2())>1e-15)
// value_jac = 2.* center * value/((dist*dist-center.two_norm2())*(dist*dist-center.two_norm2()));
// out[i] = (center.two_norm()<=dist) ? value_jac : 0.;
// out[ref_geo.size(dim)] -= out[i];
}
}
}
}
//! \brief Polynomial order of the shape functions
unsigned int order () const
{
return 0;
}
private:
const ReferenceElement<double,dim>& ref_geo;
};
}
#endif
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_BUMP_LOCALCOEFFICIENTS_HH
#define DUNE_BUMP_LOCALCOEFFICIENTS_HH
#include <cstddef>
#include <iostream>
#include <vector>
#include <dune/localfunctions/common/localkey.hh>
namespace Dune
{
/** \brief Local coefficients for simplex P1 elements
\nosubgrouping
\implements Dune::LocalCoefficientsVirtualImp
*/
template <int dimension>
class BumpLocalCoefficients
{
public:
//! \brief Standard constructor
enum {dim = dimension};
BumpLocalCoefficients(const GeometryType& type)
: ref_geo(ReferenceElements<double, dim>::general(type))
{
assert(type.dim()==dim);
fill_default();
}
//! number of coefficients
std::size_t size () const
{
unsigned int dofs = 0;
for (unsigned int codim_i=0; codim_i<(dim+1); codim_i++)
dofs += ref_geo.size(codim_i);
return dofs;
}
//! get i'th index
const LocalKey& localKey (std::size_t i) const
{
return li[i];
}
private:
void fill_default()
{
li.resize(size());
unsigned int i = 0;
/* loop over all codimensions but 0 */
for (unsigned int codim_i=dim; codim_i>0; codim_i--)
{
/* loop over all entities in codimension_i */
for (unsigned int subentity_i = 0; subentity_i < ref_geo.size(codim_i); ++subentity_i)
{
li[i] = LocalKey(subentity_i,codim_i,0);
i++;
}
}
li[size()-1] = LocalKey(0,0,0);
}
const ReferenceElement<double,dim>& ref_geo;
std::vector<LocalKey> li;
};
}
#endif
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_BUMPLOCALINTERPOLATION_HH
#define DUNE_BUMPLOCALINTERPOLATION_HH
#include <vector>
namespace Dune
{
template<class LB>
class BumpLocalInterpolation
{
public:
enum {dim = LB::Traits::dimDomain};
BumpLocalInterpolation(const GeometryType& type)
: ref_geo(ReferenceElements<double, dim>::general(type))
{}
template<typename F, typename C>
void interpolate (const F& f, std::vector<C>& out) const
{
typename LB::Traits::DomainType x;
typename LB::Traits::RangeType y;
typedef typename LB::Traits::DomainFieldType D;
// const int dofs = ref_geo.size(dim)+1;
// out.resize(dofs);
// for (int i=0; i<dofs-1; i++)
// {
// x = ref_geo.position(i,dim);
// f.evaluate(x,y);
// out[i] = y;
// }
// x[dofs-1] = ref_geo.position(0,0);
// f.evaluate(x,y);
// out[dofs-1] = y;
unsigned int dofs = 0;
for (unsigned int codim_i=0; codim_i<(dim+1); codim_i++)
dofs += ref_geo.size(codim_i);
out.resize(dofs);
unsigned int i = 0;
/* loop over all codimensions but 0 */
for (unsigned int codim_i=dim; codim_i>0; codim_i--)
{
/* loop over all entities in codimension_i */
for (unsigned int subentity_i = 0; subentity_i < ref_geo.size(codim_i); ++subentity_i)
{
x = ref_geo.position(subentity_i,codim_i);
f.evaluate(x,y);
out[i] = y;
i++;
}
}
x = ref_geo.position(0,0);
f.evaluate(x,y);
out[dofs-1] = y;
}
private:
const ReferenceElement<double,dim>& ref_geo;
};
}
#endif
dune_add_test(SOURCES bumplocal.cc)
\ No newline at end of file
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "config.h"
#include <iostream>
#include <ostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/localfunctions/bump/bumplocalbasis.hh>
#include <dune/localfunctions/bump/bumplocalcoefficients.hh>
#include <dune/localfunctions/bump/bumplocalinterpolation.hh>
#include <dune/localfunctions/bump/bumplocal.hh>
#include <dune/localfunctions/test/test-localfe.hh>
void testBumpLocalBasis(int N)
{
Dune::GeometryType type;
type.makeLine();
Dune::BumpLocalBasis<double,double,1> basis(type);
std::vector<double> vertexdata;
Dune::FieldVector<double,1> x;
std::vector<Dune::FieldVector<double,1>> y;
vertexdata.resize(N+1);
for (int i = 0; i <= N; ++i)
{
x = (double)i/N;
basis.evaluateFunction(x,y);
vertexdata.push_back(y[0]);
std::cout << x << "\t";
for (int j = 0; j < y.size(); ++j)
{
std::cout << y[j] << "\t";
}
std::cout << std::accumulate(y.begin(),y.end(),y[0]*0.) << "\t";
std::cout << std::endl;
}
}
void testBumpLocalCoefficients() {
Dune::GeometryType type;
type.makeLine();
Dune::BumpLocalCoefficients<1> coeff(type);
for (int i = 0; i < coeff.size(); ++i)
{
std::cout << coeff.localKey(i) << std::endl;
}
}
// void testBumpLocalInterpolation() {
// Dune::GeometryType type;
// type.makeLine();
// Dune::BumpLocalBasis<double,double,1> basis(type);
// Dune::BumpLocalInterpolation<Dune::BumpLocalBasis<double,double,1>> interpolation(type);
// std::vector<Dune::FieldVector<double,1>> y;
// interpolation.interpolate(basis,y);
// for (int i = 0; i < basis.size(); ++i)
// {
// std::cout << y[i] << std::endl;
// }
// }
void testBumpLocal() {
{
std::cout << "***** Line *****" << std::endl;
Dune::GeometryType type;
type.makeLine();
Dune::BumpLocalFiniteElement<double,double,1> bumpfe(type);
testFE(bumpfe);
}
{
std::cout << "***** Triangle *****" << std::endl;
Dune::GeometryType type;
type.makeTriangle();
Dune::BumpLocalFiniteElement<double,double,2> bumpfe(type);
testFE(bumpfe);
}
// {
// std::cout << "***** Quadrilateral *****" << std::endl;
// Dune::GeometryType type;
// type.makeQuadrilateral();
// Dune::BumpLocalFiniteElement<double,double,2> bumpfe(type);
// testFE(bumpfe);
// }
{
std::cout << "***** Tetrahedron *****" << std::endl;
Dune::GeometryType type;
type.makeTetrahedron();
Dune::BumpLocalFiniteElement<double,double,3> bumpfe(type);
testFE(bumpfe);
}
{
std::cout << "***** Pyramid *****" << std::endl;
Dune::GeometryType type;
type.makePyramid();
Dune::BumpLocalFiniteElement<double,double,3> bumpfe(type);
testFE(bumpfe);
}
{
std::cout << "***** Prism *****" << std::endl;
Dune::GeometryType type;
type.makePrism();
Dune::BumpLocalFiniteElement<double,double,3> bumpfe(type);
testFE(bumpfe);
}
{
std::cout << "***** Hexahedron *****" << std::endl;
Dune::GeometryType type;
type.makeHexahedron();
Dune::BumpLocalFiniteElement<double,double,3> bumpfe(type);
testFE(bumpfe);
}
}
int main(int argc, char** argv)
{
try
{
Dune::MPIHelper::instance(argc, argv);
std::cout << "Local Basis:" << std::endl;
testBumpLocalBasis(5);
std::cout << "Local Coefficients:" << std::endl;
testBumpLocalCoefficients();
// std::cout << "Local Interpolation:" << std::endl;
// testBumpLocalInterpolation();
std::cout << "Test Local FE:" << std::endl;
testBumpLocal();
} catch (std::exception &e) {
std::cerr << e.what() << std::endl;
return 1;
} catch (...) {
std::cerr << "Generic exception!" << std::endl;
return 2;
}
return 0;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment