Skip to content

GitLab

Commit 6fd5f942 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Work on pdd.

parent f1d2d34b
Hide whitespace changes
Inline Side-by-side
Showing with 1883 additions and 2009 deletions
AMDiS/src/BasisFunction.h
View file @ 6fd5f942
......@@ -278,17 +278,22 @@ namespace AMDiS {
{}
/// Returns local dof indices of the element for the given fe space.
virtual const DegreeOfFreedom *getLocalIndices(const Element*,
const DOFAdmin*,
DegreeOfFreedom*) const
virtual const DegreeOfFreedom *getLocalIndices(const Element *el,
const DOFAdmin *admin,
DegreeOfFreedom *dofPtr) const
{
return NULL;
}
/// Returns local dof indices of the element for the given fe space.
virtual void getLocalIndicesVec(const Element*,
const DOFAdmin*,
Vector<DegreeOfFreedom>*) const
virtual void getLocalIndicesVec(const Element *el,
const DOFAdmin *admin,
Vector<DegreeOfFreedom> *ve) const
{}
virtual void getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin,
std::vector<const DegreeOfFreedom*>& vec) const
{}
......
AMDiS/src/ElementDofIterator.cc
View file @ 6fd5f942
......@@ -2,16 +2,20 @@
#include "Mesh.h"
#include "DOFAdmin.h"
#include "Element.h"
#include "BasisFunction.h"
namespace AMDiS {
void ElementDofIterator::reset(Element* element)
void ElementDofIterator::reset(Element* el)
{
FUNCNAME("ElementDofIterator::reset()");
TEST_EXIT_DBG(el->getMesh() == mesh)
("Mesh and element does not fit together!\n");
TEST_EXIT_DBG(el)("No element!\n");
element = el;
dofs = element->getDOF();
mesh = element->getMesh();
dim = mesh->getDim();
// Start with vertices.
pos = 0;
......@@ -31,6 +35,9 @@ namespace AMDiS {
node0 = mesh->getNode(posIndex);
// Get number of vertices in this dimension.
nElements = Global::getGeo(posIndex, mesh->getDim());
if (inOrder)
orderPosition = basisFcts->orderOfPositionIndices(element, posIndex, 0);
}
bool ElementDofIterator::next()
......@@ -73,15 +80,28 @@ namespace AMDiS {
// Get first dof index position for the geo index position.
node0 = mesh->getNode(posIndex);
if (inOrder)
orderPosition = basisFcts->orderOfPositionIndices(element, posIndex, 0);
} else {
// That's all, we jave traversed all dofs of the mesh element.
return false;
}
} else {
if (inOrder)
orderPosition =
basisFcts->orderOfPositionIndices(element, posIndex, elementPos);
}
}
return true;
}
bool ElementDofIterator::nextStrict()
{
dofPos = nDofs;
return next();
}
}
AMDiS/src/ElementDofIterator.h
View file @ 6fd5f942
......@@ -24,6 +24,7 @@
#include "AMDiS_fwd.h"
#include "Global.h"
#include "Mesh.h"
namespace AMDiS {
......@@ -44,8 +45,12 @@ namespace AMDiS {
{
public:
/// Constructor.
ElementDofIterator(const DOFAdmin* dofAdmin)
: admin(dofAdmin)
ElementDofIterator(const FiniteElemSpace* feSpace, bool inOrderPos = false)
: admin(feSpace->getAdmin()),
basisFcts(feSpace->getBasisFcts()),
mesh(feSpace->getMesh()),
dim(mesh->getDim()),
inOrder(inOrderPos)
{}
/// Start a new traverse with the given element.
......@@ -54,35 +59,50 @@ namespace AMDiS {
/// Go to next dof. Returns false, if there is dof anymore.
bool next();
bool nextStrict();
/// Returns the dof index of the current dof.
const DegreeOfFreedom getDof()
inline const DegreeOfFreedom getDof()
{
return dofs[node0 + elementPos][n0 + dofPos];
if (inOrder)
return dofs[node0 + elementPos][n0 + orderPosition[dofPos]];
else
return dofs[node0 + elementPos][n0 + dofPos];
}
/// Returns a pointer to the current dof.
const DegreeOfFreedom* getDofPtr()
inline const DegreeOfFreedom* getDofPtr()
{
return &dofs[node0 + elementPos][n0 + dofPos];
if (inOrder)
return &dofs[node0 + elementPos][n0 + orderPosition[dofPos]];
else
return &dofs[node0 + elementPos][n0 + dofPos];
}
/// Returns \ref pos, the current position (vertex, edge, face) of the traverse.
int getCurrentPos()
inline int getCurrentPos()
{
return pos;
}
/// Returns \ref elementPos, the number of vertex, edge or face that is traversed.
int getCurrentElementPos()
inline int getCurrentElementPos()
{
return elementPos;
}
inline GeoIndex getPosIndex()
{
return posIndex;
}
protected:
/// The dof admin for which dofs should be traversed.
const DOFAdmin* admin;
const BasisFunction* basisFcts;
/// Pointer to the dofs that should be traversed.
const DegreeOfFreedom **dofs;
......@@ -92,6 +112,12 @@ namespace AMDiS {
/// Dimension of the mesh.
int dim;
bool inOrder;
int* orderPosition;
Element* element;
/// Current position (i.e., vertex, edge, face) of the traverse.
int pos;
......
AMDiS/src/InteriorBoundary.h
View file @ 6fd5f942
......@@ -70,13 +70,16 @@ namespace AMDiS {
* the classical domain decomposition parallelization.
*/
class InteriorBoundary {
public:
typedef std::map<int, std::vector<AtomicBoundary> > RankToBoundMap;
public:
InteriorBoundary() {}
AtomicBoundary& getNewAtomicBoundary(int rank);
public:
std::map<int, std::vector<AtomicBoundary> > boundary;
RankToBoundMap boundary;
};
}
......
AMDiS/src/Lagrange.cc
View file @ 6fd5f942
......@@ -44,7 +44,6 @@ namespace AMDiS {
// set function pointer
setFunctionPointer();
}
Lagrange::~Lagrange()
......@@ -60,11 +59,9 @@ namespace AMDiS {
Lagrange* Lagrange::getLagrange(int dim, int degree)
{
std::list<Lagrange*>::iterator it;
for (it = allBasFcts.begin(); it != allBasFcts.end(); it++) {
if (((*it)->dim == dim) && ((*it)->degree == degree)) {
for (it = allBasFcts.begin(); it != allBasFcts.end(); it++)
if (((*it)->dim == dim) && ((*it)->degree == degree))
return (*it);
}
}
Lagrange* newLagrange = new Lagrange(dim, degree);
allBasFcts.push_back(newLagrange);
......@@ -113,7 +110,7 @@ namespace AMDiS {
grdPhi = &grdPhiDimDegree[dim][degree];
d2Phi = &D2PhiDimDegree[dim][degree];
switch(degree) {
switch (degree) {
case 0:
refineInter_fct = refineInter0;
coarseRestr_fct = coarseRestr0;
......@@ -125,7 +122,7 @@ namespace AMDiS {
coarseInter_fct = NULL; // not yet implemented
break;
case 2:
switch(dim) {
switch (dim) {
case 1:
refineInter_fct = refineInter2_1d;
coarseRestr_fct = NULL; // not yet implemented
......@@ -145,7 +142,7 @@ namespace AMDiS {
}
break;
case 3:
switch(dim) {
switch (dim) {
case 1:
refineInter_fct = refineInter3_1d;
coarseRestr_fct = coarseRestr3_1d;
......@@ -165,7 +162,7 @@ namespace AMDiS {
}
break;
case 4:
switch(dim) {
switch (dim) {
case 1:
refineInter_fct = refineInter4_1d;
coarseRestr_fct = coarseRestr4_1d;
......@@ -194,11 +191,10 @@ namespace AMDiS {
if (static_cast<int>(baryDimDegree[dim][degree].size()) == 0) {
ndofDimDegree[dim][degree] = new DimVec<int>(dim, DEFAULT_VALUE, 0);
if (degree != 0) {
if (degree != 0)
(*ndofDimDegree[dim][degree])[VERTEX] = 1;
} else {
(*ndofDimDegree[dim][degree])[VERTEX] = 0;
}
else
(*ndofDimDegree[dim][degree])[VERTEX] = 0;
for (int i = 1; i < dim + 1; i++) {
nBasFcts = getNumberOfDOFs(i, degree);
......@@ -224,7 +220,10 @@ namespace AMDiS {
GeoIndex position, int positionIndex, int nodeIndex,
int** vertices)
{
FUNCNAME("Lagrange::setVertices()");
TEST_EXIT_DBG((*vertices) == NULL)("vertices != NULL\n");
int dimOfPosition = DIM_OF_INDEX(position, dim);
*vertices = new int[dimOfPosition + 1];
......@@ -261,14 +260,15 @@ namespace AMDiS {
int nodeIndex)
: vertices(NULL)
{
FUNCNAME("Lagrange::Phi::Phi")
// get relevant vertices
Lagrange::setVertices(owner->getDim(),
owner->getDegree(),
position,
positionIndex,
nodeIndex,
&vertices);
FUNCNAME("Lagrange::Phi::Phi()");
// get relevant vertices
Lagrange::setVertices(owner->getDim(),
owner->getDegree(),
position,
positionIndex,
nodeIndex,
&vertices);
// set function pointer
switch (owner->getDegree()) {
......@@ -354,7 +354,7 @@ namespace AMDiS {
case CENTER:
switch (owner->getDim()) {
case 1:
if(nodeIndex == 1)
if (nodeIndex == 1)
func = phi4e1;
else
func = phi4e0;
......@@ -506,7 +506,8 @@ namespace AMDiS {
}
}
Lagrange::GrdPhi::~GrdPhi() {
Lagrange::GrdPhi::~GrdPhi()
{
delete [] vertices;
}
......@@ -596,7 +597,7 @@ namespace AMDiS {
func = D2Phi4v;
break;
case EDGE:
if(nodeIndex == 1)
if (nodeIndex == 1)
func = D2Phi4e1;
else
func = D2Phi4e0;
......@@ -633,19 +634,16 @@ namespace AMDiS {
}
}
Lagrange::D2Phi::~D2Phi() {
Lagrange::D2Phi::~D2Phi()
{
delete [] vertices;
}
void Lagrange::createCoords(int* coordInd,
int numCoords,
int dimIndex,
int rest,
void Lagrange::createCoords(int* coordInd, int numCoords, int dimIndex, int rest,
DimVec<double>* vec)
{
if (vec == NULL) {
if (vec == NULL)
vec = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
}
if (dimIndex == numCoords - 1) {
(*vec)[coordInd[dimIndex]] = double(rest) / degree;
......@@ -670,7 +668,7 @@ namespace AMDiS {
int *coordInd = new int[i + 1]; // indices of relevant coords
for (int k = 0; k < i + 1; k++) {
coordInd[k] = Global::getReferenceElement(dim)->
getVertexOfPosition(INDEX_OF_DIM(i,dim), j, k);
getVertexOfPosition(INDEX_OF_DIM(i, dim), j, k);
}
createCoords(coordInd, i + 1, 0, degree);
delete [] coordInd;
......@@ -715,37 +713,32 @@ namespace AMDiS {
int dimOfPosition = DIM_OF_INDEX(position, dim);
// vertex
if (dimOfPosition == 0) {
return &sortedVertex;
}
if (dimOfPosition == 0)
return &sortedVertex;
// edge
if ((dimOfPosition == 1) && (degree == 2)) {
return &sortedEdgeDeg2;
}
if ((dimOfPosition == 1) && (degree == 2))
return &sortedEdgeDeg2;
int vertex[3];
int** dof = const_cast<int**>(el->getDOF());
int verticesOfPosition = dimOfPosition + 1;
for (int i = 0; i < verticesOfPosition; i++) {
for (int i = 0; i < verticesOfPosition; i++)
vertex[i] = Global::getReferenceElement(dim)->
getVertexOfPosition(position, positionIndex, i);
}
getVertexOfPosition(position, positionIndex, i);
if (dimOfPosition == 1) {
if (degree == 3) {
if (dof[vertex[0]][0] < dof[vertex[1]][0]) {
if (dof[vertex[0]][0] < dof[vertex[1]][0])
return sortedEdgeDeg3[0];
} else {
return sortedEdgeDeg3[1];
}
else
return sortedEdgeDeg3[1];
} else { // degree == 4
if (dof[vertex[0]][0] < dof[vertex[1]][0]) {
if (dof[vertex[0]][0] < dof[vertex[1]][0])
return sortedEdgeDeg4[0];
} else {
return sortedEdgeDeg4[1];
}
else
return sortedEdgeDeg4[1];
}
}
......@@ -774,9 +767,8 @@ namespace AMDiS {
}
// center
if ((dimOfPosition == 3) && (degree == 4)) {
return &sortedCenterDeg4;
}
if ((dimOfPosition == 3) && (degree == 4))
return &sortedCenterDeg4;
ERROR_EXIT("should not be reached\n");
return NULL;
......@@ -949,7 +941,7 @@ namespace AMDiS {
const int *indi = orderOfPositionIndices(el, posIndex, i);
for (int k = 0; k < nrDOFs; k++)
result[j++] = dof[node0][n0 + indi[k]];
result[j++] = dof[node0][n0 + indi[k]];
}
}
}
......@@ -967,6 +959,36 @@ namespace AMDiS {
getLocalIndices(el, admin, &((*indices)[0]));
}
void Lagrange::getLocalDofPtrVec(const Element *el,
const DOFAdmin *admin,
std::vector<const DegreeOfFreedom*>& vec) const
{
if (static_cast<int>(vec.size()) < nBasFcts)
vec.resize(nBasFcts);
const DegreeOfFreedom **dof = el->getDOF();
GeoIndex posIndex;
int nrDOFs, n0, node0, num = 0;
for (int pos = 0, j = 0; pos <= dim; pos++) {
posIndex = INDEX_OF_DIM(pos, dim);
nrDOFs = admin->getNumberOfDOFs(posIndex);
if (nrDOFs) {
n0 = admin->getNumberOfPreDOFs(posIndex);
node0 = admin->getMesh()->getNode(posIndex);
num = Global::getGeo(posIndex, dim);
for (int i = 0; i < num; node0++, i++) {
const int *indi = orderOfPositionIndices(el, posIndex, i);
for (int k = 0; k < nrDOFs; k++)
vec[j++] = &(dof[node0][n0 + indi[k]]);
}
}
}
}
void Lagrange::l2ScpFctBas(Quadrature *q,
AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
DOFVector<WorldVector<double> >* fh)
......@@ -989,26 +1011,27 @@ namespace AMDiS {
const Parametric *parametric;
TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
TEST_EXIT_DBG(fh->getFESpace())("no fe_space in DOF_REAL_VEC %s\n",
fh->getName().c_str());
TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)("wrong basis fcts for fh\n");
TEST_EXIT_DBG(fh->getFESpace())
("no fe_space in DOF_REAL_VEC %s\n", fh->getName().c_str());
TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)
("wrong basis fcts for fh\n");
Mesh* mesh = fh->getFESpace()->getMesh();
int n_phi = nBasFcts;
if (!q) {
if (!q)
q = Quadrature::provideQuadrature(dim, 2 * degree - 2);
}
const FastQuadrature *quad_fast = FastQuadrature::provideFastQuadrature(this, *q, INIT_PHI);
const FastQuadrature *quad_fast =
FastQuadrature::provideFastQuadrature(this, *q, INIT_PHI);
if ((parametric = mesh->getParametric())) {
if ((parametric = mesh->getParametric()))
ERROR_EXIT("not yet implemented\n");
} else {
else
wdetf_qp = new double[q->getNumPoints()];
}
ElInfo *el_info = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
ElInfo *el_info = stack.traverseFirst(mesh, -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
int numPoints = q->getNumPoints();
......@@ -1024,6 +1047,7 @@ namespace AMDiS {
for (j = 0; j < n_phi; j++) {
for (val = iq = 0; iq < numPoints; iq++)
val += quad_fast->getPhi(iq,j)*wdetf_qp[iq];
(*fh)[dof[j]] += val;
}
......@@ -1455,16 +1479,17 @@ namespace AMDiS {
RCNeighbourList* list,
int n, BasisFunction* basFct)
{
FUNCNAME("Lagrange::refineInter3_3d");
FUNCNAME("Lagrange::refineInter3_3d()");
if (n < 1)
return;
Element *el;
const DegreeOfFreedom *cd;
DegreeOfFreedom pd[20], cdi;
int i, typ, lr_set, node0, n0;
const DOFAdmin *admin;
if (n < 1)
return;
el = list->getElement(0);
typ = list->getType(0);
......@@ -1515,44 +1540,41 @@ namespace AMDiS {
cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
if (typ == 0)
{
/****************************************************************************/
/* parent of el_type 0 */
/****************************************************************************/
if (typ == 0) {
/****************************************************************************/
/* parent of el_type 0 */
/****************************************************************************/
(*drv)[cd[8]] =
(0.0625*(*drv)[pd[0]] + 0.3125*((*drv)[pd[1]] - (*drv)[pd[4]])
+ 0.9375*(*drv)[pd[5]]);
(*drv)[cd[9]] = (*drv)[pd[5]];
(*drv)[cd[17]] =
(0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
+ 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]])
- 0.125*(*drv)[pd[6]] + 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
(*drv)[cd[18]] =
(0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
+ 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]]
+ 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
}
else
{
/****************************************************************************/
/* parent of el_type 0 */
/****************************************************************************/
(*drv)[cd[8]] =
(0.0625*(*drv)[pd[0]] + 0.3125*((*drv)[pd[1]] - (*drv)[pd[4]])
+ 0.9375*(*drv)[pd[5]]);
(*drv)[cd[9]] = (*drv)[pd[5]];
(*drv)[cd[17]] =
(0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
+ 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]]
+ 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
(*drv)[cd[18]] =
(0.0625*((*drv)[pd[0]] - (*drv)[pd[1]]) +
0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[6]]
+ 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
}
(*drv)[cd[8]] =
(0.0625*(*drv)[pd[0]] + 0.3125*((*drv)[pd[1]] - (*drv)[pd[4]])
+ 0.9375*(*drv)[pd[5]]);
(*drv)[cd[9]] = (*drv)[pd[5]];
(*drv)[cd[17]] =
(0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
+ 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]])
- 0.125*(*drv)[pd[6]] + 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
(*drv)[cd[18]] =
(0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
+ 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]]
+ 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
} else {
/****************************************************************************/
/* parent of el_type 0 */
/****************************************************************************/
(*drv)[cd[8]] =
(0.0625*(*drv)[pd[0]] + 0.3125*((*drv)[pd[1]] - (*drv)[pd[4]])
+ 0.9375*(*drv)[pd[5]]);
(*drv)[cd[9]] = (*drv)[pd[5]];
(*drv)[cd[17]] =
(0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
+ 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]]
+ 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
(*drv)[cd[18]] =
(0.0625*((*drv)[pd[0]] - (*drv)[pd[1]]) +
0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[6]]
+ 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
}
/****************************************************************************/
/* adjust neighbour values */
......@@ -1561,125 +1583,121 @@ namespace AMDiS {
node0 = drv->getFESpace()->getMesh()->getNode(FACE);
n0 = admin->getNumberOfPreDOFs(FACE);
for (i = 1; i < n; i++)
{
el = list->getElement(i);
typ = list->getType(i);
basFct->getLocalIndices(el, admin, pd);
lr_set = 0;
if (list->getNeighbourElement(i, 0) && list->getNeighbourNr(i, 0) < i)
lr_set = 1;
if (list->getNeighbourElement(i, 1) && list->getNeighbourNr(i, 1) < i)
lr_set += 2;
TEST_EXIT_DBG(lr_set)("no values set on both neighbours\n");
/****************************************************************************/
/* values on child[0] */