Skip to content

GitLab

Commit 561814e2 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* This and that

parent 145602ce
Hide whitespace changes
Inline Side-by-side
Showing with 239 additions and 184 deletions
AMDiS/src/DOFVector.cc
View file @ 561814e2
......@@ -149,12 +149,12 @@ namespace AMDiS {
DOFVector<WorldVector<double> > *result = grad;
if(!result) {
if(vec && vec->getFESpace() != feSpace) {
if (!result) {
if (vec && vec->getFESpace() != feSpace) {
DELETE vec;
vec = NULL;
}
if(!vec) {
if (!vec) {
vec = NEW DOFVector<WorldVector<double> >(feSpace, "gradient");
}
result = vec;
......@@ -165,36 +165,27 @@ namespace AMDiS {
DOFVector<double> volume(feSpace, "volume");
volume.set(0.0);
Mesh *mesh = feSpace->getMesh();
int dim = mesh->getDim();
const BasisFunction *basFcts = feSpace->getBasisFcts();
DOFAdmin *admin = feSpace->getAdmin();
int numPreDOFs = admin->getNumberOfPreDOFs(0);
int nPreDOFs = feSpace->getAdmin()->getNumberOfPreDOFs(0);
DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
WorldVector<double> grd;
// traverse mesh
Mesh *mesh = feSpace->getMesh();
TraverseStack stack;
Flag fillFlag =
Mesh::CALL_LEAF_EL | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
ElInfo *elInfo = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
while (elInfo) {
double det = elInfo->getDet();
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const double *localUh = getLocalVector(elInfo->getElement(), NULL);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
WorldVector<double> grd;
basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
for (int i = 0; i < dim + 1; i++) {
DegreeOfFreedom dofIndex = dof[i][numPreDOFs];
DegreeOfFreedom dofIndex = dof[i][nPreDOFs];
(*result)[dofIndex] += grd * det;
volume[dofIndex] += det;
}
......@@ -207,7 +198,7 @@ namespace AMDiS {
for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt) {
if (*volIt != 0.0) {
*grdIt *= 1.0/(*volIt);
*grdIt *= 1.0 / (*volIt);
}
}
......@@ -232,12 +223,9 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
Element *el = elInfo->getElement();
int dim = elInfo->getMesh()->getDim();
int dow = Global::getGeo(WORLD);
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num();
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static WorldVector<double> *grd = NULL;
WorldVector<double> *result;
......@@ -247,22 +235,22 @@ namespace AMDiS {
if (grd)
delete [] grd;
grd = new WorldVector<double>[numPoints];
for (int i = 0; i < numPoints; i++) {
grd = new WorldVector<double>[nPoints];
for (int i = 0; i < nPoints; i++) {
grd[i].set(0.0);
}
result = grd;
}
double *localVec = new double[nBasFcts];
getLocalVector(el, localVec);
double *localVec = localVectors[myRank];
getLocalVector(elInfo->getElement(), localVec);
DimVec<double> grd1(dim, NO_INIT);
DimVec<double> &grd1 = *grdTmp[myRank];
int parts = Global::getGeo(PARTS, dim);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (quadFast) {
for (int i = 0; i < numPoints; i++) {
for (int i = 0; i < nPoints; i++) {
for (int j = 0; j < dim + 1; j++) {
grd1[j] = 0.0;
}
......@@ -281,10 +269,10 @@ namespace AMDiS {
}
}
} else {
// DimVec<double> grdPhi(dim, DEFAULT_VALUE, 0.0);
DimVec<double>* grdPhi = grdPhis[omp_get_thread_num()];
const BasisFunction *basFcts = feSpace->getBasisFcts();
DimVec<double>* grdPhi = grdPhis[myRank];
for (int i = 0; i < numPoints; i++) {
for (int i = 0; i < nPoints; i++) {
for (int j = 0; j < dim + 1; j++) {
grd1[j] = 0.0;
}
......@@ -305,8 +293,6 @@ namespace AMDiS {
}
}
delete [] localVec;
return const_cast<const WorldVector<double>*>(result);
}
......@@ -330,31 +316,25 @@ namespace AMDiS {
Element *el = elInfo->getElement();
int dim = elInfo->getMesh()->getDim();
int myRank = omp_get_thread_num();
int dow = Global::getGeo(WORLD);
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numPoints = quadrature->getNumPoints();
int i, j, k, l, iq;
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>[numPoints];
for(i = 0; i < numPoints; i++) {
vec = new WorldMatrix<double>[nPoints];
for (int i = 0; i < nPoints; i++) {
vec[i].set(0.0);
}
result = vec;
}
double *localVec = new double[nBasFcts];
double *localVec = localVectors[myRank];
getLocalVector(el, localVec);
DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0);
......@@ -362,56 +342,53 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (quadFast) {
for (iq = 0; iq < numPoints; iq++) {
for (k = 0; k < parts; k++)
for (l = 0; l < parts; l++)
for (int iq = 0; iq < nPoints; iq++) {
for (int k = 0; k < parts; k++)
for (int l = 0; l < parts; l++)
D2Tmp[k][l] = 0.0;
for (i = 0; i < nBasFcts; i++) {
for (k = 0; k < parts; k++)
for (l = 0; l < parts; l++)
D2Tmp[k][l] += localVec[i]* quadFast->getSecDer(iq, i, k, l);
for (int i = 0; i < nBasFcts; i++) {
for (int k = 0; k < parts; k++)
for (int l = 0; l < parts; l++)
D2Tmp[k][l] += localVec[i] * quadFast->getSecDer(iq, i, k, l);
}
for (i = 0; i < dow; i++)
for (j = 0; j < dow; j++) {
for (int i = 0; i < dow; i++)
for (int j = 0; j < dow; j++) {
result[iq][i][j] = 0.0;
for (k = 0; k < parts; k++)
for (l = 0; l < parts; l++)
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];
}
}
} else {
// DimMat<double> D2Phi(dim, NO_INIT);
const BasisFunction *basFcts = feSpace->getBasisFcts();
DimMat<double>* D2Phi = D2Phis[omp_get_thread_num()];
for (iq = 0; iq < numPoints; iq++) {
for (k = 0; k < parts; k++)
for (l = 0; l < parts; l++)
for (int iq = 0; iq < nPoints; iq++) {
for (int k = 0; k < parts; k++)
for (int l = 0; l < parts; l++)
D2Tmp[k][l] = 0.0;
for (i = 0; i < nBasFcts; i++) {
for (int i = 0; i < nBasFcts; i++) {
WARNING("not tested after index correction\n");
(*(basFcts->getD2Phi(i)))(quad->getLambda(iq), *D2Phi);
for (k = 0; k < parts; k++)
for (l = 0; l < parts; l++)
for (int k = 0; k < parts; k++)
for (int l = 0; l < parts; l++)
D2Tmp[k][l] += localVec[i] * (*D2Phi)[k][l];
}
for (i = 0; i < dow; i++)
for (j = 0; j < dow; j++) {
for (int i = 0; i < dow; i++)
for (int j = 0; j < dow; j++) {
result[iq][i][j] = 0.0;
for (k = 0; k < parts; k++)
for (l = 0; l < parts; l++)
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];
}
}
}
delete [] localVec;
return const_cast<const WorldMatrix<double>*>(result);
}
......
AMDiS/src/DOFVector.h
View file @ 561814e2
......@@ -231,6 +231,17 @@ namespace AMDiS {
*/
std::vector<DegreeOfFreedom*> localIndices;
/** \brief
* Are used to store temporary local values of an element.
* Thread safe.
*/
std::vector<T*> localVectors;
/** \brief
* Temporarly used in \ref getGrdAtQPs. Thread safe.
*/
std::vector<DimVec<double>*> grdTmp;
/** \brief
* Temporarly used in \ref getGrdAtQPs. Thread safe.
*/
......@@ -240,6 +251,11 @@ namespace AMDiS {
* Temporarly used in \ref getD2AtQPs. Thread safe.
*/
std::vector<DimMat<double>*> D2Phis;
/** \brief
* Dimension of the mesh this DOFVectorBase belongs to
*/
int dim;
};
......
AMDiS/src/DOFVector.hh
View file @ 561814e2
......@@ -29,15 +29,19 @@ namespace AMDiS {
boundaryManager(NULL)
{
nBasFcts = feSpace->getBasisFcts()->getNumber();
int dim = feSpace->getMesh()->getDim();
dim = feSpace->getMesh()->getDim();
localIndices.resize(omp_get_overall_max_threads());
localVectors.resize(omp_get_overall_max_threads());
grdPhis.resize(omp_get_overall_max_threads());
grdTmp.resize(omp_get_overall_max_threads());
D2Phis.resize(omp_get_overall_max_threads());
for (int i = 0; i < omp_get_overall_max_threads(); i++) {
localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts);
localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts);
localVectors[i] = GET_MEMORY(T, this->nBasFcts);
grdPhis[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
grdTmp[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
D2Phis[i] = NEW DimMat<double>(dim, NO_INIT);
}
}
......@@ -47,7 +51,9 @@ namespace AMDiS {
{
for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
FREE_MEMORY(localIndices[i], DegreeOfFreedom, this->nBasFcts);
FREE_MEMORY(localVectors[i], T, this->nBasFcts);
DELETE grdPhis[i];
DELETE grdTmp[i];
DELETE D2Phis[i];
}
}
......@@ -540,26 +546,19 @@ namespace AMDiS {
return 0;
}
// das hat Andreas eingefuegt: Integral...
template<typename T>
double DOFVector<T>::Int(Quadrature* q) const
{
FUNCNAME("DOFVector::Int");
FUNCNAME("DOFVector::Int()");
Mesh* mesh = feSpace->getMesh();
int deg;
dim = mesh->getDim();
if (!q)
{
deg = 2*feSpace->getBasisFcts()->getDegree();
q = Quadrature::provideQuadrature(dim, deg);
}
if (!q) {
int deg = 2 * feSpace->getBasisFcts()->getDegree();
q = Quadrature::provideQuadrature(this->dim, deg);
}
quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *q, INIT_PHI);
norm = 0.0;
traverseVector = const_cast<DOFVector<T>*>(this);
......@@ -591,28 +590,20 @@ namespace AMDiS {
return 0;
}
// ... und die L1-Norm ...
template<typename T>
double DOFVector<T>::L1Norm(Quadrature* q) const
{
FUNCNAME("DOFVector::L1Norm");
FUNCNAME("DOFVector::L1Norm()");
Mesh* mesh = feSpace->getMesh();
int deg;
dim = mesh->getDim();
if (!q)
{
deg = 2*feSpace->getBasisFcts()->getDegree();
q = Quadrature::provideQuadrature(dim, deg);
}
if (!q) {
int deg = 2 * feSpace->getBasisFcts()->getDegree();
q = Quadrature::provideQuadrature(this->dim, deg);
}
quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),*q,INIT_PHI);
norm = 0.0;
traverseVector = const_cast<DOFVector<T>*>(this);
Mesh* mesh = feSpace->getMesh();
mesh->traverse(-1,
Mesh::CALL_LEAF_EL|
......@@ -643,23 +634,15 @@ namespace AMDiS {
}
// bis hierhin gehen Andreas Ergaenzungen...
template<typename T>
double DOFVector<T>::L2NormSquare(Quadrature* q) const
{
FUNCNAME("DOFVector::L2NormSquare");
Mesh* mesh = feSpace->getMesh();
FUNCNAME("DOFVector::L2NormSquare()");
int deg;
dim = mesh->getDim();
if (!q)
{
deg = 2*feSpace->getBasisFcts()->getDegree();
q = Quadrature::provideQuadrature(dim, deg);
}
if (!q) {
int deg = 2 * feSpace->getBasisFcts()->getDegree();
q = Quadrature::provideQuadrature(this->dim, deg);
}
quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),
*q,
......@@ -667,8 +650,8 @@ namespace AMDiS {
norm = 0.0;
traverseVector = const_cast<DOFVector<T>*>(this);
mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS|Mesh::FILL_DET, L2Norm_fct);
Mesh* mesh = feSpace->getMesh();
mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, L2Norm_fct);
return norm;
}
......@@ -722,24 +705,19 @@ namespace AMDiS {
template<typename T>
double DOFVector<T>::H1NormSquare(Quadrature *q) const
{
FUNCNAME("DOFVector::H1NormSquare");
int deg;
Mesh *mesh = feSpace->getMesh();
dim = mesh->getDim();
FUNCNAME("DOFVector::H1NormSquare()");
if (!q)
{
deg = 2*feSpace->getBasisFcts()->getDegree()-2;
q = Quadrature::provideQuadrature(dim, deg);
}
quad_fast =
FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),
*q,
INIT_GRD_PHI);
if (!q) {
int deg = 2 * feSpace->getBasisFcts()->getDegree() - 2;
q = Quadrature::provideQuadrature(this->dim, deg);
}
quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),
*q, INIT_GRD_PHI);
norm = 0.0;
traverseVector = const_cast<DOFVector<T>*>(this);
traverseVector = const_cast<DOFVector<T>*>(this);
Mesh *mesh = feSpace->getMesh();
mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA,
......@@ -752,10 +730,9 @@ namespace AMDiS {
void DOFVector<T>::compressDOFIndexed(int first, int last,
std::vector<DegreeOfFreedom> &newDOF)
{
int i, j;
for(i = first; i <= last; i++) {
if((j = newDOF[i]) >= 0) {
vec[j] = vec[i];
for (int i = first; i <= last; i++) {
if (newDOF[i] >= 0) {
vec[newDOF[i]] = vec[i];
}
}
}
......@@ -1252,15 +1229,9 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
Element *el = elInfo->getElement();
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numPoints = quadrature->getNumPoints();
int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static T *localvec = NULL;
T *result;
if (vecAtQPs) {
......@@ -1268,17 +1239,17 @@ namespace AMDiS {
} else {
if (localvec)
delete [] localvec;
localvec = new T[numPoints];
for (int i = 0; i < numPoints; i++) {
localvec = new T[nPoints];
for (int i = 0; i < nPoints; i++) {
localvec[i] = 0.0;
}
result = localvec;
}
T *localVec = new T[nBasFcts];
getLocalVector(el, localVec);
T *localVec = localVectors[omp_get_thread_num()];
getLocalVector(elInfo->getElement(), localVec);
for (int i = 0; i < numPoints; i++) {
for (int i = 0; i < nPoints; i++) {
result[i] = 0.0;
for (int j = 0; j < nBasFcts; j++) {
result[i] += localVec[j] *
......@@ -1287,9 +1258,7 @@ namespace AMDiS {
((*(basFcts->getPhi(j)))(quad->getLambda(i))));
}
}
delete [] localVec;
return const_cast<const T*>(result);
}
......@@ -1299,23 +1268,18 @@ namespace AMDiS {
{
FUNCNAME("DOFVector::DoubleWell()");
Mesh* mesh = feSpace->getMesh();
int deg = 0;
dim = mesh->getDim();
if (!q) {
deg = 2 * feSpace->getBasisFcts()->getDegree();
q = Quadrature::provideQuadrature(dim, deg);
int deg = 2 * feSpace->getBasisFcts()->getDegree();
q = Quadrature::provideQuadrature(this->dim, deg);
}
quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),
*q, INIT_PHI);
norm = 0.0;
traverseVector = const_cast<DOFVector<T>*>(this);
mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS|Mesh::FILL_DET, DoubleWell_fct);
Mesh* mesh = feSpace->getMesh();
mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET,
DoubleWell_fct);
return norm;
}
......@@ -1323,15 +1287,13 @@ namespace AMDiS {
template<typename T>
int DOFVector<T>::DoubleWell_fct(ElInfo *elinfo)
{
double det = elinfo->getDet();
const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL,
quad_fast, NULL);
int numPoints = quad_fast->getNumPoints();
const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
int nPoints = quad_fast->getNumPoints();
double normT = 0.0;
for (int iq = 0; iq < numPoints; iq++) {
normT += quad_fast->getWeight(iq)*sqr(uh_vec[iq])*sqr(1.0 - uh_vec[iq]);
for (int iq = 0; iq < nPoints; iq++) {
normT += quad_fast->getWeight(iq) * sqr(uh_vec[iq]) * sqr(1.0 - uh_vec[iq]);
}
norm += det * normT;
norm += elinfo->getDet() * normT;
return 0;
}
......
AMDiS/src/ElInfo3d.cc
View file @ 561814e2
......@@ -609,4 +609,63 @@ namespace AMDiS {
}
}
void ElInfo3d::getRefSimplexCoords(DimMat<double> *coords) const
{
(*coords)[0][0] = 1.0;
(*coords)[1][0] = 0.0;
(*coords)[2][0] = 0.0;
(*coords)[3][0] = 0.0;
(*coords)[0][1] = 0.0;
(*coords)[1][1] = 1.0;
(*coords)[2][1] = 0.0;
(*coords)[3][1] = 0.0;
(*coords)[0][2] = 0.0;
(*coords)[1][2] = 0.0;
(*coords)[2][2] = 1.0;
(*coords)[3][2] = 0.0;
(*coords)[0][3] = 0.0;
(*coords)[1][3] = 0.0;
(*coords)[2][3] = 0.0;
(*coords)[3][3] = 1.0;
}
void ElInfo3d::getSubElementCoords(DimMat<double> *coords,
int iChild) const
{
FUNCNAME("ElInfo3d::getSubElementCoords()");
ERROR_EXIT("Not yet\n");
double c0 = ((*coords)[0][0] + (*coords)[0][1]) * 0.5;
double c1 = ((*coords)[1][0] + (*coords)[1][1]) * 0.5;