diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h index 9b55d279ebc60f826368b8d8f4226ec5298c4122..6a746378872a9a83af10abb8dc0efe664772a405 100644 --- a/AMDiS/src/BasisFunction.h +++ b/AMDiS/src/BasisFunction.h @@ -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 {} diff --git a/AMDiS/src/ElementDofIterator.cc b/AMDiS/src/ElementDofIterator.cc index 78c1f3b662c2b8cbcc0f74782c73e46aab3d5a24..d8dc65cacf34a95aa17891da03bd15abc6ff9316 100644 --- a/AMDiS/src/ElementDofIterator.cc +++ b/AMDiS/src/ElementDofIterator.cc @@ -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(); + } } diff --git a/AMDiS/src/ElementDofIterator.h b/AMDiS/src/ElementDofIterator.h index 7b571f79cbcf11416806d3cd92275604160f1759..a09e1f305089d4d9c9d78735e69feadd204a0f4f 100644 --- a/AMDiS/src/ElementDofIterator.h +++ b/AMDiS/src/ElementDofIterator.h @@ -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; diff --git a/AMDiS/src/InteriorBoundary.h b/AMDiS/src/InteriorBoundary.h index 59f5065cb4317ded13b8324f0eab347c3604738a..077d6d40ea20933c9506582cd2fcdb734e8ec2bb 100644 --- a/AMDiS/src/InteriorBoundary.h +++ b/AMDiS/src/InteriorBoundary.h @@ -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; }; } diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc index 5a32a601bacc199f578fb6b94eab287c2574db07..034a17783f4a537bdd521d83159dddf37c0653a3 100644 --- a/AMDiS/src/Lagrange.cc +++ b/AMDiS/src/Lagrange.cc @@ -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] */ - /****************************************************************************/ - - cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); - - switch(lr_set) - { - case 1: - (*drv)[cd[12]] = - (0.0625*((*drv)[pd[0]]+(*drv)[pd[1]]-(*drv)[pd[4]]-(*drv)[pd[5]]) - + 0.25*(-(*drv)[pd[6]] - (*drv)[pd[10]]) - + 0.5*((*drv)[pd[7]] + (*drv)[pd[11]] + (*drv)[pd[19]])); - (*drv)[cd[13]] = (*drv)[pd[19]]; - (*drv)[cd[16]] = - (0.0625*((*drv)[pd[0]]+(*drv)[pd[1]]-(*drv)[pd[4]]-(*drv)[pd[5]]) - + 0.125*(-(*drv)[pd[6]]-(*drv)[pd[8]]-(*drv)[pd[10]]-(*drv)[pd[12]]) - + 0.5*((*drv)[pd[16]] + (*drv)[pd[17]]) - + 0.25*((*drv)[pd[18]] + (*drv)[pd[19]])); - (*drv)[cd[18]] = - (0.0625*(-(*drv)[pd[0]] + (*drv)[pd[1]]) - + 0.1875*((*drv)[pd[4]] - (*drv)[pd[5]]) + 0.375*(*drv)[pd[6]] - - 0.125*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]); - break; - case 2: - (*drv)[cd[14]] = - (0.0625*((*drv)[pd[0]]+(*drv)[pd[1]]-(*drv)[pd[4]]-(*drv)[pd[5]]) - + 0.25*(-(*drv)[pd[8]] - (*drv)[pd[12]]) - + 0.5*((*drv)[pd[9]] + (*drv)[pd[13]] + (*drv)[pd[18]])); - (*drv)[cd[15]] = (*drv)[pd[18]]; - (*drv)[cd[16]] = - (0.0625*((*drv)[pd[0]]+(*drv)[pd[1]]-(*drv)[pd[4]]-(*drv)[pd[5]]) - + 0.125*(-(*drv)[pd[6]]-(*drv)[pd[8]]-(*drv)[pd[10]]-(*drv)[pd[12]]) - + 0.5*((*drv)[pd[16]] + (*drv)[pd[17]]) - + 0.25*((*drv)[pd[18]] + (*drv)[pd[19]])); - (*drv)[cd[17]] = - (0.0625*(-(*drv)[pd[0]] + (*drv)[pd[1]]) - + 0.1875*((*drv)[pd[4]] - (*drv)[pd[5]]) + 0.375*(*drv)[pd[8]] - - 0.125*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]); - break; - case 3: - (*drv)[cd[16]] = - (0.0625*((*drv)[pd[0]]+(*drv)[pd[1]]-(*drv)[pd[4]]-(*drv)[pd[5]]) - + 0.125*(-(*drv)[pd[6]]-(*drv)[pd[8]]-(*drv)[pd[10]]-(*drv)[pd[12]]) - + 0.5*((*drv)[pd[16]] + (*drv)[pd[17]]) - + 0.25*((*drv)[pd[18]] + (*drv)[pd[19]])); - break; - } + for (i = 1; i < n; i++) { + el = list->getElement(i); + typ = list->getType(i); + basFct->getLocalIndices(el, admin, pd); - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ - - cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); - - if (typ == 0) { - switch(lr_set) - { - case 1: - cdi = el->getChild(1)->getDOF(node0+1, n0); - TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n"); - (*drv)[cdi] = - (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]]); - break; - case 2: - cdi = el->getChild(1)->getDOF(node0+2, n0); - TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n"); - (*drv)[cdi] = - (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]]); - break; - } - } else { - switch(lr_set) - { - case 1: - cdi = el->getChild(1)->getDOF(node0+2, n0); - TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n"); - (*drv)[cdi] = - (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]]); - break; - case 2: - cdi = el->getChild(1)->getDOF(node0+1, n0); - TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n"); - (*drv)[cdi] = - (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]]); - break; - } + 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] */ + /****************************************************************************/ + + cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); + + switch (lr_set) { + case 1: + (*drv)[cd[12]] = + (0.0625*((*drv)[pd[0]]+(*drv)[pd[1]]-(*drv)[pd[4]]-(*drv)[pd[5]]) + + 0.25*(-(*drv)[pd[6]] - (*drv)[pd[10]]) + + 0.5*((*drv)[pd[7]] + (*drv)[pd[11]] + (*drv)[pd[19]])); + (*drv)[cd[13]] = (*drv)[pd[19]]; + (*drv)[cd[16]] = + (0.0625*((*drv)[pd[0]]+(*drv)[pd[1]]-(*drv)[pd[4]]-(*drv)[pd[5]]) + + 0.125*(-(*drv)[pd[6]]-(*drv)[pd[8]]-(*drv)[pd[10]]-(*drv)[pd[12]]) + + 0.5*((*drv)[pd[16]] + (*drv)[pd[17]]) + + 0.25*((*drv)[pd[18]] + (*drv)[pd[19]])); + (*drv)[cd[18]] = + (0.0625*(-(*drv)[pd[0]] + (*drv)[pd[1]]) + + 0.1875*((*drv)[pd[4]] - (*drv)[pd[5]]) + 0.375*(*drv)[pd[6]] + - 0.125*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]); + break; + case 2: + (*drv)[cd[14]] = + (0.0625*((*drv)[pd[0]]+(*drv)[pd[1]]-(*drv)[pd[4]]-(*drv)[pd[5]]) + + 0.25*(-(*drv)[pd[8]] - (*drv)[pd[12]]) + + 0.5*((*drv)[pd[9]] + (*drv)[pd[13]] + (*drv)[pd[18]])); + (*drv)[cd[15]] = (*drv)[pd[18]]; + (*drv)[cd[16]] = + (0.0625*((*drv)[pd[0]]+(*drv)[pd[1]]-(*drv)[pd[4]]-(*drv)[pd[5]]) + + 0.125*(-(*drv)[pd[6]]-(*drv)[pd[8]]-(*drv)[pd[10]]-(*drv)[pd[12]]) + + 0.5*((*drv)[pd[16]] + (*drv)[pd[17]]) + + 0.25*((*drv)[pd[18]] + (*drv)[pd[19]])); + (*drv)[cd[17]] = + (0.0625*(-(*drv)[pd[0]] + (*drv)[pd[1]]) + + 0.1875*((*drv)[pd[4]] - (*drv)[pd[5]]) + 0.375*(*drv)[pd[8]] + - 0.125*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]); + break; + case 3: + (*drv)[cd[16]] = + (0.0625*((*drv)[pd[0]]+(*drv)[pd[1]]-(*drv)[pd[4]]-(*drv)[pd[5]]) + + 0.125*(-(*drv)[pd[6]]-(*drv)[pd[8]]-(*drv)[pd[10]]-(*drv)[pd[12]]) + + 0.5*((*drv)[pd[16]] + (*drv)[pd[17]]) + + 0.25*((*drv)[pd[18]] + (*drv)[pd[19]])); + break; + } + + /****************************************************************************/ + /* values on child[1] */ + /****************************************************************************/ + + cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); + + if (typ == 0) { + switch (lr_set) { + case 1: + cdi = el->getChild(1)->getDOF(node0+1, n0); + TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n"); + (*drv)[cdi] = + (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]]); + break; + case 2: + cdi = el->getChild(1)->getDOF(node0+2, n0); + TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n"); + (*drv)[cdi] = + (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]]); + break; + } + } else { + switch (lr_set) { + case 1: + cdi = el->getChild(1)->getDOF(node0+2, n0); + TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n"); + (*drv)[cdi] = + (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]]); + break; + case 2: + cdi = el->getChild(1)->getDOF(node0+1, n0); + TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n"); + (*drv)[cdi] = + (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]]); + break; } } + } } void Lagrange::refineInter4_1d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::refineInter4_1d"); + FUNCNAME("Lagrange::refineInter4_1d()"); if (n < 1) return; @@ -1695,7 +1713,7 @@ namespace AMDiS { /****************************************************************************/ const DegreeOfFreedom *cdof = - basFct->getLocalIndices(el->getChild(0), admin, NULL); + basFct->getLocalIndices(el->getChild(0), admin, NULL); (*drv)[cdof[2]] = (*drv)[pdof[10]]; (*drv)[cdof[3]] = @@ -2096,447 +2114,440 @@ namespace AMDiS { /****************************************************************************/ cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); - if (typ == 0) - { - /****************************************************************************/ - /* parent of el_type 0 */ - /****************************************************************************/ - - (*drv)[cd[10]] = - (-0.0390625*(*drv)[pd[0]] + 0.2734375*(*drv)[pd[1]] - + 0.21875*(*drv)[pd[4]] - 0.546875*(*drv)[pd[5]] - + 1.09375*(*drv)[pd[6]]); - (*drv)[cd[11]] = (*drv)[pd[6]]; - (*drv)[cd[12]] = - (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] - - 0.15625*(*drv)[pd[4]] + 0.703125*(*drv)[pd[5]] - + 0.46875*(*drv)[pd[6]]); - (*drv)[cd[25]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) - - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[7]] - + 0.3125*((*drv)[pd[13]] - (*drv)[pd[31]]) - + 0.9375*(*drv)[pd[32]]); - (*drv)[cd[26]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] - (*drv)[pd[13]]) - + 0.375*((*drv)[pd[14]] - (*drv)[pd[31]] + (*drv)[pd[32]]) + if (typ == 0) { + /****************************************************************************/ + /* parent of el_type 0 */ + /****************************************************************************/ + + (*drv)[cd[10]] = + (-0.0390625*(*drv)[pd[0]] + 0.2734375*(*drv)[pd[1]] + + 0.21875*(*drv)[pd[4]] - 0.546875*(*drv)[pd[5]] + + 1.09375*(*drv)[pd[6]]); + (*drv)[cd[11]] = (*drv)[pd[6]]; + (*drv)[cd[12]] = + (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] + - 0.15625*(*drv)[pd[4]] + 0.703125*(*drv)[pd[5]] + + 0.46875*(*drv)[pd[6]]); + (*drv)[cd[25]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) + - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[7]] + + 0.3125*((*drv)[pd[13]] - (*drv)[pd[31]]) + + 0.9375*(*drv)[pd[32]]); + (*drv)[cd[26]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] - (*drv)[pd[13]]) + + 0.375*((*drv)[pd[14]] - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[33]]); + (*drv)[cd[27]] = (*drv)[pd[32]]; + (*drv)[cd[28]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) + - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[10]] + + 0.3125*((*drv)[pd[16]] - (*drv)[pd[28]]) + + 0.9375*(*drv)[pd[29]]); + (*drv)[cd[29]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] - (*drv)[pd[16]]) + + 0.375*((*drv)[pd[17]] - (*drv)[pd[28]] + (*drv)[pd[29]]) + + 0.75*(*drv)[pd[30]]); + (*drv)[cd[30]] = (*drv)[pd[29]]; + (*drv)[cd[34]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] + -(*drv)[pd[13]] - (*drv)[pd[16]]) + + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] + + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] + -(*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[34]]); + } else { + /****************************************************************************/ + /* parent of el_type 1|2 */ + /****************************************************************************/ + + (*drv)[cd[10]] = + (-0.0390625*(*drv)[pd[0]] + 0.2734375*(*drv)[pd[1]] + + 0.21875*(*drv)[pd[4]] - 0.546875*(*drv)[pd[5]] + + 1.09375*(*drv)[pd[6]]); + (*drv)[cd[11]] = (*drv)[pd[6]]; + (*drv)[cd[12]] = + (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] + - 0.15625*(*drv)[pd[4]] + 0.703125*(*drv)[pd[5]] + + 0.46875*(*drv)[pd[6]]); + (*drv)[cd[25]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) + - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[10]] + + 0.3125*((*drv)[pd[16]] - (*drv)[pd[28]]) + + 0.9375*(*drv)[pd[29]]); + (*drv)[cd[26]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] - (*drv)[pd[16]]) + + 0.375*((*drv)[pd[17]] - (*drv)[pd[28]] + (*drv)[pd[29]]) + + 0.75*(*drv)[pd[30]]); + (*drv)[cd[27]] = (*drv)[pd[29]]; + (*drv)[cd[28]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) + - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[7]] + + 0.3125*((*drv)[pd[13]] - (*drv)[pd[31]]) + + 0.9375*(*drv)[pd[32]]); + (*drv)[cd[29]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] - (*drv)[pd[13]]) + + 0.375*((*drv)[pd[14]] - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[33]]); + (*drv)[cd[30]] = (*drv)[pd[32]]; + (*drv)[cd[34]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] + - (*drv)[pd[13]] - (*drv)[pd[16]]) + + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] + + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] + - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[34]]); + } + + /****************************************************************************/ + /* adjust neighbour values */ + /****************************************************************************/ + + 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] */ + /****************************************************************************/ + + cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); + + switch (lr_set) { + case 1: + (*drv)[cd[16]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) + + 0.015625*(*drv)[pd[5]] + + 0.1875*((*drv)[pd[7]] + (*drv)[pd[13]] + - (*drv)[pd[31]] - (*drv)[pd[32]]) + + 0.375*(-(*drv)[pd[8]] - (*drv)[pd[14]]) + + 0.5*((*drv)[pd[9]] + (*drv)[pd[15]]) + 0.75*(*drv)[pd[33]]); - (*drv)[cd[27]] = (*drv)[pd[32]]; + (*drv)[cd[17]] = (*drv)[pd[33]]; + (*drv)[cd[18]] = + (0.0234375*((*drv)[pd[0]] + (*drv)[pd[1]]) + + 0.09375*(-(*drv)[pd[4]] - (*drv)[pd[6]]) + + 0.140625*(*drv)[pd[5]] + + 0.0625*(-(*drv)[pd[7]] - (*drv)[pd[13]]) + + 0.5625*((*drv)[pd[31]] + (*drv)[pd[32]])); + (*drv)[cd[22]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) + + 0.015625*(*drv)[pd[5]] + + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] + (*drv)[pd[13]] + - (*drv)[pd[14]] - (*drv)[pd[31]] - (*drv)[pd[32]]) + + 0.0625*((*drv)[pd[10]] + (*drv)[pd[16]] + - (*drv)[pd[28]] - (*drv)[pd[29]]) + + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[33]]) + + 0.5*((*drv)[pd[23]] + (*drv)[pd[26]] + (*drv)[pd[34]])); + (*drv)[cd[23]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) + + 0.015625*(*drv)[pd[5]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[13]] + - (*drv)[pd[31]] - (*drv)[pd[32]]) + + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] + (*drv)[pd[16]] + - (*drv)[pd[17]] - (*drv)[pd[28]] - (*drv)[pd[29]]) + + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[30]]) + + 0.5*((*drv)[pd[24]] + (*drv)[pd[27]] + (*drv)[pd[34]])); + (*drv)[cd[24]] = (*drv)[pd[34]]; (*drv)[cd[28]] = (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) - - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[10]] - + 0.3125*((*drv)[pd[16]] - (*drv)[pd[28]]) - + 0.9375*(*drv)[pd[29]]); + - 0.234375*(*drv)[pd[5]] + + 0.3125*((*drv)[pd[7]] - (*drv)[pd[32]]) + + 0.0625*(*drv)[pd[13]] + 0.9375*(*drv)[pd[31]]); (*drv)[cd[29]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] - (*drv)[pd[16]]) - + 0.375*((*drv)[pd[17]] - (*drv)[pd[28]] + (*drv)[pd[29]]) - + 0.75*(*drv)[pd[30]]); - (*drv)[cd[30]] = (*drv)[pd[29]]; + (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] + - 0.03125*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + + 0.09375*(*drv)[pd[6]] + + 0.125*(-(*drv)[pd[7]] + (*drv)[pd[13]] - (*drv)[pd[14]]) + + 0.375*((*drv)[pd[8]] + (*drv)[pd[31]] - (*drv)[pd[32]]) + + 0.75*(*drv)[pd[33]]); + (*drv)[cd[30]] = (*drv)[pd[31]]; (*drv)[cd[34]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] - -(*drv)[pd[13]] - (*drv)[pd[16]]) - + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] - + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] - -(*drv)[pd[31]] + (*drv)[pd[32]]) - + 0.75*(*drv)[pd[34]]); - } - else - { - /****************************************************************************/ - /* parent of el_type 1|2 */ - /****************************************************************************/ - - (*drv)[cd[10]] = - (-0.0390625*(*drv)[pd[0]] + 0.2734375*(*drv)[pd[1]] - + 0.21875*(*drv)[pd[4]] - 0.546875*(*drv)[pd[5]] - + 1.09375*(*drv)[pd[6]]); - (*drv)[cd[11]] = (*drv)[pd[6]]; - (*drv)[cd[12]] = (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] - - 0.15625*(*drv)[pd[4]] + 0.703125*(*drv)[pd[5]] - + 0.46875*(*drv)[pd[6]]); + - 0.03125*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + + 0.09375*(*drv)[pd[6]] + + 0.0625*(-(*drv)[pd[7]] - (*drv)[pd[10]] + + (*drv)[pd[13]] + (*drv)[pd[16]]) + - 0.125*(*drv)[pd[22]] + 0.375*(*drv)[pd[25]] + + 0.1875*((*drv)[pd[28]] - (*drv)[pd[29]] + + (*drv)[pd[31]] - (*drv)[pd[32]]) + + 0.75*(*drv)[pd[34]]); + break; + case 2: + (*drv)[cd[19]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) + + 0.015625*(*drv)[pd[5]] + + 0.1875*((*drv)[pd[10]] + (*drv)[pd[16]] + - (*drv)[pd[28]] - (*drv)[pd[29]]) + + 0.375*(-(*drv)[pd[11]] - (*drv)[pd[17]]) + + 0.5*((*drv)[pd[12]] + (*drv)[pd[18]]) + + 0.75*(*drv)[pd[30]]); + (*drv)[cd[20]] = (*drv)[pd[30]]; + (*drv)[cd[21]] = + (0.0234375*((*drv)[pd[0]] + (*drv)[pd[1]]) + + 0.09375*(-(*drv)[pd[4]] - (*drv)[pd[6]]) + + 0.140625*(*drv)[pd[5]] + + 0.0625*(-(*drv)[pd[10]] - (*drv)[pd[16]]) + + 0.5625*((*drv)[pd[28]] + (*drv)[pd[29]])); + (*drv)[cd[22]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) + + 0.015625*(*drv)[pd[5]] + + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] + (*drv)[pd[13]] + - (*drv)[pd[14]] - (*drv)[pd[31]] - (*drv)[pd[32]]) + + 0.0625*((*drv)[pd[10]] + (*drv)[pd[16]] + - (*drv)[pd[28]] - (*drv)[pd[29]]) + + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[33]]) + + 0.5*((*drv)[pd[23]] + (*drv)[pd[26]] + (*drv)[pd[34]])); + (*drv)[cd[23]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) + + 0.015625*(*drv)[pd[5]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[13]] + - (*drv)[pd[31]] - (*drv)[pd[32]]) + + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] + (*drv)[pd[16]] + - (*drv)[pd[17]] - (*drv)[pd[28]] - (*drv)[pd[29]]) + + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[30]]) + + 0.5*((*drv)[pd[24]] + (*drv)[pd[27]] + (*drv)[pd[34]])); + (*drv)[cd[24]] = (*drv)[pd[34]]; (*drv)[cd[25]] = (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) - - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[10]] - + 0.3125*((*drv)[pd[16]] - (*drv)[pd[28]]) - + 0.9375*(*drv)[pd[29]]); + - 0.234375*(*drv)[pd[5]] + + 0.3125*((*drv)[pd[10]] - (*drv)[pd[29]]) + + 0.0625*(*drv)[pd[16]] + 0.9375*(*drv)[pd[28]]); (*drv)[cd[26]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] - (*drv)[pd[16]]) - + 0.375*((*drv)[pd[17]] - (*drv)[pd[28]] + (*drv)[pd[29]]) + (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] + - 0.03125*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + + 0.09375*(*drv)[pd[6]] + + 0.125*(-(*drv)[pd[10]] + (*drv)[pd[16]] - (*drv)[pd[17]]) + + 0.375*((*drv)[pd[11]] + (*drv)[pd[28]] - (*drv)[pd[29]]) + 0.75*(*drv)[pd[30]]); - (*drv)[cd[27]] = (*drv)[pd[29]]; - (*drv)[cd[28]] = + (*drv)[cd[27]] = (*drv)[pd[28]]; + (*drv)[cd[34]] = + (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] + - 0.03125*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + + 0.09375*(*drv)[pd[6]] + + 0.0625*(-(*drv)[pd[7]] - (*drv)[pd[10]] + + (*drv)[pd[13]] + (*drv)[pd[16]]) + - 0.125*(*drv)[pd[22]] + 0.375*(*drv)[pd[25]] + + 0.1875*((*drv)[pd[28]] - (*drv)[pd[29]] + + (*drv)[pd[31]] - (*drv)[pd[32]]) + + 0.75*(*drv)[pd[34]]); + break; + case 3: + (*drv)[cd[22]] = (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) - - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[7]] - + 0.3125*((*drv)[pd[13]] - (*drv)[pd[31]]) - + 0.9375*(*drv)[pd[32]]); - (*drv)[cd[29]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] - (*drv)[pd[13]]) - + 0.375*((*drv)[pd[14]] - (*drv)[pd[31]] + (*drv)[pd[32]]) - + 0.75*(*drv)[pd[33]]); - (*drv)[cd[30]] = (*drv)[pd[32]]; + + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) + + 0.015625*(*drv)[pd[5]] + + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] + (*drv)[pd[13]] + - (*drv)[pd[14]] - (*drv)[pd[31]] - (*drv)[pd[32]]) + + 0.0625*((*drv)[pd[10]] + (*drv)[pd[16]] + - (*drv)[pd[28]] - (*drv)[pd[29]]) + + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[33]]) + + 0.5*((*drv)[pd[23]] + (*drv)[pd[26]] + (*drv)[pd[34]])); + (*drv)[cd[23]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) + + 0.015625*(*drv)[pd[5]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[13]] + - (*drv)[pd[31]] - (*drv)[pd[32]]) + + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] + (*drv)[pd[16]] + - (*drv)[pd[17]] - (*drv)[pd[28]] - (*drv)[pd[29]]) + + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[30]]) + + 0.5*((*drv)[pd[24]] + (*drv)[pd[27]] + (*drv)[pd[34]])); + (*drv)[cd[24]] = (*drv)[pd[34]]; (*drv)[cd[34]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] - - (*drv)[pd[13]] - (*drv)[pd[16]]) - + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] - + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] - - (*drv)[pd[31]] + (*drv)[pd[32]]) + (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] + - 0.03125*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + + 0.09375*(*drv)[pd[6]] + + 0.0625*(-(*drv)[pd[7]] - (*drv)[pd[10]] + + (*drv)[pd[13]] + (*drv)[pd[16]]) + - 0.125*(*drv)[pd[22]] + 0.375*(*drv)[pd[25]] + + 0.1875*((*drv)[pd[28]] - (*drv)[pd[29]] + + (*drv)[pd[31]] - (*drv)[pd[32]]) + 0.75*(*drv)[pd[34]]); + break; } - /****************************************************************************/ - /* adjust neighbour values */ - /****************************************************************************/ - - 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] */ - /****************************************************************************/ - - cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); - - switch(lr_set) - { - case 1: - (*drv)[cd[16]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) - + 0.015625*(*drv)[pd[5]] - + 0.1875*((*drv)[pd[7]] + (*drv)[pd[13]] - - (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.375*(-(*drv)[pd[8]] - (*drv)[pd[14]]) - + 0.5*((*drv)[pd[9]] + (*drv)[pd[15]]) - + 0.75*(*drv)[pd[33]]); - (*drv)[cd[17]] = (*drv)[pd[33]]; - (*drv)[cd[18]] = - (0.0234375*((*drv)[pd[0]] + (*drv)[pd[1]]) - + 0.09375*(-(*drv)[pd[4]] - (*drv)[pd[6]]) - + 0.140625*(*drv)[pd[5]] - + 0.0625*(-(*drv)[pd[7]] - (*drv)[pd[13]]) - + 0.5625*((*drv)[pd[31]] + (*drv)[pd[32]])); - (*drv)[cd[22]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) - + 0.015625*(*drv)[pd[5]] - + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] + (*drv)[pd[13]] - - (*drv)[pd[14]] - (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.0625*((*drv)[pd[10]] + (*drv)[pd[16]] - - (*drv)[pd[28]] - (*drv)[pd[29]]) - + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[33]]) - + 0.5*((*drv)[pd[23]] + (*drv)[pd[26]] + (*drv)[pd[34]])); - (*drv)[cd[23]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) - + 0.015625*(*drv)[pd[5]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[13]] - - (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] + (*drv)[pd[16]] - - (*drv)[pd[17]] - (*drv)[pd[28]] - (*drv)[pd[29]]) - + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[30]]) - + 0.5*((*drv)[pd[24]] + (*drv)[pd[27]] + (*drv)[pd[34]])); - (*drv)[cd[24]] = (*drv)[pd[34]]; - (*drv)[cd[28]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) - - 0.234375*(*drv)[pd[5]] - + 0.3125*((*drv)[pd[7]] - (*drv)[pd[32]]) - + 0.0625*(*drv)[pd[13]] + 0.9375*(*drv)[pd[31]]); - (*drv)[cd[29]] = - (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] - - 0.03125*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - + 0.09375*(*drv)[pd[6]] - + 0.125*(-(*drv)[pd[7]] + (*drv)[pd[13]] - (*drv)[pd[14]]) - + 0.375*((*drv)[pd[8]] + (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.75*(*drv)[pd[33]]); - (*drv)[cd[30]] = (*drv)[pd[31]]; - (*drv)[cd[34]] = - (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] - - 0.03125*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - + 0.09375*(*drv)[pd[6]] - + 0.0625*(-(*drv)[pd[7]] - (*drv)[pd[10]] - + (*drv)[pd[13]] + (*drv)[pd[16]]) - - 0.125*(*drv)[pd[22]] + 0.375*(*drv)[pd[25]] - + 0.1875*((*drv)[pd[28]] - (*drv)[pd[29]] - + (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.75*(*drv)[pd[34]]); - break; - case 2: - (*drv)[cd[19]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) - + 0.015625*(*drv)[pd[5]] - + 0.1875*((*drv)[pd[10]] + (*drv)[pd[16]] - - (*drv)[pd[28]] - (*drv)[pd[29]]) - + 0.375*(-(*drv)[pd[11]] - (*drv)[pd[17]]) - + 0.5*((*drv)[pd[12]] + (*drv)[pd[18]]) - + 0.75*(*drv)[pd[30]]); - (*drv)[cd[20]] = (*drv)[pd[30]]; - (*drv)[cd[21]] = - (0.0234375*((*drv)[pd[0]] + (*drv)[pd[1]]) - + 0.09375*(-(*drv)[pd[4]] - (*drv)[pd[6]]) - + 0.140625*(*drv)[pd[5]] - + 0.0625*(-(*drv)[pd[10]] - (*drv)[pd[16]]) - + 0.5625*((*drv)[pd[28]] + (*drv)[pd[29]])); - (*drv)[cd[22]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) - + 0.015625*(*drv)[pd[5]] - + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] + (*drv)[pd[13]] - - (*drv)[pd[14]] - (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.0625*((*drv)[pd[10]] + (*drv)[pd[16]] - - (*drv)[pd[28]] - (*drv)[pd[29]]) - + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[33]]) - + 0.5*((*drv)[pd[23]] + (*drv)[pd[26]] + (*drv)[pd[34]])); - (*drv)[cd[23]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) - + 0.015625*(*drv)[pd[5]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[13]] - - (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] + (*drv)[pd[16]] - - (*drv)[pd[17]] - (*drv)[pd[28]] - (*drv)[pd[29]]) - + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[30]]) - + 0.5*((*drv)[pd[24]] + (*drv)[pd[27]] + (*drv)[pd[34]])); - (*drv)[cd[24]] = (*drv)[pd[34]]; - (*drv)[cd[25]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) - - 0.234375*(*drv)[pd[5]] - + 0.3125*((*drv)[pd[10]] - (*drv)[pd[29]]) - + 0.0625*(*drv)[pd[16]] + 0.9375*(*drv)[pd[28]]); - (*drv)[cd[26]] = - (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] - - 0.03125*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - + 0.09375*(*drv)[pd[6]] - + 0.125*(-(*drv)[pd[10]] + (*drv)[pd[16]] - (*drv)[pd[17]]) - + 0.375*((*drv)[pd[11]] + (*drv)[pd[28]] - (*drv)[pd[29]]) - + 0.75*(*drv)[pd[30]]); - (*drv)[cd[27]] = (*drv)[pd[28]]; - (*drv)[cd[34]] = - (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] - - 0.03125*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - + 0.09375*(*drv)[pd[6]] - + 0.0625*(-(*drv)[pd[7]] - (*drv)[pd[10]] - + (*drv)[pd[13]] + (*drv)[pd[16]]) - - 0.125*(*drv)[pd[22]] + 0.375*(*drv)[pd[25]] - + 0.1875*((*drv)[pd[28]] - (*drv)[pd[29]] - + (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.75*(*drv)[pd[34]]); - break; - case 3: - (*drv)[cd[22]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) - + 0.015625*(*drv)[pd[5]] - + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] + (*drv)[pd[13]] - - (*drv)[pd[14]] - (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.0625*((*drv)[pd[10]] + (*drv)[pd[16]] - - (*drv)[pd[28]] - (*drv)[pd[29]]) - + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[33]]) - + 0.5*((*drv)[pd[23]] + (*drv)[pd[26]] + (*drv)[pd[34]])); - (*drv)[cd[23]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.03125*((*drv)[pd[4]] + (*drv)[pd[6]]) - + 0.015625*(*drv)[pd[5]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[13]] - - (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] + (*drv)[pd[16]] - - (*drv)[pd[17]] - (*drv)[pd[28]] - (*drv)[pd[29]]) - + 0.25*(-(*drv)[pd[22]] - (*drv)[pd[25]] + (*drv)[pd[30]]) - + 0.5*((*drv)[pd[24]] + (*drv)[pd[27]] + (*drv)[pd[34]])); - (*drv)[cd[24]] = (*drv)[pd[34]]; - (*drv)[cd[34]] = - (0.0234375*(*drv)[pd[0]] - 0.0390625*(*drv)[pd[1]] - - 0.03125*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - + 0.09375*(*drv)[pd[6]] - + 0.0625*(-(*drv)[pd[7]] - (*drv)[pd[10]] - + (*drv)[pd[13]] + (*drv)[pd[16]]) - - 0.125*(*drv)[pd[22]] + 0.375*(*drv)[pd[25]] - + 0.1875*((*drv)[pd[28]] - (*drv)[pd[29]] - + (*drv)[pd[31]] - (*drv)[pd[32]]) - + 0.75*(*drv)[pd[34]]); - break; - } + /****************************************************************************/ + /* values on child[1] */ + /****************************************************************************/ - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ - - cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); - - if (typ == 0) - { - switch(lr_set) - { - case 1: - (*drv)[cd[25]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) - - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[7]] - + 0.3125*((*drv)[pd[13]] - (*drv)[pd[31]]) - + 0.9375*(*drv)[pd[32]]); - (*drv)[cd[26]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] - (*drv)[pd[13]]) - + 0.375*((*drv)[pd[14]] - (*drv)[pd[31]] + (*drv)[pd[32]]) - + 0.75*(*drv)[pd[33]]); - (*drv)[cd[27]] = (*drv)[pd[32]]; - (*drv)[cd[34]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] - - (*drv)[pd[13]] - (*drv)[pd[16]]) - + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] - + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] - - (*drv)[pd[31]] + (*drv)[pd[32]]) - + 0.75*(*drv)[pd[34]]); - break; - case 2: - (*drv)[cd[28]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) - - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[10]] - + 0.3125*((*drv)[pd[16]] - (*drv)[pd[28]]) - + 0.9375*(*drv)[pd[29]]); - (*drv)[cd[29]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] - (*drv)[pd[16]]) - + 0.375*((*drv)[pd[17]] - (*drv)[pd[28]] + (*drv)[pd[29]]) - + 0.75*(*drv)[pd[30]]); - (*drv)[cd[30]] = (*drv)[pd[29]]; - (*drv)[cd[34]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] - - (*drv)[pd[13]] - (*drv)[pd[16]]) - + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] - + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] - - (*drv)[pd[31]] + (*drv)[pd[32]]) - + 0.75*(*drv)[pd[34]]); - break; - case 3: - (*drv)[cd[34]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] - - (*drv)[pd[13]] - (*drv)[pd[16]]) - + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] - + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] - - (*drv)[pd[31]] + (*drv)[pd[32]]) - + 0.75*(*drv)[pd[34]]); - break; - } - } else { - switch(lr_set) { - case 1: - (*drv)[cd[28]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) - - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[7]] - + 0.3125*((*drv)[pd[13]] - (*drv)[pd[31]]) - + 0.9375*(*drv)[pd[32]]); - (*drv)[cd[29]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] - (*drv)[pd[13]]) - + 0.375*((*drv)[pd[14]] - (*drv)[pd[31]] + (*drv)[pd[32]]) - + 0.75*(*drv)[pd[33]]); - (*drv)[cd[30]] = (*drv)[pd[32]]; - (*drv)[cd[34]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] - - (*drv)[pd[13]] - (*drv)[pd[16]]) - + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] - + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] - - (*drv)[pd[31]] + (*drv)[pd[32]]) - + 0.75*(*drv)[pd[34]]); - break; - case 2: - (*drv)[cd[25]] = - (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) - + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) - - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[10]] - + 0.3125*((*drv)[pd[16]] - (*drv)[pd[28]]) - + 0.9375*(*drv)[pd[29]]); - (*drv)[cd[26]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] - (*drv)[pd[16]]) - + 0.375*((*drv)[pd[17]] - (*drv)[pd[28]] + (*drv)[pd[29]]) - + 0.75*(*drv)[pd[30]]); - (*drv)[cd[27]] = (*drv)[pd[29]]; - (*drv)[cd[34]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] - - (*drv)[pd[13]] - (*drv)[pd[16]]) - + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] - + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] - - (*drv)[pd[31]] + (*drv)[pd[32]]) - + 0.75*(*drv)[pd[34]]); - break; - case 3: - (*drv)[cd[34]] = - (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] - + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] - - 0.03125*(*drv)[pd[6]] - + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] - - (*drv)[pd[13]] - (*drv)[pd[16]]) - + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] - + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] - - (*drv)[pd[31]] + (*drv)[pd[32]]) - + 0.75*(*drv)[pd[34]]); - break; - } + cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); + + if (typ == 0) { + switch (lr_set) { + case 1: + (*drv)[cd[25]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) + - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[7]] + + 0.3125*((*drv)[pd[13]] - (*drv)[pd[31]]) + + 0.9375*(*drv)[pd[32]]); + (*drv)[cd[26]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] - (*drv)[pd[13]]) + + 0.375*((*drv)[pd[14]] - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[33]]); + (*drv)[cd[27]] = (*drv)[pd[32]]; + (*drv)[cd[34]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] + - (*drv)[pd[13]] - (*drv)[pd[16]]) + + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] + + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] + - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[34]]); + break; + case 2: + (*drv)[cd[28]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) + - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[10]] + + 0.3125*((*drv)[pd[16]] - (*drv)[pd[28]]) + + 0.9375*(*drv)[pd[29]]); + (*drv)[cd[29]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] - (*drv)[pd[16]]) + + 0.375*((*drv)[pd[17]] - (*drv)[pd[28]] + (*drv)[pd[29]]) + + 0.75*(*drv)[pd[30]]); + (*drv)[cd[30]] = (*drv)[pd[29]]; + (*drv)[cd[34]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] + - (*drv)[pd[13]] - (*drv)[pd[16]]) + + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] + + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] + - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[34]]); + break; + case 3: + (*drv)[cd[34]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] + - (*drv)[pd[13]] - (*drv)[pd[16]]) + + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] + + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] + - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[34]]); + break; + } + } else { + switch (lr_set) { + case 1: + (*drv)[cd[28]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) + - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[7]] + + 0.3125*((*drv)[pd[13]] - (*drv)[pd[31]]) + + 0.9375*(*drv)[pd[32]]); + (*drv)[cd[29]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.125*((*drv)[pd[7]] - (*drv)[pd[8]] - (*drv)[pd[13]]) + + 0.375*((*drv)[pd[14]] - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[33]]); + (*drv)[cd[30]] = (*drv)[pd[32]]; + (*drv)[cd[34]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] + - (*drv)[pd[13]] - (*drv)[pd[16]]) + + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] + + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] + - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[34]]); + break; + case 2: + (*drv)[cd[25]] = + (0.0390625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) + + 0.15625*((*drv)[pd[4]] + (*drv)[pd[6]]) + - 0.234375*(*drv)[pd[5]] + 0.0625*(*drv)[pd[10]] + + 0.3125*((*drv)[pd[16]] - (*drv)[pd[28]]) + + 0.9375*(*drv)[pd[29]]); + (*drv)[cd[26]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.125*((*drv)[pd[10]] - (*drv)[pd[11]] - (*drv)[pd[16]]) + + 0.375*((*drv)[pd[17]] - (*drv)[pd[28]] + (*drv)[pd[29]]) + + 0.75*(*drv)[pd[30]]); + (*drv)[cd[27]] = (*drv)[pd[29]]; + (*drv)[cd[34]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] + - (*drv)[pd[13]] - (*drv)[pd[16]]) + + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] + + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] + - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[34]]); + break; + case 3: + (*drv)[cd[34]] = + (-0.0390625*(*drv)[pd[0]] + 0.0234375*(*drv)[pd[1]] + + 0.09375*(*drv)[pd[4]] - 0.046875*(*drv)[pd[5]] + - 0.03125*(*drv)[pd[6]] + + 0.0625*((*drv)[pd[7]] + (*drv)[pd[10]] + - (*drv)[pd[13]] - (*drv)[pd[16]]) + + 0.375*(*drv)[pd[22]] - 0.125*(*drv)[pd[25]] + + 0.1875*(-(*drv)[pd[28]] + (*drv)[pd[29]] + - (*drv)[pd[31]] + (*drv)[pd[32]]) + + 0.75*(*drv)[pd[34]]); + break; } } + } } @@ -2586,16 +2597,11 @@ namespace AMDiS { { FUNCNAME("Lagrange::coarseRestr2_2d()"); - if (n < 1) - return; + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); - if (!drv->getFESpace()) { - ERROR("no fe_space in dof_real_vec\n"); - return; - } else if (!drv->getFESpace()->getBasisFcts()) { - ERROR("no basis functions in fe_space %s\n", drv->getFESpace()->getName().c_str()); + if (n < 1) return; - } Element *el = list->getElement(0); const DOFAdmin *admin = drv->getFESpace()->getAdmin(); @@ -2656,7 +2662,14 @@ namespace AMDiS { void Lagrange::coarseRestr2_3d(DOFIndexed<double> *drv, RCNeighbourList *list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseRestr2_3d"); + FUNCNAME("Lagrange::coarseRestr2_3d()"); + + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); + + if (n < 1) + return; + Element *el; const DegreeOfFreedom *cdof; DegreeOfFreedom pdof[10], cdofi; @@ -2664,18 +2677,8 @@ namespace AMDiS { int node0, n0; const DOFAdmin *admin; - if (n < 1) return; el = list->getElement(0); - if (!drv->getFESpace()) { - ERROR("no fe_space in dof_real_vec\n"); - return; - } else if (!drv->getFESpace()->getBasisFcts()) { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - admin = drv->getFESpace()->getAdmin(); basFct->getLocalIndices(el, admin, pdof); @@ -2720,53 +2723,55 @@ namespace AMDiS { /* adjust neighbour values */ /****************************************************************************/ - for (i = 1; i < n; i++) - { - el = list->getElement(i); - basFct->getLocalIndices(el, admin, pdof); - - 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] */ - /****************************************************************************/ - - cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL); - - switch(lr_set) - { - case 1: - cdofi = el->getChild(0)->getDOF(node0+4, n0); - (*drv)[pdof[0]] += (-0.125*(*drv)[cdofi]); - (*drv)[pdof[1]] += (-0.125*(*drv)[cdofi]); - (*drv)[pdof[4]] += (0.25*(*drv)[cdofi]); - (*drv)[pdof[5]] += (0.5*(*drv)[cdofi]); - (*drv)[pdof[7]] += (0.5*(*drv)[cdofi]); - break; - case 2: - cdofi = el->getChild(0)->getDOF(node0+5, n0); - (*drv)[pdof[0]] += (-0.125*(*drv)[cdofi]); - (*drv)[pdof[1]] += (-0.125*(*drv)[cdofi]); - (*drv)[pdof[4]] += (0.25*(*drv)[cdofi]); - (*drv)[pdof[6]] += (0.5*(*drv)[cdofi]); - (*drv)[pdof[8]] += (0.5*(*drv)[cdofi]); - break; - } - } - } + for (i = 1; i < n; i++) { + el = list->getElement(i); + basFct->getLocalIndices(el, admin, pdof); + lr_set = 0; + if (list->getNeighbourElement(i,0) && list->getNeighbourNr(i,0) < i) + lr_set = 1; - void Lagrange::coarseRestr3_1d(DOFIndexed<double> *drv, RCNeighbourList *list, - int n, BasisFunction* basFct) + 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] */ + /****************************************************************************/ + + cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL); + + switch (lr_set) { + case 1: + cdofi = el->getChild(0)->getDOF(node0+4, n0); + (*drv)[pdof[0]] += (-0.125*(*drv)[cdofi]); + (*drv)[pdof[1]] += (-0.125*(*drv)[cdofi]); + (*drv)[pdof[4]] += (0.25*(*drv)[cdofi]); + (*drv)[pdof[5]] += (0.5*(*drv)[cdofi]); + (*drv)[pdof[7]] += (0.5*(*drv)[cdofi]); + break; + case 2: + cdofi = el->getChild(0)->getDOF(node0+5, n0); + (*drv)[pdof[0]] += (-0.125*(*drv)[cdofi]); + (*drv)[pdof[1]] += (-0.125*(*drv)[cdofi]); + (*drv)[pdof[4]] += (0.25*(*drv)[cdofi]); + (*drv)[pdof[6]] += (0.5*(*drv)[cdofi]); + (*drv)[pdof[8]] += (0.5*(*drv)[cdofi]); + break; + } + } + } + + + void Lagrange::coarseRestr3_1d(DOFIndexed<double> *drv, RCNeighbourList *list, + int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseRestr3_1d"); + FUNCNAME("Lagrange::coarseRestr3_1d()"); + + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); + Element *el; int node, n0; const DegreeOfFreedom *cdof; @@ -2776,15 +2781,6 @@ namespace AMDiS { if (n < 1) return; el = list->getElement(0); - if (!drv->getFESpace()) { - ERROR("no fe_space in dof_real_vec\n"); - return; - } else if (!drv->getFESpace()->getBasisFcts()) { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - admin = drv->getFESpace()->getAdmin(); basFct->getLocalIndices(el, admin, pdof); @@ -2877,7 +2873,11 @@ namespace AMDiS { void Lagrange::coarseRestr3_2d(DOFIndexed<double> *drv, RCNeighbourList *list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseRestr3_2d"); + FUNCNAME("Lagrange::coarseRestr3_2d()"); + + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); + const Element *el; int node, n0; const DegreeOfFreedom *cdof; @@ -2887,17 +2887,6 @@ namespace AMDiS { if (n < 1) return; el = list->getElement(0); - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } admin = drv->getFESpace()->getAdmin(); basFct->getLocalIndices(el, admin, pdof); @@ -2986,29 +2975,23 @@ namespace AMDiS { void Lagrange::coarseRestr3_3d(DOFIndexed<double> *drv, RCNeighbourList *list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseRestr3_3d"); + FUNCNAME("Lagrange::coarseRestr3_3d()"); + + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); + + if (n < 1) + return; + const Element *el; const DegreeOfFreedom *cd; DegreeOfFreedom pd[19], cdi; int node0, n0, i, typ, lr_set; const DOFAdmin *admin; - if (n < 1) return; el = list->getElement(0); typ = list->getType(0); - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - admin = drv->getFESpace()->getAdmin(); basFct->getLocalIndices(el, admin, pd); @@ -3061,27 +3044,26 @@ 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)[pd[0]] += 0.0625*((*drv)[cd[8]] + (*drv)[cd[17]] + (*drv)[cd[18]]); - (*drv)[pd[1]] += - 0.3125*(*drv)[cd[8]] + 0.0625*(-(*drv)[cd[17]] - (*drv)[cd[18]]); - (*drv)[pd[4]] += - -0.3125*(*drv)[cd[8]] + 0.1875*(-(*drv)[cd[17]] - (*drv)[cd[18]]); - (*drv)[pd[5]] += - ((*drv)[cd[9]] + 0.9375*(*drv)[cd[8]] - + 0.1875*((*drv)[cd[17]] + (*drv)[cd[18]])); - (*drv)[pd[6]] += -0.125*(*drv)[cd[17]]; - (*drv)[pd[8]] += -0.125*(*drv)[cd[18]]; - (*drv)[pd[10]] += 0.375*(*drv)[cd[17]]; - (*drv)[pd[12]] += 0.375*(*drv)[cd[18]]; - (*drv)[pd[18]] += 0.75*(*drv)[cd[18]]; - (*drv)[pd[19]] += 0.750000*(*drv)[cd[17]]; - } else { + (*drv)[pd[0]] += 0.0625*((*drv)[cd[8]] + (*drv)[cd[17]] + (*drv)[cd[18]]); + (*drv)[pd[1]] += + 0.3125*(*drv)[cd[8]] + 0.0625*(-(*drv)[cd[17]] - (*drv)[cd[18]]); + (*drv)[pd[4]] += + -0.3125*(*drv)[cd[8]] + 0.1875*(-(*drv)[cd[17]] - (*drv)[cd[18]]); + (*drv)[pd[5]] += + ((*drv)[cd[9]] + 0.9375*(*drv)[cd[8]] + + 0.1875*((*drv)[cd[17]] + (*drv)[cd[18]])); + (*drv)[pd[6]] += -0.125*(*drv)[cd[17]]; + (*drv)[pd[8]] += -0.125*(*drv)[cd[18]]; + (*drv)[pd[10]] += 0.375*(*drv)[cd[17]]; + (*drv)[pd[12]] += 0.375*(*drv)[cd[18]]; + (*drv)[pd[18]] += 0.75*(*drv)[cd[18]]; + (*drv)[pd[19]] += 0.750000*(*drv)[cd[17]]; + } else { /****************************************************************************/ /* parent of el_type 1|2 */ /****************************************************************************/ @@ -3109,159 +3091,156 @@ 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] */ - /****************************************************************************/ - - cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); - cdi = cd[16]; - - switch(lr_set) - { - case 1: - (*drv)[pd[0]] += 0.0625*((*drv)[cd[12]] + (*drv)[cdi] - (*drv)[cd[18]]); - (*drv)[pd[1]] += 0.0625*((*drv)[cd[12]] + (*drv)[cdi] + (*drv)[cd[18]]); - (*drv)[pd[4]] += - (0.0625*(-(*drv)[cd[12]] - (*drv)[cdi]) + 0.1875*(*drv)[cd[18]]); - (*drv)[pd[5]] += - (0.0625*(-(*drv)[cd[12]] - (*drv)[cdi]) - 0.1875*(*drv)[cd[18]]); - (*drv)[pd[6]] += - (-0.25*(*drv)[cd[12]] - 0.125*(*drv)[cdi] + 0.375*(*drv)[cd[18]]); - (*drv)[pd[7]] += 0.5*(*drv)[cd[12]]; - (*drv)[pd[8]] += -0.125*(*drv)[cdi]; - (*drv)[pd[10]] += - -0.25*(*drv)[cd[12]] + 0.125*(-(*drv)[cdi] - (*drv)[cd[18]]); - (*drv)[pd[11]] += 0.5*(*drv)[cd[12]]; - (*drv)[pd[12]] += -0.125*(*drv)[cdi]; - (*drv)[pd[16]] += 0.5*(*drv)[cdi]; - (*drv)[pd[17]] += 0.5*(*drv)[cdi]; - (*drv)[pd[18]] += 0.25*(*drv)[cdi]; - (*drv)[pd[19]] = - ((*drv)[cd[13]] + 0.5*(*drv)[cd[12]] - + 0.25*(*drv)[cdi] + 0.75*(*drv)[cd[18]]); - break; - case 2: - (*drv)[pd[0]] += 0.0625*((*drv)[cd[14]] + (*drv)[cdi] - (*drv)[cd[17]]); - (*drv)[pd[1]] += 0.0625*((*drv)[cd[14]] + (*drv)[cdi] + (*drv)[cd[17]]); - (*drv)[pd[4]] += - 0.0625*(-(*drv)[cd[14]] - (*drv)[cdi]) + 0.1875*(*drv)[cd[17]]; - (*drv)[pd[5]] += - 0.0625*(-(*drv)[cd[14]] - (*drv)[cdi]) - 0.1875*(*drv)[cd[17]]; - (*drv)[pd[6]] += -0.125*(*drv)[cdi]; - (*drv)[pd[8]] += - -0.25*(*drv)[cd[14]] - 0.125*(*drv)[cdi] + 0.375*(*drv)[cd[17]]; - (*drv)[pd[9]] += 0.5*(*drv)[cd[14]]; - (*drv)[pd[10]] += -0.125*(*drv)[cdi]; - (*drv)[pd[12]] += - -0.25*(*drv)[cd[14]] + 0.125*(-(*drv)[cdi] - (*drv)[cd[17]]); - (*drv)[pd[13]] += 0.5*(*drv)[cd[14]]; - (*drv)[pd[16]] += 0.5*(*drv)[cdi]; - (*drv)[pd[17]] += 0.5*(*drv)[cdi]; - (*drv)[pd[18]] = - ((*drv)[cd[15]] + 0.5*(*drv)[cd[14]] - + 0.25*(*drv)[cdi] + 0.75*(*drv)[cd[17]]); - (*drv)[pd[19]] += 0.25*(*drv)[cdi]; - break; - case 3: - (*drv)[pd[0]] += 0.0625*(*drv)[cdi]; - (*drv)[pd[1]] += 0.0625*(*drv)[cdi]; - (*drv)[pd[4]] += -0.0625*(*drv)[cdi]; - (*drv)[pd[5]] += -0.0625*(*drv)[cdi]; - (*drv)[pd[6]] += -0.125*(*drv)[cdi]; - (*drv)[pd[8]] += -0.125*(*drv)[cdi]; - (*drv)[pd[10]] += -0.125*(*drv)[cdi]; - (*drv)[pd[12]] += -0.125*(*drv)[cdi]; - (*drv)[pd[16]] += 0.5*(*drv)[cdi]; - (*drv)[pd[17]] += 0.5*(*drv)[cdi]; - (*drv)[pd[18]] += 0.25*(*drv)[cdi]; - (*drv)[pd[19]] += 0.25*(*drv)[cdi]; - break; - } + for (i = 1; i < n; i++) { + el = list->getElement(i); + typ = list->getType(i); + basFct->getLocalIndices(el, admin, pd); - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ - - cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); - - if (typ == 0) - { - switch(lr_set) - { - case 1: - cdi = el->getChild(1)->getDOF(node0+1, n0); - TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n"); - (*drv)[pd[0]] += 0.0625*(*drv)[cdi]; - (*drv)[pd[1]] += -0.0625*(*drv)[cdi]; - (*drv)[pd[4]] += -0.1875*(*drv)[cdi]; - (*drv)[pd[5]] += 0.1875*(*drv)[cdi]; - (*drv)[pd[6]] += -0.125*(*drv)[cdi]; - (*drv)[pd[10]] += 0.375*(*drv)[cdi]; - (*drv)[pd[19]] += 0.75*(*drv)[cdi]; - break; - case 2: - cdi = el->getChild(1)->getDOF(node0+2, n0); - TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n"); - (*drv)[pd[0]] += 0.0625*(*drv)[cdi]; - (*drv)[pd[1]] += -0.0625*(*drv)[cdi]; - (*drv)[pd[4]] += -0.1875*(*drv)[cdi]; - (*drv)[pd[5]] += 0.1875*(*drv)[cdi]; - (*drv)[pd[8]] += -0.125*(*drv)[cdi]; - (*drv)[pd[12]] += 0.375*(*drv)[cdi]; - (*drv)[pd[18]] += 0.75*(*drv)[cdi]; - break; - } - } - else - { - switch(lr_set) - { - case 1: - cdi = el->getChild(1)->getDOF(node0+2, n0); - TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n"); - (*drv)[pd[0]] += 0.0625*(*drv)[cdi]; - (*drv)[pd[1]] += -0.0625*(*drv)[cdi]; - (*drv)[pd[4]] += -0.1875*(*drv)[cdi]; - (*drv)[pd[5]] += 0.1875*(*drv)[cdi]; - (*drv)[pd[6]] += -0.125*(*drv)[cdi]; - (*drv)[pd[10]] += 0.375*(*drv)[cdi]; - (*drv)[pd[19]] += 0.75*(*drv)[cdi]; - break; - case 2: - cdi = el->getChild(1)->getDOF(node0+1, n0); - TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n"); - (*drv)[pd[0]] += 0.0625*(*drv)[cdi]; - (*drv)[pd[1]] += -0.0625*(*drv)[cdi]; - (*drv)[pd[4]] += -0.1875*(*drv)[cdi]; - (*drv)[pd[5]] += 0.1875*(*drv)[cdi]; - (*drv)[pd[8]] += -0.125*(*drv)[cdi]; - (*drv)[pd[12]] += 0.375*(*drv)[cdi]; - (*drv)[pd[18]] += 0.75*(*drv)[cdi]; - break; - } - } + 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] */ + /****************************************************************************/ + + cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); + cdi = cd[16]; + + switch (lr_set) { + case 1: + (*drv)[pd[0]] += 0.0625*((*drv)[cd[12]] + (*drv)[cdi] - (*drv)[cd[18]]); + (*drv)[pd[1]] += 0.0625*((*drv)[cd[12]] + (*drv)[cdi] + (*drv)[cd[18]]); + (*drv)[pd[4]] += + (0.0625*(-(*drv)[cd[12]] - (*drv)[cdi]) + 0.1875*(*drv)[cd[18]]); + (*drv)[pd[5]] += + (0.0625*(-(*drv)[cd[12]] - (*drv)[cdi]) - 0.1875*(*drv)[cd[18]]); + (*drv)[pd[6]] += + (-0.25*(*drv)[cd[12]] - 0.125*(*drv)[cdi] + 0.375*(*drv)[cd[18]]); + (*drv)[pd[7]] += 0.5*(*drv)[cd[12]]; + (*drv)[pd[8]] += -0.125*(*drv)[cdi]; + (*drv)[pd[10]] += + -0.25*(*drv)[cd[12]] + 0.125*(-(*drv)[cdi] - (*drv)[cd[18]]); + (*drv)[pd[11]] += 0.5*(*drv)[cd[12]]; + (*drv)[pd[12]] += -0.125*(*drv)[cdi]; + (*drv)[pd[16]] += 0.5*(*drv)[cdi]; + (*drv)[pd[17]] += 0.5*(*drv)[cdi]; + (*drv)[pd[18]] += 0.25*(*drv)[cdi]; + (*drv)[pd[19]] = + ((*drv)[cd[13]] + 0.5*(*drv)[cd[12]] + + 0.25*(*drv)[cdi] + 0.75*(*drv)[cd[18]]); + break; + case 2: + (*drv)[pd[0]] += 0.0625*((*drv)[cd[14]] + (*drv)[cdi] - (*drv)[cd[17]]); + (*drv)[pd[1]] += 0.0625*((*drv)[cd[14]] + (*drv)[cdi] + (*drv)[cd[17]]); + (*drv)[pd[4]] += + 0.0625*(-(*drv)[cd[14]] - (*drv)[cdi]) + 0.1875*(*drv)[cd[17]]; + (*drv)[pd[5]] += + 0.0625*(-(*drv)[cd[14]] - (*drv)[cdi]) - 0.1875*(*drv)[cd[17]]; + (*drv)[pd[6]] += -0.125*(*drv)[cdi]; + (*drv)[pd[8]] += + -0.25*(*drv)[cd[14]] - 0.125*(*drv)[cdi] + 0.375*(*drv)[cd[17]]; + (*drv)[pd[9]] += 0.5*(*drv)[cd[14]]; + (*drv)[pd[10]] += -0.125*(*drv)[cdi]; + (*drv)[pd[12]] += + -0.25*(*drv)[cd[14]] + 0.125*(-(*drv)[cdi] - (*drv)[cd[17]]); + (*drv)[pd[13]] += 0.5*(*drv)[cd[14]]; + (*drv)[pd[16]] += 0.5*(*drv)[cdi]; + (*drv)[pd[17]] += 0.5*(*drv)[cdi]; + (*drv)[pd[18]] = + ((*drv)[cd[15]] + 0.5*(*drv)[cd[14]] + + 0.25*(*drv)[cdi] + 0.75*(*drv)[cd[17]]); + (*drv)[pd[19]] += 0.25*(*drv)[cdi]; + break; + case 3: + (*drv)[pd[0]] += 0.0625*(*drv)[cdi]; + (*drv)[pd[1]] += 0.0625*(*drv)[cdi]; + (*drv)[pd[4]] += -0.0625*(*drv)[cdi]; + (*drv)[pd[5]] += -0.0625*(*drv)[cdi]; + (*drv)[pd[6]] += -0.125*(*drv)[cdi]; + (*drv)[pd[8]] += -0.125*(*drv)[cdi]; + (*drv)[pd[10]] += -0.125*(*drv)[cdi]; + (*drv)[pd[12]] += -0.125*(*drv)[cdi]; + (*drv)[pd[16]] += 0.5*(*drv)[cdi]; + (*drv)[pd[17]] += 0.5*(*drv)[cdi]; + (*drv)[pd[18]] += 0.25*(*drv)[cdi]; + (*drv)[pd[19]] += 0.25*(*drv)[cdi]; + break; + } + + /****************************************************************************/ + /* values on child[1] */ + /****************************************************************************/ + + cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); + + if (typ == 0) { + switch (lr_set) { + case 1: + cdi = el->getChild(1)->getDOF(node0+1, n0); + TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n"); + (*drv)[pd[0]] += 0.0625*(*drv)[cdi]; + (*drv)[pd[1]] += -0.0625*(*drv)[cdi]; + (*drv)[pd[4]] += -0.1875*(*drv)[cdi]; + (*drv)[pd[5]] += 0.1875*(*drv)[cdi]; + (*drv)[pd[6]] += -0.125*(*drv)[cdi]; + (*drv)[pd[10]] += 0.375*(*drv)[cdi]; + (*drv)[pd[19]] += 0.75*(*drv)[cdi]; + break; + case 2: + cdi = el->getChild(1)->getDOF(node0+2, n0); + TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n"); + (*drv)[pd[0]] += 0.0625*(*drv)[cdi]; + (*drv)[pd[1]] += -0.0625*(*drv)[cdi]; + (*drv)[pd[4]] += -0.1875*(*drv)[cdi]; + (*drv)[pd[5]] += 0.1875*(*drv)[cdi]; + (*drv)[pd[8]] += -0.125*(*drv)[cdi]; + (*drv)[pd[12]] += 0.375*(*drv)[cdi]; + (*drv)[pd[18]] += 0.75*(*drv)[cdi]; + break; + } + } else { + switch (lr_set) { + case 1: + cdi = el->getChild(1)->getDOF(node0+2, n0); + TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n"); + (*drv)[pd[0]] += 0.0625*(*drv)[cdi]; + (*drv)[pd[1]] += -0.0625*(*drv)[cdi]; + (*drv)[pd[4]] += -0.1875*(*drv)[cdi]; + (*drv)[pd[5]] += 0.1875*(*drv)[cdi]; + (*drv)[pd[6]] += -0.125*(*drv)[cdi]; + (*drv)[pd[10]] += 0.375*(*drv)[cdi]; + (*drv)[pd[19]] += 0.75*(*drv)[cdi]; + break; + case 2: + cdi = el->getChild(1)->getDOF(node0+1, n0); + TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n"); + (*drv)[pd[0]] += 0.0625*(*drv)[cdi]; + (*drv)[pd[1]] += -0.0625*(*drv)[cdi]; + (*drv)[pd[4]] += -0.1875*(*drv)[cdi]; + (*drv)[pd[5]] += 0.1875*(*drv)[cdi]; + (*drv)[pd[8]] += -0.125*(*drv)[cdi]; + (*drv)[pd[12]] += 0.375*(*drv)[cdi]; + (*drv)[pd[18]] += 0.75*(*drv)[cdi]; + break; + } } + } } void Lagrange::coarseRestr4_1d(DOFIndexed<double> *drv, RCNeighbourList *list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseRestr4_1d"); + FUNCNAME("Lagrange::coarseRestr4_1d()"); + + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); + const Element *el; DegreeOfFreedom pdof[5]; const DegreeOfFreedom *cdof; @@ -3270,18 +3249,6 @@ namespace AMDiS { if (n < 1) return; el = list->getElement(0); - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - admin = drv->getFESpace()->getAdmin(); basFct->getLocalIndices(el, admin, pdof); @@ -3434,7 +3401,11 @@ namespace AMDiS { void Lagrange::coarseRestr4_2d(DOFIndexed<double> *drv, RCNeighbourList *list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseRestr4_2d"); + FUNCNAME("Lagrange::coarseRestr4_2d()"); + + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); + const Element *el; DegreeOfFreedom pdof[15]; const DegreeOfFreedom *cdof; @@ -3443,18 +3414,6 @@ namespace AMDiS { if (n < 1) return; el = list->getElement(0); - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - admin = drv->getFESpace()->getAdmin(); basFct->getLocalIndices(el, admin, pdof); @@ -3607,7 +3566,11 @@ namespace AMDiS { void Lagrange::coarseRestr4_3d(DOFIndexed<double> *drv, RCNeighbourList *list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseRestr4_3d"); + FUNCNAME("Lagrange::coarseRestr4_3d()"); + + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); + DegreeOfFreedom pd[32]; const DegreeOfFreedom *cd; const Element *el; @@ -3618,18 +3581,6 @@ namespace AMDiS { el = list->getElement(0); typ = list->getType(0); - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - admin = drv->getFESpace()->getAdmin(); basFct->getLocalIndices(el, admin, pd); @@ -3751,64 +3702,63 @@ 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)[pd[0]] += - (0.0390625*(-(*drv)[cd[10]] - (*drv)[cd[25]] - (*drv)[cd[26]] - - (*drv)[cd[28]] - (*drv)[cd[29]] - (*drv)[cd[34]]) - + 0.0234375*(*drv)[cd[12]]); - (*drv)[pd[1]] += - (0.2734375*(*drv)[cd[10]] - + 0.0390625*(-(*drv)[cd[12]] - (*drv)[cd[25]] - (*drv)[cd[28]]) - + 0.0234375*((*drv)[cd[26]] + (*drv)[cd[29]] + (*drv)[cd[34]])); - (*drv)[pd[4]] += - (0.21875*(*drv)[cd[10]] - + 0.15625*(-(*drv)[cd[12]] + (*drv)[cd[25]] + (*drv)[cd[28]]) - + 0.09375*((*drv)[cd[26]] + (*drv)[cd[29]] + (*drv)[cd[34]])); - (*drv)[pd[5]] += - (-0.546875*(*drv)[cd[10]] + 0.703125*(*drv)[cd[12]] - + 0.234375*(-(*drv)[cd[25]] - (*drv)[cd[28]]) - + 0.046875*(-(*drv)[cd[26]] - (*drv)[cd[29]] - (*drv)[cd[34]])); - (*drv)[pd[6]] += - ((*drv)[cd[11]] + 1.09375*(*drv)[cd[10]] + 0.46875*(*drv)[cd[12]] - + 0.15625*((*drv)[cd[25]] + (*drv)[cd[28]]) - + 0.03125*(-(*drv)[cd[26]] - (*drv)[cd[29]] - (*drv)[cd[34]])); - (*drv)[pd[7]] += - (0.0625*((*drv)[cd[25]] + (*drv)[cd[34]]) + 0.125*(*drv)[cd[26]]); - (*drv)[pd[8]] += -0.125*(*drv)[cd[26]]; - (*drv)[pd[10]] += - (0.0625*((*drv)[cd[28]] + (*drv)[cd[34]]) + 0.125*(*drv)[cd[29]]); - (*drv)[pd[11]] += -0.125*(*drv)[cd[29]]; - (*drv)[pd[13]] += - (0.3125*(*drv)[cd[25]] - 0.125*(*drv)[cd[26]] - - 0.0625*(*drv)[cd[34]]); - (*drv)[pd[14]] += 0.375*(*drv)[cd[26]]; - (*drv)[pd[16]] += - (0.3125*(*drv)[cd[28]] - 0.125*(*drv)[cd[29]] - - 0.0625*(*drv)[cd[34]]); - (*drv)[pd[17]] += 0.375*(*drv)[cd[29]]; - (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; - (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; - (*drv)[pd[28]] += - (-0.3125*(*drv)[cd[28]] - 0.375*(*drv)[cd[29]] - - 0.1875*(*drv)[cd[34]]); - (*drv)[pd[29]] += - ((*drv)[cd[30]] + 0.9375*(*drv)[cd[28]] + 0.375*(*drv)[cd[29]] - + 0.1875*(*drv)[cd[34]]); - (*drv)[pd[30]] += 0.75*(*drv)[cd[29]]; - (*drv)[pd[31]] += - (-0.3125*(*drv)[cd[25]] - 0.375*(*drv)[cd[26]] - - 0.1875*(*drv)[cd[34]]); - (*drv)[pd[32]] += - ((*drv)[cd[27]] + 0.9375*(*drv)[cd[25]] - + 0.375*(*drv)[cd[26]] + 0.1875*(*drv)[cd[34]]); - (*drv)[pd[33]] += 0.75*(*drv)[cd[26]]; - (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; - } else { + (*drv)[pd[0]] += + (0.0390625*(-(*drv)[cd[10]] - (*drv)[cd[25]] - (*drv)[cd[26]] + - (*drv)[cd[28]] - (*drv)[cd[29]] - (*drv)[cd[34]]) + + 0.0234375*(*drv)[cd[12]]); + (*drv)[pd[1]] += + (0.2734375*(*drv)[cd[10]] + + 0.0390625*(-(*drv)[cd[12]] - (*drv)[cd[25]] - (*drv)[cd[28]]) + + 0.0234375*((*drv)[cd[26]] + (*drv)[cd[29]] + (*drv)[cd[34]])); + (*drv)[pd[4]] += + (0.21875*(*drv)[cd[10]] + + 0.15625*(-(*drv)[cd[12]] + (*drv)[cd[25]] + (*drv)[cd[28]]) + + 0.09375*((*drv)[cd[26]] + (*drv)[cd[29]] + (*drv)[cd[34]])); + (*drv)[pd[5]] += + (-0.546875*(*drv)[cd[10]] + 0.703125*(*drv)[cd[12]] + + 0.234375*(-(*drv)[cd[25]] - (*drv)[cd[28]]) + + 0.046875*(-(*drv)[cd[26]] - (*drv)[cd[29]] - (*drv)[cd[34]])); + (*drv)[pd[6]] += + ((*drv)[cd[11]] + 1.09375*(*drv)[cd[10]] + 0.46875*(*drv)[cd[12]] + + 0.15625*((*drv)[cd[25]] + (*drv)[cd[28]]) + + 0.03125*(-(*drv)[cd[26]] - (*drv)[cd[29]] - (*drv)[cd[34]])); + (*drv)[pd[7]] += + (0.0625*((*drv)[cd[25]] + (*drv)[cd[34]]) + 0.125*(*drv)[cd[26]]); + (*drv)[pd[8]] += -0.125*(*drv)[cd[26]]; + (*drv)[pd[10]] += + (0.0625*((*drv)[cd[28]] + (*drv)[cd[34]]) + 0.125*(*drv)[cd[29]]); + (*drv)[pd[11]] += -0.125*(*drv)[cd[29]]; + (*drv)[pd[13]] += + (0.3125*(*drv)[cd[25]] - 0.125*(*drv)[cd[26]] + - 0.0625*(*drv)[cd[34]]); + (*drv)[pd[14]] += 0.375*(*drv)[cd[26]]; + (*drv)[pd[16]] += + (0.3125*(*drv)[cd[28]] - 0.125*(*drv)[cd[29]] + - 0.0625*(*drv)[cd[34]]); + (*drv)[pd[17]] += 0.375*(*drv)[cd[29]]; + (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; + (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; + (*drv)[pd[28]] += + (-0.3125*(*drv)[cd[28]] - 0.375*(*drv)[cd[29]] + - 0.1875*(*drv)[cd[34]]); + (*drv)[pd[29]] += + ((*drv)[cd[30]] + 0.9375*(*drv)[cd[28]] + 0.375*(*drv)[cd[29]] + + 0.1875*(*drv)[cd[34]]); + (*drv)[pd[30]] += 0.75*(*drv)[cd[29]]; + (*drv)[pd[31]] += + (-0.3125*(*drv)[cd[25]] - 0.375*(*drv)[cd[26]] + - 0.1875*(*drv)[cd[34]]); + (*drv)[pd[32]] += + ((*drv)[cd[27]] + 0.9375*(*drv)[cd[25]] + + 0.375*(*drv)[cd[26]] + 0.1875*(*drv)[cd[34]]); + (*drv)[pd[33]] += 0.75*(*drv)[cd[26]]; + (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; + } else { /****************************************************************************/ /* parent of el_type 1|2 */ /****************************************************************************/ @@ -3870,618 +3820,548 @@ namespace AMDiS { /* adjust neighbour values */ /****************************************************************************/ - 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] */ - /****************************************************************************/ - - cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); - - switch(lr_set) - { - case 1: - (*drv)[pd[0]] += - (0.0390625*(-(*drv)[cd[16]] - (*drv)[cd[22]] - - (*drv)[cd[23]] - (*drv)[cd[28]]) - + 0.0234375*((*drv)[cd[18]] + (*drv)[cd[29]] + (*drv)[cd[34]])); - (*drv)[pd[1]] += - (0.0390625*(-(*drv)[cd[16]] - (*drv)[cd[22]] - (*drv)[cd[23]] - - (*drv)[cd[28]] - (*drv)[cd[29]] - (*drv)[cd[34]]) - + 0.0234375*(*drv)[cd[18]]); - (*drv)[pd[4]] += - (0.03125*((*drv)[cd[16]] + (*drv)[cd[22]] + (*drv)[cd[23]] - - (*drv)[cd[29]] - (*drv)[cd[34]]) - - 0.09375*(*drv)[cd[18]] + 0.15625*(*drv)[cd[28]]); - (*drv)[pd[5]] += - (0.015625*((*drv)[cd[16]] + (*drv)[cd[22]] + (*drv)[cd[23]]) - + 0.140625*(*drv)[cd[18]] - 0.234375*(*drv)[cd[28]] - + 0.046875*(-(*drv)[cd[29]] - (*drv)[cd[34]])); - (*drv)[pd[6]] += - (0.03125*((*drv)[cd[16]] + (*drv)[cd[22]] + (*drv)[cd[23]]) - + 0.09375*(-(*drv)[cd[18]] + (*drv)[cd[29]] + (*drv)[cd[34]]) - + 0.15625*(*drv)[cd[28]]); - (*drv)[pd[7]] += - (0.1875*(*drv)[cd[16]] - + 0.0625*(-(*drv)[cd[18]] + (*drv)[cd[23]] - (*drv)[cd[34]]) - + 0.125*((*drv)[cd[22]] - (*drv)[cd[29]]) - + 0.3125*(*drv)[cd[28]]); - (*drv)[pd[8]] += - (0.375*(-(*drv)[cd[16]] + (*drv)[cd[29]]) - - 0.125*(*drv)[cd[22]]); - (*drv)[pd[9]] += 0.5*(*drv)[cd[16]]; - (*drv)[pd[10]] += - (0.0625*((*drv)[cd[22]] - (*drv)[cd[34]]) + 0.125*(*drv)[cd[23]]); - (*drv)[pd[11]] += -0.125*(*drv)[cd[23]]; - (*drv)[pd[13]] += - (0.1875*(*drv)[cd[16]] - + 0.0625*(-(*drv)[cd[18]] + (*drv)[cd[23]] - + (*drv)[cd[28]] + (*drv)[cd[34]]) - + 0.125*((*drv)[cd[22]] + (*drv)[cd[29]])); - (*drv)[pd[14]] += - (-0.375*(*drv)[cd[16]] - + 0.125*(-(*drv)[cd[22]] - (*drv)[cd[29]])); - (*drv)[pd[15]] += 0.5*(*drv)[cd[16]]; - (*drv)[pd[16]] += - (0.0625*((*drv)[cd[22]] + (*drv)[cd[34]]) + 0.125*(*drv)[cd[23]]); - (*drv)[pd[17]] += -0.125*(*drv)[cd[23]]; - (*drv)[pd[22]] += - (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) - 0.125*(*drv)[cd[34]]); - (*drv)[pd[23]] += 0.5*(*drv)[cd[22]]; - (*drv)[pd[24]] += 0.5*(*drv)[cd[23]]; - (*drv)[pd[25]] += - (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) + 0.375*(*drv)[cd[34]]); - (*drv)[pd[26]] += 0.5*(*drv)[cd[22]]; - (*drv)[pd[27]] += 0.5*(*drv)[cd[23]]; - (*drv)[pd[28]] += - (-0.0625*(*drv)[cd[22]] - 0.125*(*drv)[cd[23]] - + 0.1875*(*drv)[cd[34]]); - (*drv)[pd[29]] += - (-0.0625*(*drv)[cd[22]] - 0.125*(*drv)[cd[23]] - - 0.1875*(*drv)[cd[34]]); - (*drv)[pd[30]] += 0.25*(*drv)[cd[23]]; - (*drv)[pd[31]] = - ((*drv)[cd[30]] + 0.1875*(-(*drv)[cd[16]] + (*drv)[cd[34]]) - + 0.5625*(*drv)[cd[18]] - 0.125*(*drv)[cd[22]] - - 0.0625*(*drv)[cd[23]] + 0.9375*(*drv)[cd[28]] - + 0.375*(*drv)[cd[29]]); - (*drv)[pd[32]] = - (0.1875*(-(*drv)[cd[16]] - (*drv)[cd[34]]) - + 0.5625*(*drv)[cd[18]] - 0.125*(*drv)[cd[22]] - - 0.0625*(*drv)[cd[23]] - 0.3125*(*drv)[cd[28]] - - 0.375*(*drv)[cd[29]]); - (*drv)[pd[33]] = - ((*drv)[cd[17]] + 0.75*((*drv)[cd[16]] + (*drv)[cd[29]]) - + 0.25*(*drv)[cd[22]]); - (*drv)[pd[34]] = - ((*drv)[cd[24]] + 0.5*((*drv)[cd[22]] + (*drv)[cd[23]]) - + 0.75*(*drv)[cd[34]]); - break; - case 2: - (*drv)[pd[0]] += - (0.0390625*(-(*drv)[cd[19]] - (*drv)[cd[22]] - - (*drv)[cd[23]] - (*drv)[cd[25]]) - + 0.0234375*((*drv)[cd[21]] + (*drv)[cd[26]] + (*drv)[cd[34]])); - (*drv)[pd[1]] += - (0.0390625*(-(*drv)[cd[19]] - (*drv)[cd[22]] - (*drv)[cd[23]] - - (*drv)[cd[25]] - (*drv)[cd[26]] - (*drv)[cd[34]]) - + 0.0234375*(*drv)[cd[21]]); - (*drv)[pd[4]] += - (0.03125*((*drv)[cd[19]] + (*drv)[cd[22]] + (*drv)[cd[23]] - - (*drv)[cd[26]] - (*drv)[cd[34]]) - - 0.09375*(*drv)[cd[21]] + 0.15625*(*drv)[cd[25]]); - (*drv)[pd[5]] += - (0.015625*((*drv)[cd[19]] + (*drv)[cd[22]] + (*drv)[cd[23]]) - + 0.140625*(*drv)[cd[21]] - 0.234375*(*drv)[cd[25]] - + 0.046875*(-(*drv)[cd[26]] - (*drv)[cd[34]])); - (*drv)[pd[6]] += - (0.03125*((*drv)[cd[19]] + (*drv)[cd[22]] + (*drv)[cd[23]]) - + 0.09375*(-(*drv)[cd[21]] + (*drv)[cd[26]] + (*drv)[cd[34]]) - + 0.15625*(*drv)[cd[25]]); - (*drv)[pd[7]] += - (0.125*(*drv)[cd[22]] + 0.0625*((*drv)[cd[23]] - (*drv)[cd[34]])); - (*drv)[pd[8]] += -0.125*(*drv)[cd[22]]; - (*drv)[pd[10]] += - (0.1875*(*drv)[cd[19]] - + 0.0625*(-(*drv)[cd[21]] + (*drv)[cd[22]] - (*drv)[cd[34]]) - + 0.125*((*drv)[cd[23]] - (*drv)[cd[26]]) - + 0.3125*(*drv)[cd[25]]); - (*drv)[pd[11]] += - (0.375*(-(*drv)[cd[19]] + (*drv)[cd[26]]) - 0.125*(*drv)[cd[23]]); - (*drv)[pd[12]] += 0.5*(*drv)[cd[19]]; - (*drv)[pd[13]] += - (0.125*(*drv)[cd[22]] + 0.0625*((*drv)[cd[23]] + (*drv)[cd[34]])); - (*drv)[pd[14]] += -0.125*(*drv)[cd[22]]; - (*drv)[pd[16]] += - (0.1875*(*drv)[cd[19]] - + 0.0625*(-(*drv)[cd[21]] + (*drv)[cd[22]] - + (*drv)[cd[25]] + (*drv)[cd[34]]) - + 0.125*((*drv)[cd[23]] + (*drv)[cd[26]])); - (*drv)[pd[17]] += - (-0.375*(*drv)[cd[19]] - + 0.125*(-(*drv)[cd[23]] - (*drv)[cd[26]])); - (*drv)[pd[18]] += 0.5*(*drv)[cd[19]]; - (*drv)[pd[22]] += - (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) - 0.125*(*drv)[cd[34]]); - (*drv)[pd[23]] += (0.5*(*drv)[cd[22]]); - (*drv)[pd[24]] += (0.5*(*drv)[cd[23]]); - (*drv)[pd[25]] += - (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) + 0.375*(*drv)[cd[34]]); - (*drv)[pd[26]] += (0.5*(*drv)[cd[22]]); - (*drv)[pd[27]] += (0.5*(*drv)[cd[23]]); - (*drv)[pd[28]] = - ((*drv)[cd[27]] + 0.1875*(-(*drv)[cd[19]] + (*drv)[cd[34]]) - + 0.5625*(*drv)[cd[21]] - 0.0625*(*drv)[cd[22]] - - 0.125*(*drv)[cd[23]] + 0.9375*(*drv)[cd[25]] - + 0.375*(*drv)[cd[26]]); - (*drv)[pd[29]] = - (0.1875*(-(*drv)[cd[19]] - (*drv)[cd[34]]) - + 0.5625*(*drv)[cd[21]] - 0.0625*(*drv)[cd[22]] - - 0.125*(*drv)[cd[23]] - 0.3125*(*drv)[cd[25]] - - 0.375*(*drv)[cd[26]]); - (*drv)[pd[30]] = - ((*drv)[cd[20]] + 0.75*((*drv)[cd[19]] + (*drv)[cd[26]]) - + 0.25*(*drv)[cd[23]]); - (*drv)[pd[31]] += - (-0.125*(*drv)[cd[22]] - 0.0625*(*drv)[cd[23]] - + 0.1875*(*drv)[cd[34]]); - (*drv)[pd[32]] += - (-0.125*(*drv)[cd[22]] - 0.0625*(*drv)[cd[23]] - - 0.1875*(*drv)[cd[34]]); - (*drv)[pd[33]] += 0.25*(*drv)[cd[22]]; - (*drv)[pd[34]] = - ((*drv)[cd[24]] + 0.5*((*drv)[cd[22]] + (*drv)[cd[23]]) - + 0.75*(*drv)[cd[34]]); - break; - case 3: - (*drv)[pd[0]] += - (0.0390625*(-(*drv)[cd[22]] - (*drv)[cd[23]]) - + 0.0234375*(*drv)[cd[34]]); - (*drv)[pd[1]] += - (0.0390625*(-(*drv)[cd[22]] - (*drv)[cd[23]] - (*drv)[cd[34]])); - (*drv)[pd[4]] += - (0.03125*((*drv)[cd[22]] + (*drv)[cd[23]] - (*drv)[cd[34]])); - (*drv)[pd[5]] += - (0.015625*((*drv)[cd[22]] + (*drv)[cd[23]]) - - 0.046875*(*drv)[cd[34]]); - (*drv)[pd[6]] += - (0.03125*((*drv)[cd[22]] + (*drv)[cd[23]]) - + 0.09375*(*drv)[cd[34]]); - (*drv)[pd[7]] += - (0.125*(*drv)[cd[22]] + 0.0625*((*drv)[cd[23]] - (*drv)[cd[34]])); - (*drv)[pd[8]] += -0.125*(*drv)[cd[22]]; - (*drv)[pd[10]] += - (0.0625*((*drv)[cd[22]] - (*drv)[cd[34]]) + 0.125*(*drv)[cd[23]]); - (*drv)[pd[11]] += -0.125*(*drv)[cd[23]]; - (*drv)[pd[13]] += - (0.125*(*drv)[cd[22]] + 0.0625*((*drv)[cd[23]] + (*drv)[cd[34]])); - (*drv)[pd[14]] += -0.125*(*drv)[cd[22]]; - (*drv)[pd[16]] += - (0.0625*((*drv)[cd[22]] + (*drv)[cd[34]]) + 0.125*(*drv)[cd[23]]); - (*drv)[pd[17]] += -0.125*(*drv)[cd[23]]; - (*drv)[pd[22]] += - (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) - 0.125*(*drv)[cd[34]]); - (*drv)[pd[23]] += 0.5*(*drv)[cd[22]]; - (*drv)[pd[24]] += 0.5*(*drv)[cd[23]]; - (*drv)[pd[25]] += - (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) + 0.375*(*drv)[cd[34]]); - (*drv)[pd[26]] += 0.5*(*drv)[cd[22]]; - (*drv)[pd[27]] += 0.5*(*drv)[cd[23]]; - (*drv)[pd[28]] += - (-0.0625*(*drv)[cd[22]] - 0.125*(*drv)[cd[23]] - + 0.1875*(*drv)[cd[34]]); - (*drv)[pd[29]] += - (-0.0625*(*drv)[cd[22]] - 0.125*(*drv)[cd[23]] - - 0.1875*(*drv)[cd[34]]); - (*drv)[pd[30]] += 0.25*(*drv)[cd[23]]; - (*drv)[pd[31]] += - (-0.125*(*drv)[cd[22]] - 0.0625*(*drv)[cd[23]] - + 0.1875*(*drv)[cd[34]]); - (*drv)[pd[32]] += - (-0.125*(*drv)[cd[22]] - 0.0625*(*drv)[cd[23]] - - 0.1875*(*drv)[cd[34]]); - (*drv)[pd[33]] += 0.25*(*drv)[cd[22]]; - (*drv)[pd[34]] = - ((*drv)[cd[24]] + 0.5*((*drv)[cd[22]] + (*drv)[cd[23]]) - + 0.75*(*drv)[cd[34]]); - break; - } - - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ - - cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); - - if (typ == 0) { - switch(lr_set) { - case 1: - (*drv)[pd[0]] += - (0.0390625*(-(*drv)[cd[25]] - (*drv)[cd[26]] - (*drv)[cd[34]])); - (*drv)[pd[1]] += - (-0.0390625*(*drv)[cd[25]] - + 0.0234375*((*drv)[cd[26]] + (*drv)[cd[34]])); - (*drv)[pd[4]] += - (0.15625*(*drv)[cd[25]] - + 0.09375*((*drv)[cd[26]] + (*drv)[cd[34]])); - (*drv)[pd[5]] += - (-0.234375*(*drv)[cd[25]] - + 0.046875*(-(*drv)[cd[26]] - (*drv)[cd[34]])); - (*drv)[pd[6]] += - (0.15625*(*drv)[cd[25]] - + 0.03125*(-(*drv)[cd[26]] - (*drv)[cd[34]])); - (*drv)[pd[7]] += - (0.0625*((*drv)[cd[25]] + (*drv)[cd[34]]) - + 0.125*(*drv)[cd[26]]); - (*drv)[pd[8]] += -0.125*(*drv)[cd[26]]; - (*drv)[pd[10]] += 0.0625*(*drv)[cd[34]]; - (*drv)[pd[13]] += - (0.3125*(*drv)[cd[25]] - 0.125*(*drv)[cd[26]] - - 0.0625*(*drv)[cd[34]]); - (*drv)[pd[14]] += 0.375*(*drv)[cd[26]]; - (*drv)[pd[16]] += -0.0625*(*drv)[cd[34]]; - (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; - (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; - (*drv)[pd[28]] += -0.1875*(*drv)[cd[34]]; - (*drv)[pd[29]] += 0.1875*(*drv)[cd[34]]; - (*drv)[pd[31]] += - (-0.3125*(*drv)[cd[25]] - 0.375*(*drv)[cd[26]] - - 0.1875*(*drv)[cd[34]]); - (*drv)[pd[32]] += - ((*drv)[cd[27]] + 0.9375*(*drv)[cd[25]] - + 0.375*(*drv)[cd[26]] + 0.1875*(*drv)[cd[34]]); - (*drv)[pd[33]] += 0.75*(*drv)[cd[26]]; - (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; - break; - case 2: - (*drv)[pd[0]] += - (0.0390625*(-(*drv)[cd[28]] - (*drv)[cd[29]] - (*drv)[cd[34]])); - (*drv)[pd[1]] += - (-0.0390625*(*drv)[cd[28]] - + 0.0234375*((*drv)[cd[29]] + (*drv)[cd[34]])); - (*drv)[pd[4]] += - (0.15625*(*drv)[cd[28]] - + 0.09375*((*drv)[cd[29]] + (*drv)[cd[34]])); - (*drv)[pd[5]] += - (-0.234375*(*drv)[cd[28]] - + 0.046875*(-(*drv)[cd[29]] - (*drv)[cd[34]])); - (*drv)[pd[6]] += - (0.15625*(*drv)[cd[28]] - + 0.03125*(-(*drv)[cd[29]] - (*drv)[cd[34]])); - (*drv)[pd[7]] += 0.0625*(*drv)[cd[34]]; - (*drv)[pd[10]] += - (0.0625*((*drv)[cd[28]] + (*drv)[cd[34]]) - + 0.125*(*drv)[cd[29]]); - (*drv)[pd[11]] += -0.125*(*drv)[cd[29]]; - (*drv)[pd[13]] += -0.0625*(*drv)[cd[34]]; - (*drv)[pd[16]] += - (0.3125*(*drv)[cd[28]] - 0.125*(*drv)[cd[29]] - - 0.0625*(*drv)[cd[34]]); - (*drv)[pd[17]] += 0.375*(*drv)[cd[29]]; - (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; - (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; - (*drv)[pd[28]] += - (-0.3125*(*drv)[cd[28]] - 0.375*(*drv)[cd[29]] - - 0.1875*(*drv)[cd[34]]); - (*drv)[pd[29]] += - ((*drv)[cd[30]] + 0.9375*(*drv)[cd[28]] + 0.375*(*drv)[cd[29]] - + 0.1875*(*drv)[cd[34]]); - (*drv)[pd[30]] += 0.75*(*drv)[cd[29]]; - (*drv)[pd[31]] += -0.1875*(*drv)[cd[34]]; - (*drv)[pd[32]] += 0.1875*(*drv)[cd[34]]; - (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; - break; - case 3: - (*drv)[pd[0]] += -0.0390625*(*drv)[cd[34]]; - (*drv)[pd[1]] += 0.0234375*(*drv)[cd[34]]; - (*drv)[pd[4]] += 0.09375*(*drv)[cd[34]]; - (*drv)[pd[5]] += -0.046875*(*drv)[cd[34]]; - (*drv)[pd[6]] += -0.03125*(*drv)[cd[34]]; - (*drv)[pd[7]] += 0.0625*(*drv)[cd[34]]; - (*drv)[pd[10]] += 0.0625*(*drv)[cd[34]]; - (*drv)[pd[13]] += -0.0625*(*drv)[cd[34]]; - (*drv)[pd[16]] += -0.0625*(*drv)[cd[34]]; - (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; - (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; - (*drv)[pd[28]] += -0.1875*(*drv)[cd[34]]; - (*drv)[pd[29]] += 0.1875*(*drv)[cd[34]]; - (*drv)[pd[31]] += -0.1875*(*drv)[cd[34]]; - (*drv)[pd[32]] += 0.1875*(*drv)[cd[34]]; - (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; - break; - } - } - else { - switch(lr_set) { - case 1: - (*drv)[pd[0]] += - (0.0390625*(-(*drv)[cd[28]] - (*drv)[cd[29]] - (*drv)[cd[34]])); - (*drv)[pd[1]] += - (-0.0390625*(*drv)[cd[28]] - + 0.0234375*((*drv)[cd[29]] + (*drv)[cd[34]])); - (*drv)[pd[4]] += - (0.15625*(*drv)[cd[28]] - + 0.09375*((*drv)[cd[29]] + (*drv)[cd[34]])); - (*drv)[pd[5]] += - (-0.234375*(*drv)[cd[28]] - + 0.046875*(-(*drv)[cd[29]] - (*drv)[cd[34]])); - (*drv)[pd[6]] += - (0.15625*(*drv)[cd[28]] - + 0.03125*(-(*drv)[cd[29]] - (*drv)[cd[34]])); - (*drv)[pd[7]] += - (0.0625*((*drv)[cd[28]] + (*drv)[cd[34]]) - + 0.125*(*drv)[cd[29]]); - (*drv)[pd[8]] += -0.125*(*drv)[cd[29]]; - (*drv)[pd[10]] += 0.0625*(*drv)[cd[34]]; - (*drv)[pd[13]] += - (0.3125*(*drv)[cd[28]] - 0.125*(*drv)[cd[29]] - - 0.0625*(*drv)[cd[34]]); - (*drv)[pd[14]] += 0.375*(*drv)[cd[29]]; - (*drv)[pd[16]] += -0.0625*(*drv)[cd[34]]; - (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; - (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; - (*drv)[pd[28]] += -0.1875*(*drv)[cd[34]]; - (*drv)[pd[29]] += 0.1875*(*drv)[cd[34]]; - (*drv)[pd[31]] += - (-0.3125*(*drv)[cd[28]] - 0.375*(*drv)[cd[29]] - - 0.1875*(*drv)[cd[34]]); - (*drv)[pd[32]] += - ((*drv)[cd[30]] + 0.9375*(*drv)[cd[28]] + 0.375*(*drv)[cd[29]] - + 0.1875*(*drv)[cd[34]]); - (*drv)[pd[33]] += 0.75*(*drv)[cd[29]]; - (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; - break; - case 2: - (*drv)[pd[0]] += - (0.0390625*(-(*drv)[cd[25]] - (*drv)[cd[26]] - (*drv)[cd[34]])); - (*drv)[pd[1]] += (-0.0390625*(*drv)[cd[25]] - + 0.0234375*((*drv)[cd[26]] + (*drv)[cd[34]])); - (*drv)[pd[4]] += (0.15625*(*drv)[cd[25]] - + 0.09375*((*drv)[cd[26]] + (*drv)[cd[34]])); - (*drv)[pd[5]] += (-0.234375*(*drv)[cd[25]] - + 0.046875*(-(*drv)[cd[26]] - (*drv)[cd[34]])); - (*drv)[pd[6]] += (0.15625*(*drv)[cd[25]] - + 0.03125*(-(*drv)[cd[26]] - (*drv)[cd[34]])); - (*drv)[pd[7]] += 0.0625*(*drv)[cd[34]]; - (*drv)[pd[10]] += (0.0625*((*drv)[cd[25]] + (*drv)[cd[34]]) - + 0.125*(*drv)[cd[26]]); - (*drv)[pd[11]] += -0.125*(*drv)[cd[26]]; - (*drv)[pd[13]] += -0.0625*(*drv)[cd[34]]; - (*drv)[pd[16]] += (0.3125*(*drv)[cd[25]] - 0.125*(*drv)[cd[26]] - - 0.0625*(*drv)[cd[34]]); - (*drv)[pd[17]] += 0.375*(*drv)[cd[26]]; - (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; - (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; - (*drv)[pd[28]] += (-0.3125*(*drv)[cd[25]] - 0.375*(*drv)[cd[26]] - - 0.1875*(*drv)[cd[34]]); - (*drv)[pd[29]] += ((*drv)[cd[27]] + 0.9375*(*drv)[cd[25]] - + 0.375*(*drv)[cd[26]] + 0.1875*(*drv)[cd[34]]); - (*drv)[pd[30]] += 0.75*(*drv)[cd[26]]; - (*drv)[pd[31]] += -0.1875*(*drv)[cd[34]]; - (*drv)[pd[32]] += 0.1875*(*drv)[cd[34]]; - (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; - break; - case 3: - (*drv)[pd[0]] += -0.0390625*(*drv)[cd[34]]; - (*drv)[pd[1]] += 0.0234375*(*drv)[cd[34]]; - (*drv)[pd[4]] += 0.09375*(*drv)[cd[34]]; - (*drv)[pd[5]] += -0.046875*(*drv)[cd[34]]; - (*drv)[pd[6]] += -0.03125*(*drv)[cd[34]]; - (*drv)[pd[7]] += 0.0625*(*drv)[cd[34]]; - (*drv)[pd[10]] += 0.0625*(*drv)[cd[34]]; - (*drv)[pd[13]] += -0.0625*(*drv)[cd[34]]; - (*drv)[pd[16]] += -0.0625*(*drv)[cd[34]]; - (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; - (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; - (*drv)[pd[28]] += -0.1875*(*drv)[cd[34]]; - (*drv)[pd[29]] += 0.1875*(*drv)[cd[34]]; - (*drv)[pd[31]] += -0.1875*(*drv)[cd[34]]; - (*drv)[pd[32]] += 0.1875*(*drv)[cd[34]]; - (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; - break; - } - } - } - } - - // ====== coarseInter functions ====== - - void Lagrange::coarseInter0(DOFIndexed<double> *drv, RCNeighbourList* list, - int n, BasisFunction* basFct) - { - FUNCNAME("Lagrange::coarseInter0"); - DegreeOfFreedom cdof, pdof; - const DOFAdmin *admin; - const Mesh *mesh = NULL; - - if (n < 1) return; - - Element* el = list->getElement(0); - - if (!drv->getFESpace()) { - ERROR("no fe_space in dof_real_vec\n"); - return; - } else if (!drv->getFESpace()->getBasisFcts()) { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - - admin = drv->getFESpace()->getAdmin(); - mesh = drv->getFESpace()->getMesh(); - - /****************************************************************************/ - /* values on child[0] */ - /****************************************************************************/ - - cdof = el->getChild(0)->getDOF(mesh->getNode(CENTER)+2, - admin->getNumberOfPreDOFs(CENTER)); - pdof = el->getDOF(mesh->getNode(CENTER)+2, admin->getNumberOfPreDOFs(CENTER)); - (*drv)[pdof] = (*drv)[cdof]; - } - - void Lagrange::coarseInter2_1d(DOFIndexed<double> *drv, RCNeighbourList *list, - int n, BasisFunction* basFct) - { - FUNCNAME("real_coarseInter2"); - Element *el; - //double *v = NULL; - DegreeOfFreedom cdof, pdof; - const DOFAdmin *admin; - Mesh *mesh = NULL; + for (i = 1; i < n; i++) { + el = list->getElement(i); + typ = list->getType(i); + basFct->getLocalIndices(el, admin, pd); - if (n < 1) return; - el = list->getElement(0); + lr_set = 0; + if (list->getNeighbourElement(i,0) && list->getNeighbourNr(i,0) < i) + lr_set = 1; - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } + if (list->getNeighbourElement(i,1) && list->getNeighbourNr(i,1) < i) + lr_set += 2; - admin = drv->getFESpace()->getAdmin(); - mesh = const_cast<Mesh*>(drv->getFESpace()->getMesh()); + TEST_EXIT_DBG(lr_set)("no values set on both neighbours\n"); - /****************************************************************************/ - /* values on child[0] */ - /****************************************************************************/ + /****************************************************************************/ + /* values on child[0] */ + /****************************************************************************/ - cdof = el->getChild(0)->getDOF(mesh->getNode(VERTEX)+1, - admin->getNumberOfPreDOFs(VERTEX)); - pdof = el->getDOF(mesh->getNode(CENTER), admin->getNumberOfPreDOFs(CENTER)); - (*drv)[pdof] = (*drv)[cdof]; - } + cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); - void Lagrange::coarseInter2_2d(DOFIndexed<double> *drv, RCNeighbourList* list, - int n, BasisFunction* basFct) - { - FUNCNAME("Lagrange::coarseInter2_2d"); - Element *el; - DegreeOfFreedom cdof, pdof; - const DOFAdmin *admin; - const Mesh *mesh = NULL; + switch (lr_set) { + case 1: + (*drv)[pd[0]] += + (0.0390625*(-(*drv)[cd[16]] - (*drv)[cd[22]] + - (*drv)[cd[23]] - (*drv)[cd[28]]) + + 0.0234375*((*drv)[cd[18]] + (*drv)[cd[29]] + (*drv)[cd[34]])); + (*drv)[pd[1]] += + (0.0390625*(-(*drv)[cd[16]] - (*drv)[cd[22]] - (*drv)[cd[23]] + - (*drv)[cd[28]] - (*drv)[cd[29]] - (*drv)[cd[34]]) + + 0.0234375*(*drv)[cd[18]]); + (*drv)[pd[4]] += + (0.03125*((*drv)[cd[16]] + (*drv)[cd[22]] + (*drv)[cd[23]] + - (*drv)[cd[29]] - (*drv)[cd[34]]) + - 0.09375*(*drv)[cd[18]] + 0.15625*(*drv)[cd[28]]); + (*drv)[pd[5]] += + (0.015625*((*drv)[cd[16]] + (*drv)[cd[22]] + (*drv)[cd[23]]) + + 0.140625*(*drv)[cd[18]] - 0.234375*(*drv)[cd[28]] + + 0.046875*(-(*drv)[cd[29]] - (*drv)[cd[34]])); + (*drv)[pd[6]] += + (0.03125*((*drv)[cd[16]] + (*drv)[cd[22]] + (*drv)[cd[23]]) + + 0.09375*(-(*drv)[cd[18]] + (*drv)[cd[29]] + (*drv)[cd[34]]) + + 0.15625*(*drv)[cd[28]]); + (*drv)[pd[7]] += + (0.1875*(*drv)[cd[16]] + + 0.0625*(-(*drv)[cd[18]] + (*drv)[cd[23]] - (*drv)[cd[34]]) + + 0.125*((*drv)[cd[22]] - (*drv)[cd[29]]) + + 0.3125*(*drv)[cd[28]]); + (*drv)[pd[8]] += + (0.375*(-(*drv)[cd[16]] + (*drv)[cd[29]]) + - 0.125*(*drv)[cd[22]]); + (*drv)[pd[9]] += 0.5*(*drv)[cd[16]]; + (*drv)[pd[10]] += + (0.0625*((*drv)[cd[22]] - (*drv)[cd[34]]) + 0.125*(*drv)[cd[23]]); + (*drv)[pd[11]] += -0.125*(*drv)[cd[23]]; + (*drv)[pd[13]] += + (0.1875*(*drv)[cd[16]] + + 0.0625*(-(*drv)[cd[18]] + (*drv)[cd[23]] + + (*drv)[cd[28]] + (*drv)[cd[34]]) + + 0.125*((*drv)[cd[22]] + (*drv)[cd[29]])); + (*drv)[pd[14]] += + (-0.375*(*drv)[cd[16]] + + 0.125*(-(*drv)[cd[22]] - (*drv)[cd[29]])); + (*drv)[pd[15]] += 0.5*(*drv)[cd[16]]; + (*drv)[pd[16]] += + (0.0625*((*drv)[cd[22]] + (*drv)[cd[34]]) + 0.125*(*drv)[cd[23]]); + (*drv)[pd[17]] += -0.125*(*drv)[cd[23]]; + (*drv)[pd[22]] += + (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) - 0.125*(*drv)[cd[34]]); + (*drv)[pd[23]] += 0.5*(*drv)[cd[22]]; + (*drv)[pd[24]] += 0.5*(*drv)[cd[23]]; + (*drv)[pd[25]] += + (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) + 0.375*(*drv)[cd[34]]); + (*drv)[pd[26]] += 0.5*(*drv)[cd[22]]; + (*drv)[pd[27]] += 0.5*(*drv)[cd[23]]; + (*drv)[pd[28]] += + (-0.0625*(*drv)[cd[22]] - 0.125*(*drv)[cd[23]] + + 0.1875*(*drv)[cd[34]]); + (*drv)[pd[29]] += + (-0.0625*(*drv)[cd[22]] - 0.125*(*drv)[cd[23]] + - 0.1875*(*drv)[cd[34]]); + (*drv)[pd[30]] += 0.25*(*drv)[cd[23]]; + (*drv)[pd[31]] = + ((*drv)[cd[30]] + 0.1875*(-(*drv)[cd[16]] + (*drv)[cd[34]]) + + 0.5625*(*drv)[cd[18]] - 0.125*(*drv)[cd[22]] + - 0.0625*(*drv)[cd[23]] + 0.9375*(*drv)[cd[28]] + + 0.375*(*drv)[cd[29]]); + (*drv)[pd[32]] = + (0.1875*(-(*drv)[cd[16]] - (*drv)[cd[34]]) + + 0.5625*(*drv)[cd[18]] - 0.125*(*drv)[cd[22]] + - 0.0625*(*drv)[cd[23]] - 0.3125*(*drv)[cd[28]] + - 0.375*(*drv)[cd[29]]); + (*drv)[pd[33]] = + ((*drv)[cd[17]] + 0.75*((*drv)[cd[16]] + (*drv)[cd[29]]) + + 0.25*(*drv)[cd[22]]); + (*drv)[pd[34]] = + ((*drv)[cd[24]] + 0.5*((*drv)[cd[22]] + (*drv)[cd[23]]) + + 0.75*(*drv)[cd[34]]); + break; + case 2: + (*drv)[pd[0]] += + (0.0390625*(-(*drv)[cd[19]] - (*drv)[cd[22]] + - (*drv)[cd[23]] - (*drv)[cd[25]]) + + 0.0234375*((*drv)[cd[21]] + (*drv)[cd[26]] + (*drv)[cd[34]])); + (*drv)[pd[1]] += + (0.0390625*(-(*drv)[cd[19]] - (*drv)[cd[22]] - (*drv)[cd[23]] + - (*drv)[cd[25]] - (*drv)[cd[26]] - (*drv)[cd[34]]) + + 0.0234375*(*drv)[cd[21]]); + (*drv)[pd[4]] += + (0.03125*((*drv)[cd[19]] + (*drv)[cd[22]] + (*drv)[cd[23]] + - (*drv)[cd[26]] - (*drv)[cd[34]]) + - 0.09375*(*drv)[cd[21]] + 0.15625*(*drv)[cd[25]]); + (*drv)[pd[5]] += + (0.015625*((*drv)[cd[19]] + (*drv)[cd[22]] + (*drv)[cd[23]]) + + 0.140625*(*drv)[cd[21]] - 0.234375*(*drv)[cd[25]] + + 0.046875*(-(*drv)[cd[26]] - (*drv)[cd[34]])); + (*drv)[pd[6]] += + (0.03125*((*drv)[cd[19]] + (*drv)[cd[22]] + (*drv)[cd[23]]) + + 0.09375*(-(*drv)[cd[21]] + (*drv)[cd[26]] + (*drv)[cd[34]]) + + 0.15625*(*drv)[cd[25]]); + (*drv)[pd[7]] += + (0.125*(*drv)[cd[22]] + 0.0625*((*drv)[cd[23]] - (*drv)[cd[34]])); + (*drv)[pd[8]] += -0.125*(*drv)[cd[22]]; + (*drv)[pd[10]] += + (0.1875*(*drv)[cd[19]] + + 0.0625*(-(*drv)[cd[21]] + (*drv)[cd[22]] - (*drv)[cd[34]]) + + 0.125*((*drv)[cd[23]] - (*drv)[cd[26]]) + + 0.3125*(*drv)[cd[25]]); + (*drv)[pd[11]] += + (0.375*(-(*drv)[cd[19]] + (*drv)[cd[26]]) - 0.125*(*drv)[cd[23]]); + (*drv)[pd[12]] += 0.5*(*drv)[cd[19]]; + (*drv)[pd[13]] += + (0.125*(*drv)[cd[22]] + 0.0625*((*drv)[cd[23]] + (*drv)[cd[34]])); + (*drv)[pd[14]] += -0.125*(*drv)[cd[22]]; + (*drv)[pd[16]] += + (0.1875*(*drv)[cd[19]] + + 0.0625*(-(*drv)[cd[21]] + (*drv)[cd[22]] + + (*drv)[cd[25]] + (*drv)[cd[34]]) + + 0.125*((*drv)[cd[23]] + (*drv)[cd[26]])); + (*drv)[pd[17]] += + (-0.375*(*drv)[cd[19]] + + 0.125*(-(*drv)[cd[23]] - (*drv)[cd[26]])); + (*drv)[pd[18]] += 0.5*(*drv)[cd[19]]; + (*drv)[pd[22]] += + (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) - 0.125*(*drv)[cd[34]]); + (*drv)[pd[23]] += (0.5*(*drv)[cd[22]]); + (*drv)[pd[24]] += (0.5*(*drv)[cd[23]]); + (*drv)[pd[25]] += + (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) + 0.375*(*drv)[cd[34]]); + (*drv)[pd[26]] += (0.5*(*drv)[cd[22]]); + (*drv)[pd[27]] += (0.5*(*drv)[cd[23]]); + (*drv)[pd[28]] = + ((*drv)[cd[27]] + 0.1875*(-(*drv)[cd[19]] + (*drv)[cd[34]]) + + 0.5625*(*drv)[cd[21]] - 0.0625*(*drv)[cd[22]] + - 0.125*(*drv)[cd[23]] + 0.9375*(*drv)[cd[25]] + + 0.375*(*drv)[cd[26]]); + (*drv)[pd[29]] = + (0.1875*(-(*drv)[cd[19]] - (*drv)[cd[34]]) + + 0.5625*(*drv)[cd[21]] - 0.0625*(*drv)[cd[22]] + - 0.125*(*drv)[cd[23]] - 0.3125*(*drv)[cd[25]] + - 0.375*(*drv)[cd[26]]); + (*drv)[pd[30]] = + ((*drv)[cd[20]] + 0.75*((*drv)[cd[19]] + (*drv)[cd[26]]) + + 0.25*(*drv)[cd[23]]); + (*drv)[pd[31]] += + (-0.125*(*drv)[cd[22]] - 0.0625*(*drv)[cd[23]] + + 0.1875*(*drv)[cd[34]]); + (*drv)[pd[32]] += + (-0.125*(*drv)[cd[22]] - 0.0625*(*drv)[cd[23]] + - 0.1875*(*drv)[cd[34]]); + (*drv)[pd[33]] += 0.25*(*drv)[cd[22]]; + (*drv)[pd[34]] = + ((*drv)[cd[24]] + 0.5*((*drv)[cd[22]] + (*drv)[cd[23]]) + + 0.75*(*drv)[cd[34]]); + break; + case 3: + (*drv)[pd[0]] += + (0.0390625*(-(*drv)[cd[22]] - (*drv)[cd[23]]) + + 0.0234375*(*drv)[cd[34]]); + (*drv)[pd[1]] += + (0.0390625*(-(*drv)[cd[22]] - (*drv)[cd[23]] - (*drv)[cd[34]])); + (*drv)[pd[4]] += + (0.03125*((*drv)[cd[22]] + (*drv)[cd[23]] - (*drv)[cd[34]])); + (*drv)[pd[5]] += + (0.015625*((*drv)[cd[22]] + (*drv)[cd[23]]) + - 0.046875*(*drv)[cd[34]]); + (*drv)[pd[6]] += + (0.03125*((*drv)[cd[22]] + (*drv)[cd[23]]) + + 0.09375*(*drv)[cd[34]]); + (*drv)[pd[7]] += + (0.125*(*drv)[cd[22]] + 0.0625*((*drv)[cd[23]] - (*drv)[cd[34]])); + (*drv)[pd[8]] += -0.125*(*drv)[cd[22]]; + (*drv)[pd[10]] += + (0.0625*((*drv)[cd[22]] - (*drv)[cd[34]]) + 0.125*(*drv)[cd[23]]); + (*drv)[pd[11]] += -0.125*(*drv)[cd[23]]; + (*drv)[pd[13]] += + (0.125*(*drv)[cd[22]] + 0.0625*((*drv)[cd[23]] + (*drv)[cd[34]])); + (*drv)[pd[14]] += -0.125*(*drv)[cd[22]]; + (*drv)[pd[16]] += + (0.0625*((*drv)[cd[22]] + (*drv)[cd[34]]) + 0.125*(*drv)[cd[23]]); + (*drv)[pd[17]] += -0.125*(*drv)[cd[23]]; + (*drv)[pd[22]] += + (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) - 0.125*(*drv)[cd[34]]); + (*drv)[pd[23]] += 0.5*(*drv)[cd[22]]; + (*drv)[pd[24]] += 0.5*(*drv)[cd[23]]; + (*drv)[pd[25]] += + (0.25*(-(*drv)[cd[22]] - (*drv)[cd[23]]) + 0.375*(*drv)[cd[34]]); + (*drv)[pd[26]] += 0.5*(*drv)[cd[22]]; + (*drv)[pd[27]] += 0.5*(*drv)[cd[23]]; + (*drv)[pd[28]] += + (-0.0625*(*drv)[cd[22]] - 0.125*(*drv)[cd[23]] + + 0.1875*(*drv)[cd[34]]); + (*drv)[pd[29]] += + (-0.0625*(*drv)[cd[22]] - 0.125*(*drv)[cd[23]] + - 0.1875*(*drv)[cd[34]]); + (*drv)[pd[30]] += 0.25*(*drv)[cd[23]]; + (*drv)[pd[31]] += + (-0.125*(*drv)[cd[22]] - 0.0625*(*drv)[cd[23]] + + 0.1875*(*drv)[cd[34]]); + (*drv)[pd[32]] += + (-0.125*(*drv)[cd[22]] - 0.0625*(*drv)[cd[23]] + - 0.1875*(*drv)[cd[34]]); + (*drv)[pd[33]] += 0.25*(*drv)[cd[22]]; + (*drv)[pd[34]] = + ((*drv)[cd[24]] + 0.5*((*drv)[cd[22]] + (*drv)[cd[23]]) + + 0.75*(*drv)[cd[34]]); + break; + } + + /****************************************************************************/ + /* values on child[1] */ + /****************************************************************************/ - if (n < 1) return; - el = list->getElement(0); + cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; + if (typ == 0) { + switch (lr_set) { + case 1: + (*drv)[pd[0]] += + (0.0390625*(-(*drv)[cd[25]] - (*drv)[cd[26]] - (*drv)[cd[34]])); + (*drv)[pd[1]] += + (-0.0390625*(*drv)[cd[25]] + + 0.0234375*((*drv)[cd[26]] + (*drv)[cd[34]])); + (*drv)[pd[4]] += + (0.15625*(*drv)[cd[25]] + + 0.09375*((*drv)[cd[26]] + (*drv)[cd[34]])); + (*drv)[pd[5]] += + (-0.234375*(*drv)[cd[25]] + + 0.046875*(-(*drv)[cd[26]] - (*drv)[cd[34]])); + (*drv)[pd[6]] += + (0.15625*(*drv)[cd[25]] + + 0.03125*(-(*drv)[cd[26]] - (*drv)[cd[34]])); + (*drv)[pd[7]] += + (0.0625*((*drv)[cd[25]] + (*drv)[cd[34]]) + + 0.125*(*drv)[cd[26]]); + (*drv)[pd[8]] += -0.125*(*drv)[cd[26]]; + (*drv)[pd[10]] += 0.0625*(*drv)[cd[34]]; + (*drv)[pd[13]] += + (0.3125*(*drv)[cd[25]] - 0.125*(*drv)[cd[26]] + - 0.0625*(*drv)[cd[34]]); + (*drv)[pd[14]] += 0.375*(*drv)[cd[26]]; + (*drv)[pd[16]] += -0.0625*(*drv)[cd[34]]; + (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; + (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; + (*drv)[pd[28]] += -0.1875*(*drv)[cd[34]]; + (*drv)[pd[29]] += 0.1875*(*drv)[cd[34]]; + (*drv)[pd[31]] += + (-0.3125*(*drv)[cd[25]] - 0.375*(*drv)[cd[26]] + - 0.1875*(*drv)[cd[34]]); + (*drv)[pd[32]] += + ((*drv)[cd[27]] + 0.9375*(*drv)[cd[25]] + + 0.375*(*drv)[cd[26]] + 0.1875*(*drv)[cd[34]]); + (*drv)[pd[33]] += 0.75*(*drv)[cd[26]]; + (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; + break; + case 2: + (*drv)[pd[0]] += + (0.0390625*(-(*drv)[cd[28]] - (*drv)[cd[29]] - (*drv)[cd[34]])); + (*drv)[pd[1]] += + (-0.0390625*(*drv)[cd[28]] + + 0.0234375*((*drv)[cd[29]] + (*drv)[cd[34]])); + (*drv)[pd[4]] += + (0.15625*(*drv)[cd[28]] + + 0.09375*((*drv)[cd[29]] + (*drv)[cd[34]])); + (*drv)[pd[5]] += + (-0.234375*(*drv)[cd[28]] + + 0.046875*(-(*drv)[cd[29]] - (*drv)[cd[34]])); + (*drv)[pd[6]] += + (0.15625*(*drv)[cd[28]] + + 0.03125*(-(*drv)[cd[29]] - (*drv)[cd[34]])); + (*drv)[pd[7]] += 0.0625*(*drv)[cd[34]]; + (*drv)[pd[10]] += + (0.0625*((*drv)[cd[28]] + (*drv)[cd[34]]) + + 0.125*(*drv)[cd[29]]); + (*drv)[pd[11]] += -0.125*(*drv)[cd[29]]; + (*drv)[pd[13]] += -0.0625*(*drv)[cd[34]]; + (*drv)[pd[16]] += + (0.3125*(*drv)[cd[28]] - 0.125*(*drv)[cd[29]] + - 0.0625*(*drv)[cd[34]]); + (*drv)[pd[17]] += 0.375*(*drv)[cd[29]]; + (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; + (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; + (*drv)[pd[28]] += + (-0.3125*(*drv)[cd[28]] - 0.375*(*drv)[cd[29]] + - 0.1875*(*drv)[cd[34]]); + (*drv)[pd[29]] += + ((*drv)[cd[30]] + 0.9375*(*drv)[cd[28]] + 0.375*(*drv)[cd[29]] + + 0.1875*(*drv)[cd[34]]); + (*drv)[pd[30]] += 0.75*(*drv)[cd[29]]; + (*drv)[pd[31]] += -0.1875*(*drv)[cd[34]]; + (*drv)[pd[32]] += 0.1875*(*drv)[cd[34]]; + (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; + break; + case 3: + (*drv)[pd[0]] += -0.0390625*(*drv)[cd[34]]; + (*drv)[pd[1]] += 0.0234375*(*drv)[cd[34]]; + (*drv)[pd[4]] += 0.09375*(*drv)[cd[34]]; + (*drv)[pd[5]] += -0.046875*(*drv)[cd[34]]; + (*drv)[pd[6]] += -0.03125*(*drv)[cd[34]]; + (*drv)[pd[7]] += 0.0625*(*drv)[cd[34]]; + (*drv)[pd[10]] += 0.0625*(*drv)[cd[34]]; + (*drv)[pd[13]] += -0.0625*(*drv)[cd[34]]; + (*drv)[pd[16]] += -0.0625*(*drv)[cd[34]]; + (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; + (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; + (*drv)[pd[28]] += -0.1875*(*drv)[cd[34]]; + (*drv)[pd[29]] += 0.1875*(*drv)[cd[34]]; + (*drv)[pd[31]] += -0.1875*(*drv)[cd[34]]; + (*drv)[pd[32]] += 0.1875*(*drv)[cd[34]]; + (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; + break; + } + } else { + switch (lr_set) { + case 1: + (*drv)[pd[0]] += + (0.0390625*(-(*drv)[cd[28]] - (*drv)[cd[29]] - (*drv)[cd[34]])); + (*drv)[pd[1]] += + (-0.0390625*(*drv)[cd[28]] + + 0.0234375*((*drv)[cd[29]] + (*drv)[cd[34]])); + (*drv)[pd[4]] += + (0.15625*(*drv)[cd[28]] + + 0.09375*((*drv)[cd[29]] + (*drv)[cd[34]])); + (*drv)[pd[5]] += + (-0.234375*(*drv)[cd[28]] + + 0.046875*(-(*drv)[cd[29]] - (*drv)[cd[34]])); + (*drv)[pd[6]] += + (0.15625*(*drv)[cd[28]] + + 0.03125*(-(*drv)[cd[29]] - (*drv)[cd[34]])); + (*drv)[pd[7]] += + (0.0625*((*drv)[cd[28]] + (*drv)[cd[34]]) + + 0.125*(*drv)[cd[29]]); + (*drv)[pd[8]] += -0.125*(*drv)[cd[29]]; + (*drv)[pd[10]] += 0.0625*(*drv)[cd[34]]; + (*drv)[pd[13]] += + (0.3125*(*drv)[cd[28]] - 0.125*(*drv)[cd[29]] + - 0.0625*(*drv)[cd[34]]); + (*drv)[pd[14]] += 0.375*(*drv)[cd[29]]; + (*drv)[pd[16]] += -0.0625*(*drv)[cd[34]]; + (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; + (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; + (*drv)[pd[28]] += -0.1875*(*drv)[cd[34]]; + (*drv)[pd[29]] += 0.1875*(*drv)[cd[34]]; + (*drv)[pd[31]] += + (-0.3125*(*drv)[cd[28]] - 0.375*(*drv)[cd[29]] + - 0.1875*(*drv)[cd[34]]); + (*drv)[pd[32]] += + ((*drv)[cd[30]] + 0.9375*(*drv)[cd[28]] + 0.375*(*drv)[cd[29]] + + 0.1875*(*drv)[cd[34]]); + (*drv)[pd[33]] += 0.75*(*drv)[cd[29]]; + (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; + break; + case 2: + (*drv)[pd[0]] += + (0.0390625*(-(*drv)[cd[25]] - (*drv)[cd[26]] - (*drv)[cd[34]])); + (*drv)[pd[1]] += (-0.0390625*(*drv)[cd[25]] + + 0.0234375*((*drv)[cd[26]] + (*drv)[cd[34]])); + (*drv)[pd[4]] += (0.15625*(*drv)[cd[25]] + + 0.09375*((*drv)[cd[26]] + (*drv)[cd[34]])); + (*drv)[pd[5]] += (-0.234375*(*drv)[cd[25]] + + 0.046875*(-(*drv)[cd[26]] - (*drv)[cd[34]])); + (*drv)[pd[6]] += (0.15625*(*drv)[cd[25]] + + 0.03125*(-(*drv)[cd[26]] - (*drv)[cd[34]])); + (*drv)[pd[7]] += 0.0625*(*drv)[cd[34]]; + (*drv)[pd[10]] += (0.0625*((*drv)[cd[25]] + (*drv)[cd[34]]) + + 0.125*(*drv)[cd[26]]); + (*drv)[pd[11]] += -0.125*(*drv)[cd[26]]; + (*drv)[pd[13]] += -0.0625*(*drv)[cd[34]]; + (*drv)[pd[16]] += (0.3125*(*drv)[cd[25]] - 0.125*(*drv)[cd[26]] + - 0.0625*(*drv)[cd[34]]); + (*drv)[pd[17]] += 0.375*(*drv)[cd[26]]; + (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; + (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; + (*drv)[pd[28]] += (-0.3125*(*drv)[cd[25]] - 0.375*(*drv)[cd[26]] + - 0.1875*(*drv)[cd[34]]); + (*drv)[pd[29]] += ((*drv)[cd[27]] + 0.9375*(*drv)[cd[25]] + + 0.375*(*drv)[cd[26]] + 0.1875*(*drv)[cd[34]]); + (*drv)[pd[30]] += 0.75*(*drv)[cd[26]]; + (*drv)[pd[31]] += -0.1875*(*drv)[cd[34]]; + (*drv)[pd[32]] += 0.1875*(*drv)[cd[34]]; + (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; + break; + case 3: + (*drv)[pd[0]] += -0.0390625*(*drv)[cd[34]]; + (*drv)[pd[1]] += 0.0234375*(*drv)[cd[34]]; + (*drv)[pd[4]] += 0.09375*(*drv)[cd[34]]; + (*drv)[pd[5]] += -0.046875*(*drv)[cd[34]]; + (*drv)[pd[6]] += -0.03125*(*drv)[cd[34]]; + (*drv)[pd[7]] += 0.0625*(*drv)[cd[34]]; + (*drv)[pd[10]] += 0.0625*(*drv)[cd[34]]; + (*drv)[pd[13]] += -0.0625*(*drv)[cd[34]]; + (*drv)[pd[16]] += -0.0625*(*drv)[cd[34]]; + (*drv)[pd[22]] += 0.375*(*drv)[cd[34]]; + (*drv)[pd[25]] += -0.125*(*drv)[cd[34]]; + (*drv)[pd[28]] += -0.1875*(*drv)[cd[34]]; + (*drv)[pd[29]] += 0.1875*(*drv)[cd[34]]; + (*drv)[pd[31]] += -0.1875*(*drv)[cd[34]]; + (*drv)[pd[32]] += 0.1875*(*drv)[cd[34]]; + (*drv)[pd[34]] += 0.75*(*drv)[cd[34]]; + break; + } } + } + } - admin = drv->getFESpace()->getAdmin(); - mesh = drv->getFESpace()->getMesh(); + // ====== coarseInter functions ====== - /****************************************************************************/ - /* values on child[0] */ - /****************************************************************************/ + void Lagrange::coarseInter0(DOFIndexed<double> *drv, RCNeighbourList* list, + int n, BasisFunction* basFct) + { + FUNCNAME("Lagrange::coarseInter0()"); + + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); + + if (n < 1) + return; - cdof = el->getChild(0)->getDOF(mesh->getNode(VERTEX)+2, - admin->getNumberOfPreDOFs(VERTEX)); - pdof = el->getDOF(mesh->getNode(EDGE)+2, admin->getNumberOfPreDOFs(EDGE)); + Element* el = list->getElement(0); + const DOFAdmin *admin = drv->getFESpace()->getAdmin(); + const Mesh *mesh = drv->getFESpace()->getMesh(); + + // values on child[0] + DegreeOfFreedom cdof = el->getChild(0)->getDOF(mesh->getNode(CENTER)+2, + admin->getNumberOfPreDOFs(CENTER)); + DegreeOfFreedom pdof = el->getDOF(mesh->getNode(CENTER)+2, admin->getNumberOfPreDOFs(CENTER)); + (*drv)[pdof] = (*drv)[cdof]; + } + + void Lagrange::coarseInter2_1d(DOFIndexed<double> *drv, RCNeighbourList *list, + int n, BasisFunction* basFct) + { + FUNCNAME("real_coarseInter2()"); + + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); + + if (n < 1) + return; + + Element *el = list->getElement(0); + const DOFAdmin *admin = drv->getFESpace()->getAdmin(); + Mesh *mesh = const_cast<Mesh*>(drv->getFESpace()->getMesh()); + + // values on child[0] + DegreeOfFreedom cdof = el->getChild(0)->getDOF(mesh->getNode(VERTEX)+1, + admin->getNumberOfPreDOFs(VERTEX)); + DegreeOfFreedom pdof = el->getDOF(mesh->getNode(CENTER), admin->getNumberOfPreDOFs(CENTER)); + (*drv)[pdof] = (*drv)[cdof]; + } + + void Lagrange::coarseInter2_2d(DOFIndexed<double> *drv, RCNeighbourList* list, + int n, BasisFunction* basFct) + { + FUNCNAME("Lagrange::coarseInter2_2d()"); + + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); + + if (n < 1) + return; + + Element *el = list->getElement(0); + const DOFAdmin *admin = drv->getFESpace()->getAdmin(); + const Mesh *mesh = drv->getFESpace()->getMesh(); + + // values on child[0] + DegreeOfFreedom cdof = el->getChild(0)->getDOF(mesh->getNode(VERTEX)+2, + admin->getNumberOfPreDOFs(VERTEX)); + DegreeOfFreedom pdof = el->getDOF(mesh->getNode(EDGE)+2, admin->getNumberOfPreDOFs(EDGE)); (*drv)[pdof] = (*drv)[cdof]; } void Lagrange::coarseInter2_3d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseInter2_3d"); - Element *el; - int cdof, pdof; - const DOFAdmin *admin; - const Mesh *mesh = NULL; + FUNCNAME("Lagrange::coarseInter2_3d()"); - if (n < 1) return; - el = list->getElement(0); + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - admin = drv->getFESpace()->getAdmin(); - mesh = drv->getFESpace()->getMesh(); + if (n < 1) + return; - cdof = el->getChild(0)->getDOF(mesh->getNode(VERTEX)+3, - admin->getNumberOfPreDOFs(VERTEX)); - pdof = el->getDOF(mesh->getNode(EDGE), admin->getNumberOfPreDOFs(EDGE)); + Element *el = list->getElement(0); + const DOFAdmin *admin = drv->getFESpace()->getAdmin(); + const Mesh *mesh = drv->getFESpace()->getMesh(); + + int cdof = el->getChild(0)->getDOF(mesh->getNode(VERTEX)+3, + admin->getNumberOfPreDOFs(VERTEX)); + int pdof = el->getDOF(mesh->getNode(EDGE), admin->getNumberOfPreDOFs(EDGE)); (*drv)[pdof] = (*drv)[cdof]; } void Lagrange::coarseInter3_1d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseInter3_1d"); - const Element *el, *child; - int cdof, pdof, node, n0; - const DOFAdmin *admin; - const Mesh *mesh = NULL; + FUNCNAME("Lagrange::coarseInter3_1d()"); - if (n < 1) return; - el = list->getElement(0); + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - admin = drv->getFESpace()->getAdmin(); - mesh = drv->getFESpace()->getMesh(); + if (n < 1) + return; - node = mesh->getNode(EDGE); - n0 = admin->getNumberOfPreDOFs(EDGE); + const Element *el = list->getElement(0); + const Element *child = el->getChild(0); + const DOFAdmin *admin = drv->getFESpace()->getAdmin(); + const Mesh *mesh = drv->getFESpace()->getMesh(); + int node = mesh->getNode(EDGE); + int n0 = admin->getNumberOfPreDOFs(EDGE); - /****************************************************************************/ - /* values on child[0] */ - /****************************************************************************/ + // values on child[0] - child = el->getChild(0); + int cdof, pdof; if (el->getDOF(0, 0) < el->getDOF(1, 0)) - pdof = el->getDOF(node+2, n0); + pdof = el->getDOF(node + 2, n0); else - pdof = el->getDOF(node+2, n0+1); + pdof = el->getDOF(node + 2, n0 + 1); if (child->getDOF(1, 0) < child->getDOF(2, 0)) - cdof = child->getDOF(node, n0+1); + cdof = child->getDOF(node, n0 + 1); else cdof = child->getDOF(node, n0); @@ -4495,9 +4375,7 @@ namespace AMDiS { (*drv)[el->getDOF(mesh->getNode(CENTER), admin->getNumberOfPreDOFs(CENTER))] = (*drv)[cdof]; - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ + // values on child[1] child = el->getChild(1); @@ -4513,11 +4391,10 @@ namespace AMDiS { (*drv)[pdof] = (*drv)[cdof]; - if (n <= 1) return; + if (n <= 1) + return; - /****************************************************************************/ - /* adjust neighbour values */ - /****************************************************************************/ + // adjust neighbour values el = list->getElement(1); child = el->getChild(0); @@ -4534,37 +4411,24 @@ namespace AMDiS { void Lagrange::coarseInter3_2d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseInter3_2d"); - const Element *el, *child; - int cdof, pdof, node, n0; - const DOFAdmin *admin; - const Mesh *mesh = NULL; + FUNCNAME("Lagrange::coarseInter3_2d()"); - if (n < 1) return; - el = list->getElement(0); + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - admin = drv->getFESpace()->getAdmin(); - mesh = drv->getFESpace()->getMesh(); + if (n < 1) + return; - node = mesh->getNode(EDGE); - n0 = admin->getNumberOfPreDOFs(EDGE); + const Element *el = list->getElement(0); + const Element *child = el->getChild(0); + const DOFAdmin *admin = drv->getFESpace()->getAdmin(); + const Mesh *mesh = drv->getFESpace()->getMesh(); + int node = mesh->getNode(EDGE); + int n0 = admin->getNumberOfPreDOFs(EDGE); - /****************************************************************************/ - /* values on child[0] */ - /****************************************************************************/ + // values on child[0] - child = el->getChild(0); + int cdof, pdof; if (el->getDOF(0, 0) < el->getDOF(1, 0)) pdof = el->getDOF(node+2, n0); @@ -4586,9 +4450,7 @@ namespace AMDiS { (*drv)[el->getDOF(mesh->getNode(CENTER), admin->getNumberOfPreDOFs(CENTER))] = (*drv)[cdof]; - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ + // values on child[1] child = el->getChild(1); @@ -4604,11 +4466,10 @@ namespace AMDiS { (*drv)[pdof] = (*drv)[cdof]; - if (n <= 1) return; + if (n <= 1) + return; - /****************************************************************************/ - /* adjust neighbour values */ - /****************************************************************************/ + // adjust neighbour values el = list->getElement(1); child = el->getChild(0); @@ -4625,43 +4486,29 @@ namespace AMDiS { void Lagrange::coarseInter3_3d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseInter3_3d"); - const Element *el; - int i, lr_set; - int node_e, node_f, n0_e, n0_f; - const DegreeOfFreedom **pds, **cds; - DegreeOfFreedom pd_o[20], cd, pd; - const DOFAdmin *admin; + FUNCNAME("Lagrange::coarseInter3_3d()"); - if (n < 1) return; - el = list->getElement(0); + if (n < 1) + return; - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); - admin = drv->getFESpace()->getAdmin(); + const Element *el = list->getElement(0); + const DOFAdmin *admin = drv->getFESpace()->getAdmin(); - node_e = drv->getFESpace()->getMesh()->getNode(EDGE); - n0_e = admin->getNumberOfPreDOFs(EDGE); - node_f = drv->getFESpace()->getMesh()->getNode(FACE); - n0_f = admin->getNumberOfPreDOFs(FACE); + int node_e = drv->getFESpace()->getMesh()->getNode(EDGE); + int n0_e = admin->getNumberOfPreDOFs(EDGE); + int node_f = drv->getFESpace()->getMesh()->getNode(FACE); + int n0_f = admin->getNumberOfPreDOFs(FACE); - /****************************************************************************/ - /* values on child[0] */ - /****************************************************************************/ - pds = el->getDOF(); - cds = el->getChild(0)->getDOF(); - basFct->getLocalIndices(el, admin, pd_o); + // values on child[0] + + const DegreeOfFreedom **pds = el->getDOF(); + const DegreeOfFreedom **cds = el->getChild(0)->getDOF(); + DegreeOfFreedom pd_o[20], cd, pd; + basFct->getLocalIndices(el, admin, pd_o); pd = (pds[0][0] < pds[1][0]) ? pds[node_e][n0_e] : pds[node_e][n0_e+1]; cd = cds[0][0] < cds[3][0] ? cds[node_e+2][n0_e+1] : cds[node_e+2][n0_e]; (*drv)[pd] = (*drv)[cd]; @@ -4674,124 +4521,96 @@ namespace AMDiS { cd = cds[1][0] < cds[3][0] ? cds[node_e+4][n0_e+1] : cds[node_e+4][n0_e]; (*drv)[pd] = (*drv)[cd]; - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ + // values on child[1] cds = el->getChild(1)->getDOF(); - pd = (pds[0][0] < pds[1][0]) ? pds[node_e][n0_e+1] : pds[node_e][n0_e]; cd = cds[0][0] < cds[3][0] ? cds[node_e+2][n0_e+1] : cds[node_e+2][n0_e]; (*drv)[pd] = (*drv)[cd]; - /****************************************************************************/ - /* adjust neighbour values */ - /****************************************************************************/ + // adjust neighbour values - for (i = 1; i < n; i++) - { - el = list->getElement(i); - - pds = el->getDOF(); - cds = el->getChild(0)->getDOF(); - - 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"); - - switch(lr_set) - { - case 1: - pd = el->getDOF(node_f+3, n0_f); - cd = cds[1][0] < cds[3][0] ? - cds[node_e+4][n0_e+1] : cds[node_e+4][n0_e]; - (*drv)[pd] = (*drv)[cd]; - break; - case 2: - pd = el->getDOF(node_f+2, n0_f); - cd = cds[2][0] < cds[3][0] ? - cds[node_e+5][n0_e+1] : cds[node_e+5][n0_e]; - (*drv)[pd] = (*drv)[cd]; - break; - } + for (int i = 1; i < n; i++) { + el = list->getElement(i); + + pds = el->getDOF(); + cds = el->getChild(0)->getDOF(); + + int 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"); + + switch (lr_set) { + case 1: + pd = el->getDOF(node_f+3, n0_f); + cd = cds[1][0] < cds[3][0] ? + cds[node_e+4][n0_e+1] : cds[node_e+4][n0_e]; + (*drv)[pd] = (*drv)[cd]; + break; + case 2: + pd = el->getDOF(node_f+2, n0_f); + cd = cds[2][0] < cds[3][0] ? + cds[node_e+5][n0_e+1] : cds[node_e+5][n0_e]; + (*drv)[pd] = (*drv)[cd]; + break; } + } } void Lagrange::coarseInter4_1d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseInter4_1d"); - Element *el; - DegreeOfFreedom pdof[5]; - const DegreeOfFreedom *cdof; - const DOFAdmin *admin; + FUNCNAME("Lagrange::coarseInter4_1d()"); - if (n < 1) return; - el = list->getElement(0); - - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); - admin = drv->getFESpace()->getAdmin(); + if (n < 1) + return; + Element *el = list->getElement(0); + const DOFAdmin *admin = drv->getFESpace()->getAdmin(); + DegreeOfFreedom pdof[5]; basFct->getLocalIndices(el, admin, pdof); - /****************************************************************************/ - /* values on child[0] */ - /****************************************************************************/ + // values on child[0] - cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL); + const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL); (*drv)[pdof[9]] = (*drv)[cdof[4]]; (*drv)[pdof[10]] = (*drv)[cdof[2]]; (*drv)[pdof[12]] = (*drv)[cdof[14]]; (*drv)[pdof[14]] = (*drv)[cdof[7]]; - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ + // values on child[1] cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL); (*drv)[pdof[11]] = (*drv)[cdof[7]]; (*drv)[pdof[13]] = (*drv)[cdof[14]]; - if (n <= 1) return; + if (n <= 1) + return; - /****************************************************************************/ - /* adjust neighbour values */ - /****************************************************************************/ + // adjust neighbour values el = list->getElement(1); basFct->getLocalIndices(el, admin, pdof); - /****************************************************************************/ - /* values on neighbour's child[0] */ - /****************************************************************************/ + // values on neighbour's child[0] cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL); (*drv)[pdof[12]] = (*drv)[cdof[14]]; (*drv)[pdof[14]] = (*drv)[cdof[7]]; - /****************************************************************************/ - /* values on neighbour's child[1] */ - /****************************************************************************/ + // values on neighbour's child[1] cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL); @@ -4801,73 +4620,52 @@ namespace AMDiS { void Lagrange::coarseInter4_2d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct) { - FUNCNAME("Lagrange::coarseInter4_2d"); - const Element *el; - DegreeOfFreedom pdof[15]; - const DegreeOfFreedom *cdof; - const DOFAdmin *admin; - - if (n < 1) return; - el = list->getElement(0); + FUNCNAME("Lagrange::coarseInter4_2d()"); - if (!drv->getFESpace()) - { - ERROR("no fe_space in dof_real_vec\n"); - return; - } - else if (!drv->getFESpace()->getBasisFcts()) - { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); - admin = drv->getFESpace()->getAdmin(); + if (n < 1) + return; + const Element *el = list->getElement(0); + const DOFAdmin *admin = drv->getFESpace()->getAdmin(); + DegreeOfFreedom pdof[15]; basFct->getLocalIndices(el, admin, pdof); - /****************************************************************************/ - /* values on child[0] */ - /****************************************************************************/ + // values on child[0] - cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL); + const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL); (*drv)[pdof[9]] = (*drv)[cdof[4]]; (*drv)[pdof[10]] = (*drv)[cdof[2]]; (*drv)[pdof[12]] = (*drv)[cdof[14]]; (*drv)[pdof[14]] = (*drv)[cdof[7]]; - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ + // values on child[1] cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL); (*drv)[pdof[11]] = (*drv)[cdof[7]]; (*drv)[pdof[13]] = (*drv)[cdof[14]]; - if (n <= 1) return; + if (n <= 1) + return; - /****************************************************************************/ - /* adjust neighbour values */ - /****************************************************************************/ + // adjust neighbour values el = list->getElement(1); basFct->getLocalIndices(el, admin, pdof); - /****************************************************************************/ - /* values on neighbour's child[0] */ - /****************************************************************************/ + // values on neighbour's child[0] cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL); (*drv)[pdof[12]] = (*drv)[cdof[14]]; (*drv)[pdof[14]] = (*drv)[cdof[7]]; - /****************************************************************************/ - /* values on neighbour's child[1] */ - /****************************************************************************/ + // values on neighbour's child[1] cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL); @@ -4879,35 +4677,22 @@ namespace AMDiS { { FUNCNAME("void Lagrange::coarseInter4_3d()"); - DegreeOfFreedom pd[35]; - const DegreeOfFreedom *cd; - const Element *el; - int i, typ, lr_set; - const DOFAdmin *admin; + TEST_EXIT(drv->getFESpace())("No fe_space in dof_real_vec!\n"); + TEST_EXIT(drv->getFESpace()->getBasisFcts())("No basis functions in fe_space!\n"); if (n < 1) return; - el = list->getElement(0); - typ = list->getType(0); - if (!drv->getFESpace()) { - ERROR("no fe_space in dof_real_vec\n"); - return; - } else if (!drv->getFESpace()->getBasisFcts()) { - ERROR("no basis functions in fe_space %s\n", - drv->getFESpace()->getName().c_str()); - return; - } - - admin = drv->getFESpace()->getAdmin(); + const Element* el = list->getElement(0); + int typ = list->getType(0); + const DOFAdmin *admin = drv->getFESpace()->getAdmin(); + DegreeOfFreedom pd[35]; basFct->getLocalIndices(el, admin, pd); - /****************************************************************************/ - /* values on child[0] */ - /****************************************************************************/ + // values on child[0] - cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); + const DegreeOfFreedom *cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); (*drv)[pd[4]] = (*drv)[cd[11]]; (*drv)[pd[5]] = (*drv)[cd[3]]; @@ -4917,41 +4702,32 @@ namespace AMDiS { (*drv)[pd[33]] = (*drv)[cd[17]]; (*drv)[pd[34]] = (*drv)[cd[24]]; - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ + // values on child[1] cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); if (typ == 0) { - - /****************************************************************************/ - /* parent of el_type 0 */ - /****************************************************************************/ + // parent of el_type 0 (*drv)[pd[6]] = (*drv)[cd[11]]; (*drv)[pd[29]] = (*drv)[cd[30]]; (*drv)[pd[32]] = (*drv)[cd[27]]; - } else { - /****************************************************************************/ - /* parent of el_type 1|2 */ - /****************************************************************************/ + } else { + // parent of el_type 1|2 (*drv)[pd[6]] = (*drv)[cd[11]]; (*drv)[pd[29]] = (*drv)[cd[27]]; (*drv)[pd[32]] = (*drv)[cd[30]]; } - /****************************************************************************/ - /* adjust neighbour values */ - /****************************************************************************/ + // adjust neighbour values - for (i = 1; i < n; i++) { + for (int i = 1; i < n; i++) { el = list->getElement(i); typ = list->getType(i); basFct->getLocalIndices(el, admin, pd); - lr_set = 0; + int lr_set = 0; if (list->getNeighbourElement(i,0) && list->getNeighbourNr(i,0) < i) lr_set = 1; @@ -4960,13 +4736,11 @@ namespace AMDiS { TEST_EXIT_DBG(lr_set)("no values set on both neighbours\n"); - /****************************************************************************/ - /* values on child[0] */ - /****************************************************************************/ + // values on child[0] cd = basFct->getLocalIndices(el->getChild(0), admin, NULL); - switch(lr_set) { + switch (lr_set) { case 1: (*drv)[pd[31]] = (*drv)[cd[30]]; (*drv)[pd[33]] = (*drv)[cd[17]]; @@ -4982,14 +4756,12 @@ namespace AMDiS { break; } - /****************************************************************************/ - /* values on child[1] */ - /****************************************************************************/ + // values on child[1] cd = basFct->getLocalIndices(el->getChild(1), admin, NULL); if (typ == 0) { - switch(lr_set) { + switch (lr_set) { case 1: (*drv)[pd[32]] = (*drv)[cd[27]]; break; @@ -4998,7 +4770,7 @@ namespace AMDiS { break; } } else { - switch(lr_set) { + switch (lr_set) { case 1: (*drv)[pd[32]] = (*drv)[cd[30]]; break; diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h index b0f23c249ac65e2a81444f4254240377880cb4b1..cc1d9e8261efb6b43b4c4c68b4abeecc08932027 100644 --- a/AMDiS/src/Lagrange.h +++ b/AMDiS/src/Lagrange.h @@ -119,14 +119,18 @@ namespace AMDiS { } /// Implements BasisFunction::getLocalIndices(). - const DegreeOfFreedom *getLocalIndices(const Element*, - const DOFAdmin*, - DegreeOfFreedom*) const; + const DegreeOfFreedom *getLocalIndices(const Element *el, + const DOFAdmin *admin, + DegreeOfFreedom *dofs) const; /// - void getLocalIndicesVec(const Element*, - const DOFAdmin*, - Vector<DegreeOfFreedom>*) const; + void getLocalIndicesVec(const Element *el, + const DOFAdmin *admin, + Vector<DegreeOfFreedom> *vec) const; + + void getLocalDofPtrVec(const Element *el, + const DOFAdmin *admin, + std::vector<const DegreeOfFreedom*>& vec) const; /// Implements BasisFunction::l2ScpFctBas void l2ScpFctBas(Quadrature* q, diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index ce15ab509e28ed9d81681c70277e6456b8bbd081..5a4cf2c1857122b03529737ae4f974c9a95bd90f 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -267,53 +267,54 @@ namespace AMDiS { me->setIndex(macroElements.size()); } - void Mesh::removeMacroElements(std::vector<MacroElement*>& macros) + void Mesh::removeMacroElements(std::vector<MacroElement*>& macros, + const FiniteElemSpace *feSpace) { FUNCNAME("Mesh::removeMacroElement()"); - TEST_EXIT(dim == 2)("Not yet implemented!\n"); + typedef std::map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap; + typedef std::map<const DegreeOfFreedom*, GeoIndex> DofPosMap; + + TEST_EXIT(dim == 2)("Not yet implemented for dim != 2!\n"); + TEST_EXIT(admin.size() == 1)("Not yet implemented for multiple admins!\n"); + TEST_EXIT(admin[0])("There is something wrong!\n"); + + ElementDofIterator elDofIter(feSpace); // Map that stores for each dof pointer (which may have a list of dofs) // all macro element indices that own the dof. - std::map<const DegreeOfFreedom*, std::set<MacroElement*> > dofsOwner; + DofElMap dofsOwner; + DofPosMap dofsPosIndex; // Determine all dof owner macro elements. for (std::deque<MacroElement*>::iterator macroIt = macroElements.begin(); macroIt != macroElements.end(); ++macroIt) { - Element *el = (*macroIt)->getElement(); - for (int i = 0; i < 3; i++) - dofsOwner[el->getDOF(i)].insert(*macroIt); + elDofIter.reset((*macroIt)->getElement()); + do { + dofsOwner[elDofIter.getDofPtr()].insert(*macroIt); + dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex(); + } while (elDofIter.nextStrict()); } - + // Remove all the given macro elements. for (std::vector<MacroElement*>::iterator macroIt = macros.begin(); macroIt != macros.end(); ++macroIt) { - bool found = false; - - // Remove the macro element from mesh's list of all macro elements. - for (std::deque<MacroElement*>::iterator it = macroElements.begin(); - it != macroElements.end(); - ++it) { - if (*it == *macroIt) { - macroElements.erase(it, it + 1); - found = true; - break; - } - } - - TEST_EXIT(found)("Cannot find MacroElement that should be removed!\n"); - + + std::deque<MacroElement*>::iterator mEl = + find(macroElements.begin(), macroElements.end(), *macroIt); + TEST_EXIT(mEl != macroElements.end()) + ("Cannot find MacroElement that should be removed!\n"); + macroElements.erase(mEl); + // Go through all neighbours of the macro element and remove this macro element // to be neighbour of some other macro element. for (int i = 0; i <= dim; i++) { if ((*macroIt)->getNeighbour(i)) { - for (int j = 0; j <= dim; j++) { - if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt) { + for (int j = 0; j <= dim; j++) + if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt) (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL); - } - } } else { // There is no neighbour at this edge, so we have to decrease the number // of edges in the mesh. @@ -325,30 +326,28 @@ namespace AMDiS { nElements--; // Remove this macro element from the dof owner list. - for (std::map<const DegreeOfFreedom*, std::set<MacroElement*> >::iterator dofsIt = dofsOwner.begin(); - dofsIt != dofsOwner.end(); - ++dofsIt) { + for (DofElMap::iterator dofsIt = dofsOwner.begin(); + dofsIt != dofsOwner.end(); ++dofsIt) { std::set<MacroElement*>::iterator mIt = dofsIt->second.find(*macroIt); - if (mIt != dofsIt->second.end()) { + if (mIt != dofsIt->second.end()) dofsIt->second.erase(mIt); - } } // And remove the macro element from memory delete *macroIt; } + int nRemainDofs = 0; // Check now all the dofs, that have no owner anymore and therefore have to // be removed. - for (std::map<const DegreeOfFreedom*, std::set<MacroElement*> >::iterator dofsIt = dofsOwner.begin(); - dofsIt != dofsOwner.end(); - ++dofsIt) { - if (dofsIt->second.size() == 0) { - freeDOF(const_cast<DegreeOfFreedom*>(dofsIt->first), VERTEX); - } else { + for (DofElMap::iterator dofsIt = dofsOwner.begin(); + dofsIt != dofsOwner.end(); ++dofsIt) { + if (dofsIt->second.size() == 0) + freeDOF(const_cast<DegreeOfFreedom*>(dofsIt->first), + dofsPosIndex[dofsIt->first]); + else nRemainDofs++; - } } nVertices = nRemainDofs; @@ -868,26 +867,27 @@ namespace AMDiS { { DimVec<double>* baryCoords; bool found = false; - ElementDofIterator elDofIter(feSpace->getAdmin()); TraverseStack stack; + Vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber()); + ElInfo *elInfo = stack.traverseFirst(this, -1, - Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); + Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); while (elInfo) { - elDofIter.reset(elInfo->getElement()); - int i = 0; - do { - if (elDofIter.getDofPtr() == dof) { + feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), + feSpace->getAdmin(), + &dofVec); + for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) { + if (dofVec[i] == *dof) { baryCoords = feSpace->getBasisFcts()->getCoords(i); elInfo->coordToWorld(*baryCoords, coords); found = true; - break; + break; } - i++; - } while (elDofIter.next()); - + } + if (found) break; - + elInfo = stack.traverseNext(elInfo); } diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h index 4fc0f3e71311ccf67fd91618f59151ee6db924d7..707eb55f38a4c2fbc04f176fe08bbf1357c16d29 100644 --- a/AMDiS/src/Mesh.h +++ b/AMDiS/src/Mesh.h @@ -429,7 +429,8 @@ namespace AMDiS { * that there are no global or local refinements, i.e., all macro elements have * no children. */ - void removeMacroElements(std::vector<MacroElement*>& macros); + void removeMacroElements(std::vector<MacroElement*>& macros, + const FiniteElemSpace* feSpace); /// Frees the array of DOF pointers (see \ref createDOFPtrs) void freeDOFPtrs(DegreeOfFreedom **ptrs); diff --git a/AMDiS/src/ParallelDomainProblem.cc b/AMDiS/src/ParallelDomainProblem.cc index 1427fb61bb22d395be01767bb375a509c7ea0eaa..d7c1b0bf91e4876385ef4e28501874a87ea36e97 100644 --- a/AMDiS/src/ParallelDomainProblem.cc +++ b/AMDiS/src/ParallelDomainProblem.cc @@ -85,11 +85,13 @@ namespace AMDiS { createInteriorBoundaryInfo(rankDOFs, boundaryDOFs); +#if (DEBUG != 0) + DbgTestInteriorBoundary(); +#endif // === Remove all macro elements that are not part of the rank partition. === removeMacroElements(); - // === Reset all DOFAdmins of the mesh. === @@ -107,6 +109,7 @@ namespace AMDiS { admin.setFirstHole(mapLocalGlobalDOFs.size()); } + // === Global refinements. === int globalRefinement = 0; @@ -118,11 +121,12 @@ namespace AMDiS { updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs); } - #if (DEBUG != 0) - testInteriorBoundary(); + DbgTestCommonDofs(); #endif + exit(0); + // === Create petsc matrix. === int ierr; @@ -408,20 +412,17 @@ namespace AMDiS { // === Once we have this information, we must care about their order. Eventually === - // === all the boundaries have to be in the same order on both ranks (because === - // === each boundary is shared by exactly two ranks). === + // === all the boundaries have to be in the same order on both ranks that share === + // === the bounday. === - std::vector<int*> sendBuffers; - std::vector<int*> recvBuffers; + std::vector<int*> sendBuffers, recvBuffers; MPI::Request request[myIntBoundary.boundary.size() + otherIntBoundary.boundary.size()]; int requestCounter = 0; - for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = - myIntBoundary.boundary.begin(); - rankIt != myIntBoundary.boundary.end(); - ++rankIt) { + for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin(); + rankIt != myIntBoundary.boundary.end(); ++rankIt) { int nSendInt = rankIt->second.size(); int* buffer = new int[nSendInt]; for (int i = 0; i < nSendInt; i++) @@ -432,14 +433,12 @@ namespace AMDiS { mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0); } - for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = - otherIntBoundary.boundary.begin(); - rankIt != otherIntBoundary.boundary.end(); - ++rankIt) { + for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); + rankIt != otherIntBoundary.boundary.end(); ++rankIt) { int nRecvInt = rankIt->second.size(); int *buffer = new int[nRecvInt]; recvBuffers.push_back(buffer); - + request[requestCounter++] = mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0); } @@ -449,10 +448,8 @@ namespace AMDiS { int i = 0; - for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = - otherIntBoundary.boundary.begin(); - rankIt != otherIntBoundary.boundary.end(); - ++rankIt) { + for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); + rankIt != otherIntBoundary.boundary.end(); ++rankIt) { // === We have received from rank "rankIt->first" the ordered list of element === // === indices. We now have to sort the corresponding list in this rank to === @@ -498,7 +495,7 @@ namespace AMDiS { macrosToRemove.push_back(*it); } - mesh->removeMacroElements(macrosToRemove); + mesh->removeMacroElements(macrosToRemove, feSpace); } @@ -509,7 +506,6 @@ namespace AMDiS { { FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()"); - // === Get rank information about DOFs. === // Stores to each DOF pointer the set of ranks the DOF is part of. @@ -541,8 +537,7 @@ namespace AMDiS { std::map<int, int> recvNewDofs; for (DofToRank::iterator it = boundaryDOFs.begin(); - it != boundaryDOFs.end(); - ++it) { + it != boundaryDOFs.end(); ++it) { if (it->second == mpiRank) { // If the boundary dof is a rank dof, it must be send to other ranks. @@ -575,8 +570,7 @@ namespace AMDiS { // === Send and receive the dof indices at boundary. === - std::vector<int*> sendBuffers(sendNewDofs.size()); - std::vector<int*> recvBuffers(recvNewDofs.size()); + std::vector<int*> sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size()); sendDofs.clear(); recvDofs.clear(); @@ -597,7 +591,7 @@ namespace AMDiS { dofIt = sendIt->second.begin(); dofIt != sendIt->second.end(); ++dofIt) { - sendBuffers[i][c++] = (dofIt->first)[0]; + sendBuffers[i][c++] = *(dofIt->first); sendBuffers[i][c++] = dofIt->second; sendDofs[sendIt->first].push_back(dofIt->first); @@ -639,13 +633,12 @@ namespace AMDiS { isRankDOF.clear(); for (int i = 0; i < nRankDOFs; i++) { - const_cast<DegreeOfFreedom*>(rankDOFs[i])[0] = i; + *const_cast<DegreeOfFreedom*>(rankDOFs[i]) = i; mapLocalGlobalDOFs[i] = rstart + i; mapGlobalLocalDOFs[rstart + i] = i; isRankDOF[i] = true; } - // === Change dof indices at boundary from other ranks. === // Within this small data structure we track which dof index was already changed. @@ -656,11 +649,9 @@ namespace AMDiS { // rule was applied, the dof pointer is set to false in this data structure and // is not allowed to be changed anymore. std::map<const DegreeOfFreedom*, bool> dofChanged; - for (DofToRank::iterator dofIt = boundaryDOFs.begin(); - dofIt != boundaryDOFs.end(); - ++dofIt) { + for (DofToRank::iterator dofIt = boundaryDOFs.begin(); dofIt != boundaryDOFs.end(); + ++dofIt) dofChanged[dofIt->first] = false; - } i = 0; for (std::map<int, int>::iterator recvIt = recvNewDofs.begin(); @@ -677,9 +668,8 @@ namespace AMDiS { // Iterate over all boundary dofs to find the dof, which index we have to change. - for (DofToRank::iterator dofIt = boundaryDOFs.begin(); - dofIt != boundaryDOFs.end(); - ++dofIt) { + for (DofToRank::iterator dofIt = boundaryDOFs.begin(); + dofIt != boundaryDOFs.end(); ++dofIt) { if (*(dofIt->first) == oldDof && !dofChanged[dofIt->first]) { dofChanged[dofIt->first] = true; @@ -703,24 +693,22 @@ namespace AMDiS { } - void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, - int& nOverallDOFs) + void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs) { FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()"); std::set<const DegreeOfFreedom*> rankDOFs; DofToRank newBoundaryDOFs; - + // === Get all DOFs in ranks partition. === - ElementDofIterator elDofIt(feSpace->getAdmin()); + ElementDofIterator elDofIt(feSpace); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { - Element *element = elInfo->getElement(); - + Element *element = elInfo->getElement(); elDofIt.reset(element); do { rankDOFs.insert(elDofIt.getDofPtr()); @@ -735,13 +723,12 @@ namespace AMDiS { RankToDofContainer sendNewDofs; RankToDofContainer recvNewDofs; - for (std::map<int, std::vector<AtomicBoundary> >::iterator it = - myIntBoundary.boundary.begin(); - it != myIntBoundary.boundary.end(); - ++it) { + for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin(); + it != myIntBoundary.boundary.end(); ++it) { + for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin(); - boundIt != it->second.end(); - ++boundIt) { + boundIt != it->second.end(); ++boundIt) { + const DegreeOfFreedom *dof1 = NULL; const DegreeOfFreedom *dof2 = NULL; @@ -769,36 +756,37 @@ namespace AMDiS { newBoundaryDOFs[dof1] = boundaryDOFs[dof1]; newBoundaryDOFs[dof2] = boundaryDOFs[dof2]; - - if (find(sendNewDofs[it->first].begin(), sendNewDofs[it->first].end(), dof1) == - sendNewDofs[it->first].end()) - sendNewDofs[it->first].push_back(dof1); - if (find(sendNewDofs[it->first].begin(), sendNewDofs[it->first].end(), dof2) == - sendNewDofs[it->first].end()) - sendNewDofs[it->first].push_back(dof2); - + + if (find(sendNewDofs[it->first].begin(), sendNewDofs[it->first].end(), dof1) == + sendNewDofs[it->first].end()) + sendNewDofs[it->first].push_back(dof1); + if (find(sendNewDofs[it->first].begin(), sendNewDofs[it->first].end(), dof2) == + sendNewDofs[it->first].end()) + sendNewDofs[it->first].push_back(dof2); + DofContainer boundDOFs; + addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, - boundDOFs); - addAllEdgeDOFs(boundIt->rankObject.el, - boundIt->rankObject.ithObjAtBoundary, - boundDOFs); + boundDOFs); + addAllEdgeDOFs(boundIt->rankObject.el, + boundIt->rankObject.ithObjAtBoundary, + boundDOFs); for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) { newBoundaryDOFs[boundDOFs[i]] = mpiRank; sendNewDofs[it->first].push_back(boundDOFs[i]); } + } } - for (std::map<int, std::vector<AtomicBoundary> >::iterator it = - otherIntBoundary.boundary.begin(); - it != otherIntBoundary.boundary.end(); - ++it) { + for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); + it != otherIntBoundary.boundary.end(); ++it) { + for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin(); - boundIt != it->second.end(); - ++boundIt) { + boundIt != it->second.end(); ++boundIt) { + const DegreeOfFreedom *dof1 = NULL; const DegreeOfFreedom *dof2 = NULL; @@ -812,8 +800,8 @@ namespace AMDiS { dof2 = boundIt->rankObject.el->getDOF(2); break; case 2: - dof1 = boundIt->rankObject.el->getDOF(0); - dof2 = boundIt->rankObject.el->getDOF(1); + dof1 = boundIt->rankObject.el->getDOF(1); + dof2 = boundIt->rankObject.el->getDOF(0); break; default: ERROR_EXIT("Should never happen!\n"); @@ -829,23 +817,23 @@ namespace AMDiS { newBoundaryDOFs[dof1] = boundaryDOFs[dof1]; newBoundaryDOFs[dof2] = boundaryDOFs[dof2]; - if (find(recvNewDofs[it->first].begin(), recvNewDofs[it->first].end(), dof2) == - recvNewDofs[it->first].end()) - recvNewDofs[it->first].push_back(dof2); - if (find(recvNewDofs[it->first].begin(), recvNewDofs[it->first].end(), dof1) == - recvNewDofs[it->first].end()) - recvNewDofs[it->first].push_back(dof1); - + if (find(recvNewDofs[it->first].begin(), recvNewDofs[it->first].end(), dof1) == + recvNewDofs[it->first].end()) + recvNewDofs[it->first].push_back(dof1); + if (find(recvNewDofs[it->first].begin(), recvNewDofs[it->first].end(), dof2) == + recvNewDofs[it->first].end()) + recvNewDofs[it->first].push_back(dof2); + DofContainer boundDOFs; + addAllEdgeDOFs(boundIt->rankObject.el, - boundIt->rankObject.ithObjAtBoundary, - boundDOFs); + boundIt->rankObject.ithObjAtBoundary, + boundDOFs); addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, boundDOFs); for (int i = static_cast<int>(boundDOFs.size()) - 1; i >= 0; i--) { - // for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) { rankDOFs.erase(boundDOFs[i]); newBoundaryDOFs[boundDOFs[i]] = it->first; recvNewDofs[it->first].push_back(boundDOFs[i]); @@ -875,15 +863,13 @@ namespace AMDiS { int i = 0; for (std::set<const DegreeOfFreedom*>::iterator dofIt = rankDOFs.begin(); - dofIt != rankDOFs.end(); - ++dofIt, i++) { - const_cast<DegreeOfFreedom*>(*dofIt)[0] = i; + dofIt != rankDOFs.end(); ++dofIt, i++) { + *(const_cast<DegreeOfFreedom*>(*dofIt)) = i; mapLocalGlobalDOFs[i] = rstart + i; mapGlobalLocalDOFs[rstart + i] = i; isRankDOF[i] = true; } - // === Send new DOF indices. === std::vector<int*> sendBuffers(sendNewDofs.size()); @@ -894,14 +880,12 @@ namespace AMDiS { i = 0; for (RankToDofContainer::iterator sendIt = sendNewDofs.begin(); - sendIt != sendNewDofs.end(); - ++sendIt, i++) { + sendIt != sendNewDofs.end(); ++sendIt, i++) { int nSendDofs = sendIt->second.size(); sendBuffers[i] = new int[nSendDofs]; int c = 0; for (DofContainer::iterator dofIt = sendIt->second.begin(); - dofIt != sendIt->second.end(); - ++dofIt) + dofIt != sendIt->second.end(); ++dofIt) sendBuffers[i][c++] = (*dofIt)[0] + rstart; request[requestCounter++] = @@ -910,8 +894,7 @@ namespace AMDiS { i = 0; for (RankToDofContainer::iterator recvIt = recvNewDofs.begin(); - recvIt != recvNewDofs.end(); - ++recvIt, i++) { + recvIt != recvNewDofs.end(); ++recvIt, i++) { int nRecvDofs = recvIt->second.size(); recvBuffers[i] = new int[nRecvDofs]; @@ -921,22 +904,18 @@ namespace AMDiS { MPI::Request::Waitall(requestCounter, request); - i = 0; for (RankToDofContainer::iterator sendIt = sendNewDofs.begin(); - sendIt != sendNewDofs.end(); - ++sendIt) + sendIt != sendNewDofs.end(); ++sendIt) delete [] sendBuffers[i++]; i = 0; int newDofIndex = nRankDOFs; for (RankToDofContainer::iterator recvIt = recvNewDofs.begin(); - recvIt != recvNewDofs.end(); - ++recvIt) { + recvIt != recvNewDofs.end(); ++recvIt) { int j = 0; for (DofContainer::iterator dofIt = recvIt->second.begin(); - dofIt != recvIt->second.end(); - ++dofIt) { + dofIt != recvIt->second.end(); ++dofIt) { (*const_cast<DegreeOfFreedom*>(*dofIt)) = newDofIndex; mapLocalGlobalDOFs[newDofIndex] = recvBuffers[i][j]; mapGlobalLocalDOFs[recvBuffers[i][j]] = newDofIndex; @@ -952,6 +931,7 @@ namespace AMDiS { sendDofs = sendNewDofs; recvDofs = recvNewDofs; + } @@ -997,23 +977,25 @@ namespace AMDiS { switch (ithEdge) { case 0: - if (el->getSecondChild()) + if (el->getSecondChild()) { addAllEdgeDOFs(el->getSecondChild(), 2, dofs); - else + } else { addThisEdge = true; + } break; case 1: - if (el->getFirstChild()) + if (el->getFirstChild()) { addAllEdgeDOFs(el->getFirstChild(), 2, dofs); - else + } else { addThisEdge = true; + } break; case 2: if (el->getFirstChild()) { - addAllEdgeDOFs(el->getFirstChild(), 1, dofs); addAllEdgeDOFs(el->getFirstChild(), 0, dofs); + addAllEdgeDOFs(el->getSecondChild(), 1, dofs); } else { addThisEdge = true; } @@ -1024,7 +1006,7 @@ namespace AMDiS { } if (addThisEdge) { - ElementDofIterator elDofIter(feSpace->getAdmin()); + ElementDofIterator elDofIter(feSpace, true); elDofIter.reset(el); do { if (elDofIter.getCurrentPos() == 1 && @@ -1042,13 +1024,12 @@ namespace AMDiS { { // === Determine to each dof the set of partitions the dof belongs to. === - ElementDofIterator elDofIt(feSpace->getAdmin()); + ElementDofIterator elDofIt(feSpace); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *element = elInfo->getElement(); - elDofIt.reset(element); do { // Determine to each dof the partition(s) it corresponds to. @@ -1104,24 +1085,80 @@ namespace AMDiS { } - void ParallelDomainBase::testInteriorBoundary() + void ParallelDomainBase::DbgTestInteriorBoundary() + { + FUNCNAME("ParallelDomainBase::DbgTestInteriorBoundary()"); + + std::vector<int*> sendBuffers, recvBuffers; + + MPI::Request request[myIntBoundary.boundary.size() + + otherIntBoundary.boundary.size()]; + int requestCounter = 0; + + for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin(); + rankIt != myIntBoundary.boundary.end(); ++rankIt) { + int nSendInt = rankIt->second.size(); + int* buffer = new int[nSendInt]; + for (int i = 0; i < nSendInt; i++) + buffer[i] = (rankIt->second)[i].rankObject.el->getIndex(); + sendBuffers.push_back(buffer); + + request[requestCounter++] = + mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0); + } + + for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); + rankIt != otherIntBoundary.boundary.end(); ++rankIt) { + int nRecvInt = rankIt->second.size(); + int *buffer = new int[nRecvInt]; + recvBuffers.push_back(buffer); + + request[requestCounter++] = + mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0); + } + + + MPI::Request::Waitall(requestCounter, request); + + for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++) + delete [] sendBuffers[i]; + + int bufCounter = 0; + for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); + rankIt != otherIntBoundary.boundary.end(); ++rankIt) { + + TEST_EXIT(rankIt->second.size() == otherIntBoundary.boundary[rankIt->first].size()) + ("Boundaries does not fit together!\n"); + + for (int i = 0; i < static_cast<int>(rankIt->second.size()); i++) { + int elIndex1 = recvBuffers[bufCounter][i]; + int elIndex2 = otherIntBoundary.boundary[rankIt->first][i].neighbourObject.el->getIndex(); + + TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at interior boundary!\n"); + } + + delete [] recvBuffers[bufCounter++]; + } + } + + + void ParallelDomainBase::DbgTestCommonDofs() { - FUNCNAME("ParallelDomainBase::testInteriorBoundary()"); + FUNCNAME("ParallelDomainBase::DbgTestCommonDofs()"); // Maps to each neighbour rank an array of WorldVectors. This array contain the // coordinates of all dofs this rank shares on the interior boundary with the // neighbour rank. A rank sends the coordinates to another rank, if it owns the // boundarys dof. - std::map<int, std::vector<WorldVector<double> > > sendCoords; + RankToCoords sendCoords; // A rank receives all boundary dofs that are at its interior boundaries but are // not owned by the rank. This map stores for each rank the the coordinates of dofs // this rank expectes to receive from. - std::map<int, std::vector<WorldVector<double> > > recvCoords; + RankToCoords recvCoords; WorldVector<double> coords; - for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator it = - sendDofs.begin(); + for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it) { for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) { @@ -1131,8 +1168,7 @@ namespace AMDiS { } } - for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator it = - recvDofs.begin(); + for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) { @@ -1148,16 +1184,11 @@ namespace AMDiS { MPI::Request request[(mpiSize - 1) * 2]; int requestCounter = 0; - for (std::map<int, std::vector<WorldVector<double> > >::iterator it - = sendCoords.begin(); - it != sendCoords.end(); ++it) { + for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) sendSize[it->first] = it->second.size(); - } - for (std::map<int, std::vector<WorldVector<double> > >::iterator it = recvCoords.begin(); - it != recvCoords.end(); ++it) { + for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) recvSize[it->first] = it->second.size(); - } for (int i = 0; i < mpiSize; i++) { if (i == mpiRank) @@ -1189,13 +1220,10 @@ namespace AMDiS { } } - std::map<int, double*> sendCoordsBuffer; - std::map<int, double*> recvCoordsBuffer; + std::map<int, double*> sendCoordsBuffer, recvCoordsBuffer; requestCounter = 0; - for (std::map<int, std::vector<WorldVector<double> > >::iterator it - = sendCoords.begin(); - it != sendCoords.end(); ++it) { + for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) { sendCoordsBuffer[it->first] = new double[it->second.size() * 2]; for (int i = 0; i < static_cast<int>(it->second.size()); i++) { @@ -1208,8 +1236,7 @@ namespace AMDiS { MPI_DOUBLE, it->first, 0); } - for (std::map<int, std::vector<WorldVector<double> > >::iterator it = recvCoords.begin(); - it != recvCoords.end(); ++it) { + for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) { recvCoordsBuffer[it->first] = new double[it->second.size() * 2]; request[requestCounter++] = mpiComm.Irecv(recvCoordsBuffer[it->first], @@ -1219,19 +1246,27 @@ namespace AMDiS { MPI::Request::Waitall(requestCounter, request); - for (std::map<int, std::vector<WorldVector<double> > >::iterator it - = sendCoords.begin(); - it != sendCoords.end(); ++it) { + for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) delete [] sendCoordsBuffer[it->first]; - } - for (std::map<int, std::vector<WorldVector<double> > >::iterator it = recvCoords.begin(); - it != recvCoords.end(); ++it) { + for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) { for (int i = 0; i < static_cast<int>(it->second.size()); i++) { - if ((it->second)[i][0] != recvCoordsBuffer[it->first][i * 2] || - (it->second)[i][1] != recvCoordsBuffer[it->first][i * 2 + 1]) - ERROR_EXIT("WRONG COORDS!\n"); + std::cout << "[DBG] i = " << i << std::endl; + + std::cout.precision(4); + std::cout << "[DBG] " + << "Rank " << mpiRank << " from rank " << it->first + << " expect coords (" + << (it->second)[i][0] << " , " << (it->second)[i][1] + << ") received coords (" + << recvCoordsBuffer[it->first][i * 2] << " , " + << recvCoordsBuffer[it->first][i * 2 + 1] << ")" << std::endl; + + bool c0 = (it->second)[i][0] == recvCoordsBuffer[it->first][i * 2]; + bool c1 = (it->second)[i][1] == recvCoordsBuffer[it->first][i * 2 + 1]; + + TEST_EXIT(c0 && c1)("Wrong DOFs!\n"); } delete [] recvCoordsBuffer[it->first]; diff --git a/AMDiS/src/ParallelDomainProblem.h b/AMDiS/src/ParallelDomainProblem.h index 462c189ac813315b0fd1feb3ce0fe44824483ab2..94f6cd7c63474e692e13c2a1f90bce84b1be1f1e 100644 --- a/AMDiS/src/ParallelDomainProblem.h +++ b/AMDiS/src/ParallelDomainProblem.h @@ -64,6 +64,12 @@ namespace AMDiS { /// Defines a mapping type from DOF indices to boolean values. typedef std::map<DegreeOfFreedom, bool> DofToBool; + /// Defines a mapping type from rank numbers to sets of coordinates. + typedef std::map<int, std::vector<WorldVector<double> > > RankToCoords; + + /// Forward type (it maps rank numbers to the interior boundary objects). + typedef InteriorBoundary::RankToBoundMap RankToBoundMap; + public: ParallelDomainBase(const std::string& name, ProblemIterationInterface *iterationIF, @@ -196,6 +202,8 @@ namespace AMDiS { DofContainer& rankDOFs, DofToRank& boundaryDOFs); + + void DbgTestInteriorBoundary(); /** \brief * This function is used for debugging only. It traverses all interior boundaries @@ -203,7 +211,7 @@ namespace AMDiS { * neighbours. The function fails, when dof indices on an interior boundary does * not fit together. */ - void testInteriorBoundary(); + void DbgTestCommonDofs(); protected: ///