Skip to content

GitLab

Commit 808eac8e authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Merge with Simons branch AMDiS_generic.

parent e43e8013
Hide whitespace changes
Inline Side-by-side
AMDiS/src/BasisFunction.cc
View file @ 808eac8e
......@@ -37,75 +37,12 @@ namespace AMDiS {
dow = Global::getGeo(WORLD);
}
BasisFunction::~BasisFunction()
{
delete nDOF;
}
/****************************************************************************/
/* some routines for evaluation of a finite element function, its gradient */
/* and second derivatives; all those with _fast use the preevaluated */
/* basis functions at that point. */
/****************************************************************************/
double BasisFunction::evalUh(const DimVec<double>& lambda,
const ElementVector& uh_loc) const
{
double val = 0.0;
for (int i = 0; i < nBasFcts; i++)
val += uh_loc[i] * (*(*phi)[i])(lambda);
return(val);
}
const WorldVector<double>& BasisFunction::evalUh(const DimVec<double>& lambda,
const mtl::dense_vector<WorldVector<double> >& uh_loc,
WorldVector<double>* values) const
{
static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
WorldVector<double> *val = (NULL != values) ? values : &Values;
for (int n = 0; n < dow; n++)
(*val)[n] = 0.0;
for (int i = 0; i < nBasFcts; i++) {
double phil = (*(*phi)[i])(lambda);
for (int n = 0; n < dow; n++)
(*val)[n] += uh_loc[i][n] * phil;
}
return((*val));
}
const WorldVector<double>& BasisFunction::evalGrdUh(const DimVec<double>& lambda,
const DimVec<WorldVector<double> >& grd_lambda,
const ElementVector& uh_loc,
WorldVector<double>* val) const
{
TEST_EXIT_DBG(val)("Return value is NULL\n");
mtl::dense_vector<double> grdTmp1(dim + 1);
mtl::dense_vector<double> grdTmp2(dim + 1, 0.0);
for (int i = 0; i < nBasFcts; i++) {
(*(*grdPhi)[i])(lambda, grdTmp1);
for (int j = 0; j < dim + 1; j++)
grdTmp2[j] += uh_loc[i] * grdTmp1[j];
}
for (int i = 0; i < dow; i++) {
(*val)[i] = 0.0;
for (int j = 0; j < dim + 1; j++)
(*val)[i] += grd_lambda[j][i] * grdTmp2[j];
}
return ((*val));
}
const WorldMatrix<double>& BasisFunction::evalD2Uh(const DimVec<double>& lambda,
const DimVec<WorldVector<double> >& grd_lambda,
......
AMDiS/src/BasisFunction.h
View file @ 808eac8e
......@@ -28,6 +28,7 @@
#include "Global.h"
#include "Boundary.h"
#include "MatrixVector.h"
#include "FixVec.h"
namespace AMDiS {
......@@ -311,17 +312,9 @@ namespace AMDiS {
* Evaluates elements value at barycentric coordinates lambda with local
* coefficient vector uh.
*/
double evalUh(const DimVec<double>& lambda, const ElementVector& uh) const;
template<typename T>
T evalUh(const DimVec<double>& lambda, const mtl::dense_vector<T>& uh) const;
/** \brief
* Evaluates elements value at barycentric coordinates lambda with local
* coefficient vector uh. If val is not NULL the result will be stored
* there, otherwise a pointer to a static local variable is returned which
* will be overwritten after the next call.
*/
const WorldVector<double>& evalUh(const DimVec<double>& lambda,
const mtl::dense_vector<WorldVector<double> >& uh,
WorldVector<double>* val) const;
/** \brief
* Evaluates the gradient at barycentric coordinates lambda. Lambda is the
......@@ -330,10 +323,12 @@ namespace AMDiS {
* there, otherwise a pointer to a static local variable is returned which
* will be overwritten after the next call.
*/
const WorldVector<double>& evalGrdUh(const DimVec<double>& lambda,
const DimVec<WorldVector<double> >& Lambda,
const ElementVector& uh,
WorldVector<double>* val) const;
template<typename T>
typename GradientType<T>::type& evalGrdUh(const DimVec<double>& lambda,
const DimVec<WorldVector<double> >& Lambda,
const mtl::dense_vector<T>& uh,
typename GradientType<T>::type& val) const;
/** \brief
* Evaluates the second derivative at barycentric coordinates lambda.
......@@ -377,5 +372,6 @@ namespace AMDiS {
};
}
#include "BasisFunction.hh"
#endif // AMDIS_BASISFUNCTION_H
AMDiS/src/DOFVector.cc
View file @ 808eac8e
......@@ -148,7 +148,7 @@ namespace AMDiS {
mtl::dense_vector<WorldVector<double> > uh(lambda.size());
for (int i = 0; i < lambda.size(); i++)
uh[i] = operator[](localIndices[i]);
basFcts->evalUh(lambda, uh, val);
*val = basFcts->evalUh(lambda, uh);
} else
throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
......@@ -301,322 +301,10 @@ namespace AMDiS {
template<>
DOFVector<WorldVector<double> >*
DOFVector<double>::getGradient(DOFVector<WorldVector<double> > *grad) const
{
FUNCNAME("DOFVector<double>::getGradient()");
// define result vector
static DOFVector<WorldVector<double> > *result = NULL;
if (grad) {
result = grad;
} else {
if(result && result->getFeSpace() != feSpace) {
delete result;
result = new DOFVector<WorldVector<double> >(feSpace, "gradient");
}
}
Mesh *mesh = feSpace->getMesh();
int dim = mesh->getDim();
const BasisFunction *basFcts = feSpace->getBasisFcts();
DOFAdmin *admin = feSpace->getAdmin();
// count number of nodes and dofs per node
std::vector<int> nNodeDOFs;
std::vector<int> nNodePreDofs;
std::vector<DimVec<double>*> bary;
int nNodes = 0;
int nDofs = 0;
for (int i = 0; i < dim + 1; i++) {
GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
int nPositions = mesh->getGeo(geoIndex);
int numPreDofs = admin->getNumberOfPreDofs(i);
for (int j = 0; j < nPositions; j++) {
int dofs = basFcts->getNumberOfDofs(geoIndex);
nNodeDOFs.push_back(dofs);
nDofs += dofs;
nNodePreDofs.push_back(numPreDofs);
}
nNodes += nPositions;
}
TEST_EXIT_DBG(nDofs == basFcts->getNumber())
("number of dofs != number of basis functions\n");
for (int i = 0; i < nDofs; i++)
bary.push_back(basFcts->getCoords(i));
ElementVector localUh(basFcts->getNumber());
// traverse mesh
std::vector<bool> visited(getUsedSize(), false);
TraverseStack stack;
Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
while (elInfo) {
const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
getLocalVector(elInfo->getElement(), localUh);
int localDOFNr = 0;
for (int i = 0; i < nNodes; i++) { // for all nodes
for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
if (!visited[dofIndex]) {
basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda,
localUh, &((*result)[dofIndex]));
visited[dofIndex] = true;
}
localDOFNr++;
}
}
elInfo = stack.traverseNext(elInfo);
}
return result;
}
template<>
DOFVector<WorldVector<double> >*
DOFVector<double>::getRecoveryGradient(DOFVector<WorldVector<double> > *grad) const
{
FUNCNAME("DOFVector<double>::getRecoveryGradient()");
// define result vector
static DOFVector<WorldVector<double> > *vec = NULL;
DOFVector<WorldVector<double> > *result = grad;
if (!result) {
if (vec && vec->getFeSpace() != feSpace) {
delete vec;
vec = NULL;
}
if (!vec)
vec = new DOFVector<WorldVector<double> >(feSpace, "gradient");
result = vec;
}
result->set(WorldVector<double>(DEFAULT_VALUE, 0.0));
DOFVector<double> volume(feSpace, "volume");
volume.set(0.0);
const BasisFunction *basFcts = feSpace->getBasisFcts();
int nBasisFcts = basFcts->getNumber();
DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
WorldVector<double> grd;
// traverse mesh
Mesh *mesh = feSpace->getMesh();
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
mtl::dense_vector<double> localUh(basFcts->getNumber());
std::vector<DegreeOfFreedom> localIndices(nBasisFcts);
while (elInfo) {
double det = elInfo->getDet();
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
getLocalVector(elInfo->getElement(), localUh);
basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
for (int i = 0; i < nBasisFcts; i++) {
(*result)[localIndices[i]] += grd * det;
volume[localIndices[i]] += det;
}
elInfo = stack.traverseNext(elInfo);
}
DOFVector<double>::Iterator volIt(&volume, USED_DOFS);
DOFVector<WorldVector<double> >::Iterator grdIt(result, USED_DOFS);
for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
if (*volIt != 0.0)
*grdIt *= 1.0 / (*volIt);
return result;
}
template<>
const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
WorldVector<double> *grdAtQPs) const
{
FUNCNAME("DOFVector<double>::getGrdAtQPs()");
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
if (quad && quadFast) {
TEST_EXIT_DBG(quad == quadFast->getQuadrature())
("quad != quadFast->quadrature\n");
}
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
int dow = Global::getGeo(WORLD);
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static WorldVector<double> *grd = NULL;
WorldVector<double> *result;
if (grdAtQPs) {
result = grdAtQPs;
} else {
if (grd)
delete [] grd;
grd = new WorldVector<double>[nPoints];
for (int i = 0; i < nPoints; i++)
grd[i].set(0.0);
result = grd;
}
ElementVector localVec(nBasFcts);
getLocalVector(elInfo->getElement(), localVec);
mtl::dense_vector<double> grd1(dim + 1);
int parts = Global::getGeo(PARTS, dim);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (quadFast) {
for (int i = 0; i < nPoints; i++) {
grd1 = 0.0;
for (int j = 0; j < nBasFcts; j++)
for (int k = 0; k < parts; k++)
grd1[k] += quadFast->getGradient(i, j, k) * localVec[j];
for (int l=0; l < dow; l++) {
result[i][l] = 0.0;
for (int k = 0; k < parts; k++)
result[i][l] += grdLambda[k][l] * grd1[k];
}
}
} else {
const BasisFunction *basFcts = feSpace->getBasisFcts();
mtl::dense_vector<double> grdPhi(dim + 1);
for (int i = 0; i < nPoints; i++) {
grd1 = 0.0;
for (int j = 0; j < nBasFcts; j++) {
(*(basFcts->getGrdPhi(j)))(quad->getLambda(i), grdPhi);
for (int k = 0; k < parts; k++)
grd1[k] += grdPhi[k] * localVec[j];
}
for (int l = 0; l < dow; l++) {
result[i][l] = 0.0;
for (int k = 0; k < parts; k++)
result[i][l] += grdLambda[k][l] * grd1[k];
}
}
}
return const_cast<const WorldVector<double>*>(result);
}
template<>
const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
WorldVector<double> *grdAtQPs) const
{
FUNCNAME("DOFVector<double>::getGrdAtQPs()");
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
if (quad && quadFast) {
TEST_EXIT_DBG(quad == quadFast->getQuadrature())("quad != quadFast->quadrature\n");
}
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
int dow = Global::getGeo(WORLD);
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static WorldVector<double> *grd = NULL;
WorldVector<double> *result;
if (grdAtQPs) {
result = grdAtQPs;
} else {
if (grd)
delete [] grd;
grd = new WorldVector<double>[nPoints];
for (int i = 0; i < nPoints; i++)
grd[i].set(0.0);
result = grd;
}
ElementVector localVec(nBasFcts);
getLocalVector(largeElInfo->getElement(), localVec);
const BasisFunction *basFcts = feSpace->getBasisFcts();
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
int parts = Global::getGeo(PARTS, dim);
const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();
mtl::dense_vector<double> grd1(dim + 1);
mtl::dense_vector<double> grdPhi(dim + 1);
mtl::dense_vector<double> tmp(dim + 1);
for (int i = 0; i < nPoints; i++) {
grd1 = 0.0;
for (int j = 0; j < nBasFcts; j++) {
grdPhi = 0.0;
for (int k = 0; k < nBasFcts; k++) {
(*(basFcts->getGrdPhi(k)))(quad->getLambda(i), tmp);
tmp *= m[j][k];
grdPhi += tmp;
}
for (int k = 0; k < parts; k++)
grd1[k] += grdPhi[k] * localVec[j];
}
for (int l = 0; l < dow; l++) {
result[i][l] = 0.0;
for (int k = 0; k < parts; k++)
result[i][l] += grdLambda[k][l] * grd1[k];
}
}
return const_cast<const WorldVector<double>*>(result);
}
template<>
const WorldMatrix<double> *DOFVectorBase<double>::getD2AtQPs(const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
WorldMatrix<double> *d2AtQPs) const
void DOFVectorBase<double>::getD2AtQPs( const ElInfo *elInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
mtl::dense_vector<D2Type<double>::type> &d2AtQPs) const
{
FUNCNAME("DOFVector<double>::getD2AtQPs()");
......@@ -634,28 +322,15 @@ namespace AMDiS {
int dow = Global::getGeo(WORLD);
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static WorldMatrix<double> *vec = NULL;
WorldMatrix<double> *result;
if (d2AtQPs) {
result = d2AtQPs;
} else {
if(vec) delete [] vec;
vec = new WorldMatrix<double>[nPoints];
for (int i = 0; i < nPoints; i++) {
vec[i].set(0.0);
}
result = vec;
}
ElementVector localVec(nBasFcts);
mtl::dense_vector<double> localVec(nBasFcts);
getLocalVector(el, localVec);
DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0);
int parts = Global::getGeo(PARTS, dim);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
d2AtQPs.change_dim(nPoints);
if (quadFast) {
for (int iq = 0; iq < nPoints; iq++) {
for (int k = 0; k < parts; k++)
......@@ -670,10 +345,10 @@ namespace AMDiS {
for (int i = 0; i < dow; i++)
for (int j = 0; j < dow; j++) {
result[iq][i][j] = 0.0;
d2AtQPs[iq][i][j] = 0.0;
for (int k = 0; k < parts; k++)
for (int l = 0; l < parts; l++)
result[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
d2AtQPs[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
}
}
} else {
......@@ -696,15 +371,13 @@ namespace AMDiS {
for (int i = 0; i < dow; i++)
for (int j = 0; j < dow; j++) {
result[iq][i][j] = 0.0;
d2AtQPs[iq][i][j] = 0.0;
for (int k = 0; k < parts; k++)
for (int l = 0; l < parts; l++)
result[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
d2AtQPs[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
}
}
}
return const_cast<const WorldMatrix<double>*>(result);
}
......@@ -827,7 +500,7 @@ namespace AMDiS {
for (int i = 0; i < nBasFcts; i++) {
if (vec[myLocalIndices[i]] == nul) {
coords = basFcts->getCoords(i);
vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs, NULL) * factor;
vec[myLocalIndices[i]] += vBasFcts->evalUh(*coords, vLocalCoeffs) * factor;
}
}
elInfo = stack.traverseNext(elInfo);
......@@ -916,7 +589,7 @@ namespace AMDiS {
DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
if (!visited[dofIndex]) {
basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, &grd);
basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda, localUh, grd);
for (int k = 0; k < dow; k++)
(*(*result)[k])[dofIndex] = grd[k];
......@@ -970,50 +643,6 @@ namespace AMDiS {
}
template<>
void DOFVectorBase<double>::getVecAtQPs(const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
const Quadrature *quad,
const FastQuadrature *quadFast,
mtl::dense_vector<double>& vecAtQPs) const
{
FUNCNAME("DOFVector<double>::getVecAtQPs()");
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
if (quad && quadFast) {
TEST_EXIT_DBG(quad == quadFast->getQuadrature())
("quad != quadFast->quadrature\n");
}
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
if (smallElInfo->getMesh() == feSpace->getMesh())
return getVecAtQPs(smallElInfo, quad, quadFast, vecAtQPs);
const BasisFunction *basFcts = feSpace->getBasisFcts();
int nPoints =
quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
vecAtQPs.change_dim(nPoints);
ElementVector localVec(nBasFcts);
getLocalVector(largeElInfo->getElement(), localVec);
mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
for (int i = 0; i < nPoints; i++) {
vecAtQPs[i] = 0.0;
for (int j = 0; j < nBasFcts; j++) {
double val = 0.0;
for (int k = 0; k < nBasFcts; k++)
val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
vecAtQPs[i] += localVec[j] * val;
}
}
}
template<>
void DOFVectorBase<double>::assemble(double factor, ElInfo *elInfo,
......
AMDiS/src/DOFVector.h
View file @ 808eac8e
......@@ -87,21 +87,46 @@ namespace AMDiS {
const FastQuadrature *quadFast,
mtl::dense_vector<T>& vecAtQPs) const;
const WorldVector<T> *getGrdAtQPs(const ElInfo *elInfo,
const Quadrature *quad,