* Output support for higher order Lagrange Elements

AMDiS/src/DataCollector.cc

@@ -116,6 +116,7 @@ namespace AMDiS {
       elInfo = stack.traverseNext(elInfo);
     }
 
+
     // Remove all interpolation marks and, instead, set to each
     // interpolation point its continous index starting from 0.
     int i = 0;
@@ -126,7 +127,8 @@ namespace AMDiS {
     }
 
     // Traverse elements to create interpolation values.
-    elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_);
+    counter = 0;
+    elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_ | Mesh::FILL_COORDS);
     while (elInfo) {
       if (!writeElem_ || writeElem_(elInfo))
 	addInterpData(elInfo);
@@ -268,102 +270,73 @@ namespace AMDiS {
   {
     FUNCNAME("DataCollector::addValueData()");
 
-    const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
-
-    /* vertex dofs */
-    int dof_offset = localAdmin_->getNumberOfPreDOFs(VERTEX);
-
+    const BasisFunction *basisFcts = feSpace_->getBasisFcts();
+    const DegreeOfFreedom *localDOFs = basisFcts->getLocalIndices(elInfo->getElement(), localAdmin_, NULL);
+    const int nBasisFcts = basisFcts->getNumber();
+  
+       
     for (int i = 0; i < mesh_->getGeo(VERTEX); i++) {
-      (*interpPointInd_)[dof[i][dof_offset]] = -2; // mark as vertex
+      (*interpPointInd_)[localDOFs[i]] = -2; // mark as vertex
       
       // get coords of this vertex
       WorldVector<double> vertexCoords = elInfo->getCoord(i);
       
       // search for coords at this dof
       ::std::list<WorldVector<double> >::iterator it =
-	  find((*dofCoords_)[dof[i][dof_offset]].begin(),
-	       (*dofCoords_)[dof[i][dof_offset]].end(),
+	  find((*dofCoords_)[localDOFs[i]].begin(),
+	       (*dofCoords_)[localDOFs[i]].end(),
 	       vertexCoords);
       
       // coords not yet in list?
-      if (it == (*dofCoords_)[dof[i][dof_offset]].end()) {
+      if (it == (*dofCoords_)[localDOFs[i]].end()) {
 	// add new coords to list
-	(*dofCoords_)[dof[i][dof_offset]].push_back(vertexCoords);
+	(*dofCoords_)[localDOFs[i]].push_back(vertexCoords);
       }
     }
+   
 
-    int nInterpPoints = 0;
-    const BasisFunction *basisFcts = feSpace_->getBasisFcts();
-
-    for (int i = 1; i <= dim_; i++) {
-      int num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_));
-      int node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
-      dof_offset = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_));
-
-      for (int j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) {
-	int node = node_offset + j;
-
-	for (int k = 0; k < num_dofs; k++) {
-	  int dof_index = dof_offset + k;
-
-
-	  WorldVector<double> interpolCoords;
-	  elInfo->coordToWorld((*basisFcts->getCoords(mesh_->getGeo(VERTEX) + nInterpPoints)), 
-			       &interpolCoords);
-
-	  nInterpPoints++;
-	  
-	  if ((*interpPointInd_)[dof[node][dof_index]] == -1) {
-	    // mark as interpolation point
-	    (*interpPointInd_)[dof[node][dof_index]] = -3; 
-
-	    // search for interpolation point coordinates, and insert them to the 
-	    // dof-entry, if not contained in the list
-	    ::std::list<WorldVector<double> >::iterator it =
-		find((*interpPointCoords_)[dof[node][dof_index]].begin(),
-		     (*interpPointCoords_)[dof[node][dof_index]].end(),
+    for (int i = mesh_->getGeo(VERTEX); i < nBasisFcts; i++) {
+      WorldVector<double> interpolCoords;
+      elInfo->coordToWorld(*basisFcts->getCoords(i), &interpolCoords);
+      
+      if ((*interpPointInd_)[localDOFs[i]] == -1) {
+	// mark as interpolation point
+	(*interpPointInd_)[localDOFs[i]] = -3; 
+	
+	// search for interpolation point coordinates, and insert them to the 
+	// dof-entry, if not contained in the list
+	::std::list<WorldVector<double> >::iterator it =
+		find((*interpPointCoords_)[localDOFs[i]].begin(),
+		     (*interpPointCoords_)[localDOFs[i]].end(),
 		     interpolCoords);
-	    if (it == (*interpPointCoords_)[dof[node][dof_index]].end()) {
-	      (*interpPointCoords_)[dof[node][dof_index]].push_back(interpolCoords); 
-	    }
-
-
-	    nInterpPoints_++;
-	  }
+	
+	if (it == (*interpPointCoords_)[localDOFs[i]].end()) {
+	  (*interpPointCoords_)[localDOFs[i]].push_back(interpolCoords); 
+	  nInterpPoints_++;
 	}
-      }
+      }      
     }
-  
+ 
     return(0);
   }
 
   int DataCollector::addInterpData(ElInfo *elInfo)
   {
     FUNCNAME("DataCollector::addInterpData()");
-
-    const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
-
+    
     ::std::vector<int> elemInterpPoints;
     elemInterpPoints.clear();
 
-    for (int i = 1; i <= dim_; i++) {
-      int num_dofs = localAdmin_->getNumberOfDOFs(INDEX_OF_DIM(i, dim_));
-      int node_offset = mesh_->getNode(INDEX_OF_DIM(i, dim_));
-      int dof_offset = localAdmin_->getNumberOfPreDOFs(INDEX_OF_DIM(i, dim_));
-
-      for (int j = 0; j < mesh_->getGeo(INDEX_OF_DIM(i, dim_)); j++) {
-	int node = node_offset + j;
+    const BasisFunction *basisFcts = feSpace_->getBasisFcts();
+    const DegreeOfFreedom *localDOFs = basisFcts->getLocalIndices(elInfo->getElement(), localAdmin_, NULL);
+    const int nBasisFcts = basisFcts->getNumber();
 
-	for (int k = 0; k < num_dofs; k++) {
-	  int dof_index = dof_offset + k;
-	  
-	  elemInterpPoints.push_back((*interpPointInd_)[dof[node][dof_index]]);
-	}
-      }
+    for (int i = mesh_->getGeo(VERTEX); i < nBasisFcts; i++) {
+      elemInterpPoints.push_back((*interpPointInd_)[localDOFs[i]]);
     }
 
     interpPoints_.push_back(elemInterpPoints);
-
+  
     return(0);
   }
 

AMDiS/src/DataCollector.h

@@ -294,6 +294,8 @@ namespace AMDiS {
        * Pointer to a function which decides whether an element is considered.
        */
       bool (*writeElem_)(ElInfo*);
+
+      int counter;
     };
 }
 

AMDiS/src/Element.cc

@@ -12,8 +12,11 @@ namespace AMDiS {
   {
     ElementRegion_ED* red_;
 
-    if (!elementData) return -1;
-    red_=dynamic_cast<ElementRegion_ED*>(elementData->getElementData(ELEMENT_REGION));
+    if (!elementData) 
+      return -1;
+
+    red_ = dynamic_cast<ElementRegion_ED*>(elementData->getElementData(ELEMENT_REGION));
+
     if (red_) 
       return red_->getRegion();
     
@@ -22,8 +25,9 @@ namespace AMDiS {
 
   void Element::setDOFPtrs()
   {
-    FUNCNAME("Element::setDOFPtrs");
+    FUNCNAME("Element::setDOFPtrs()");
     TEST_EXIT(mesh)("no mesh!\n");
+
     dof = mesh->createDOFPtrs();  
   }
 
@@ -33,27 +37,30 @@ namespace AMDiS {
 
     index = mesh ? mesh->getNextElementIndex() : -1; 
   
-    child[0]  = NULL;
-    child[1]  = NULL;
+    child[0] = NULL;
+    child[1] = NULL;
 
-    dof       = mesh ? mesh->createDOFPtrs() : NULL;
     newCoord = NULL;
-
     elementData = NULL;
+
+    if (mesh) {
+      setDOFPtrs();
+    } else {
+      mesh = NULL;
+    }
   }
 
   // call destructor through Mesh::freeElement !!!
   Element::~Element()
   {
-    if(child[0]) DELETE child[0];
-    if(child[1]) DELETE child[1];
-
-    //if(elementData) DELETE elementData;
+    if (child[0]) 
+      DELETE child[0];
+    if (child[1]) 
+      DELETE child[1];
 
-    if (newCoord)
-      {
-	DELETE newCoord;
-      }
+    if (newCoord) {
+      DELETE newCoord;
+    }
   }
 
   /****************************************************************************/
@@ -65,7 +72,7 @@ namespace AMDiS {
   /* CHANGE_DOFS_1 changes old dofs to NEGATIVE new dofs */
 
 #define CHANGE_DOFS_1(el)						\
-  ldof = el->dof[n0+i] + nd0;						\
+  ldof = el->dof[n0 + i] + nd0;						\
     for (j = 0; j < nd; j++) {						\
       if ((k = ldof[j]) >= 0) {						\
 	/* do it only once! (dofs are visited more than once) */	\
@@ -99,7 +106,7 @@ namespace AMDiS {
       }
     }
 
-    if(mesh->getDim() > 1) {
+    if (mesh->getDim() > 1) {
       if ((nd = admin->getNumberOfDOFs(EDGE)))  {
 	nd0 = admin->getNumberOfPreDOFs(EDGE);
 	n0 = admin->getMesh()->getNode(EDGE);
@@ -109,7 +116,7 @@ namespace AMDiS {
       }
     }
 
-    if (3==mesh->getDim()) {
+    if (mesh->getDim() == 3) {
       if ((nd = admin->getNumberOfDOFs(FACE)))  {
 	nd0 = admin->getNumberOfPreDOFs(FACE);
 	n0 = admin->getMesh()->getNode(FACE);
@@ -145,7 +152,7 @@ namespace AMDiS {
       }
     }
 
-    if(mesh->getDim() > 1) {
+    if (mesh->getDim() > 1) {
       if ((nd = admin->getNumberOfDOFs(EDGE)))  {
 	nd0 = admin->getNumberOfPreDOFs(EDGE);
 	n0 = admin->getMesh()->getNode(EDGE);
@@ -155,7 +162,7 @@ namespace AMDiS {
       }
     }
 
-    if (3==mesh->getDim()) {
+    if (mesh->getDim() == 3) {
       if ((nd = admin->getNumberOfDOFs(FACE)))  {
 	nd0 = admin->getNumberOfPreDOFs(FACE);
 	n0 = admin->getMesh()->getNode(FACE);
@@ -183,30 +190,30 @@ namespace AMDiS {
 
   int Element::oppVertex(FixVec<DegreeOfFreedom*, DIMEN> pdof) const
   {
-    int      i, j, nv = 0, ov = 0;
+    int nv = 0, ov = 0;
 
     int vertices = mesh->getGeo(VERTEX);
     int dim = mesh->getDim();
 
-    for (i = 0; i < vertices; i++)
-      {
-	if (nv < i-1)  return(-1);
-
-	for (j = 0; j < dim; j++)
-	  {
-	    if (dof[i] == pdof[j])
-	      {
-		/****************************************************************************/
-		/* i is a common vertex                                                     */
-		/****************************************************************************/
-		ov += i;
-		nv++;
-		break;
-	      }
-	  }
-
+    for (int i = 0; i < vertices; i++) {
+      if (nv < i-1)  
+	return(-1);
+
+      for (int j = 0; j < dim; j++) {
+	if (dof[i] == pdof[j]) {
+	  /****************************************************************************/
+	  /* i is a common vertex                                                     */
+	  /****************************************************************************/
+	  ov += i;
+	  nv++;
+	  break;
+	}
       }
-    if (nv != mesh->getDim()) return(-1);
+
+    }
+    
+    if (nv != mesh->getDim()) 
+      return(-1);
     /****************************************************************************/
     /*  the opposite vertex is 3(6) - (sum of indices of common vertices) in    */
     /*  2d(3d)                                                                  */
@@ -243,7 +250,7 @@ namespace AMDiS {
       newCoord=NULL;
     };
   }
-
+ 
   void Element::serialize(::std::ostream &out) 
   {
     // write children

AMDiS/src/Element.h

@@ -305,9 +305,9 @@ namespace AMDiS {
     /** \brief
      * Sets the pointer to the DOFs of the i-th node of Element
      */
-    DegreeOfFreedom* setDOF(int i, DegreeOfFreedom* p) {
-      dof[i] = p;
-      return dof[i];
+    DegreeOfFreedom* setDOF(int pos, DegreeOfFreedom* p) {
+      dof[pos] = p;
+      return dof[pos];
     };
 
     /** \brief

AMDiS/src/Lagrange.cc

@@ -360,7 +360,9 @@ namespace AMDiS {
     }
   }
 
-  Lagrange::Phi::~Phi() { DELETE [] vertices; }
+  Lagrange::Phi::~Phi() { 
+    DELETE [] vertices; 
+  }
 
 
   Lagrange::GrdPhi::GrdPhi(Lagrange* owner_, 
@@ -1003,16 +1005,17 @@ namespace AMDiS {
 
     DegreeOfFreedom* result;
 
-    if(indices) {
+    if (indices) {
       result = indices;
     } else {
-      if(localVec) FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize);
+      if (localVec) 
+	FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize);
       localVec = GET_MEMORY(DegreeOfFreedom, nBasFcts);
       localVecSize = nBasFcts;
       result = localVec;
     }
 
-    for(pos=0, j=0; pos <= dim; pos++) {
+    for (pos = 0, j = 0; pos <= dim; pos++) {
       posIndex = INDEX_OF_DIM(pos, dim);
       n0 = admin->getNumberOfPreDOFs(posIndex);
       node0 = admin->getMesh()->getNode(posIndex);
@@ -1032,33 +1035,6 @@ namespace AMDiS {
     return result;
   }
 
-  // const double* Lagrange::getVec(const Element* el, const DOFVector<double> * dv, 
-  // 			       double * d) const
-  // {
-  //   static double* localVec = NULL;
-  //   const DOFAdmin* admin = dv->getFESpace()->getAdmin();
-
-  //   int i;
-
-  //   double* result;
-  
-  //   if(d) {
-  //     result = d;
-  //   } else {
-  //     if(localVec) FREE_MEMORY(localVec, double, nBasFcts);
-  //     localVec = GET_MEMORY(double, nBasFcts);
-  //     result = localVec;
-  //   }
-
-  //   const DegreeOfFreedom *localIndices = getLocalIndices(el, admin, NULL);
-
-  //   for(i = 0; i < nBasFcts; i++) {
-  //     result[i] = (*dv)[localIndices[i]];
-  //   }
-
-  //   return result;
-  // }
-
   void Lagrange::l2ScpFctBas(Quadrature *q,
 			     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
 			     DOFVector<WorldVector<double> >* fh)
@@ -1234,12 +1210,13 @@ namespace AMDiS {
 				 RCNeighbourList* list, 
 				 int n, BasisFunction* basFct)
   {
-    FUNCNAME("Lagrange::refineInter2_2d");
+    FUNCNAME("Lagrange::refineInter2_2d()");
 
-    if (n < 1) return;
+    if (n < 1) 
+      return;
 
-    Element  *el;
-    int       node, n0;
+    Element *el;
+    int node, n0;
     DegreeOfFreedom cdof;
     const DegreeOfFreedom *pdof;
     const DOFAdmin *admin = drv->getFESpace()->getAdmin();
@@ -1267,7 +1244,7 @@ namespace AMDiS {
   
     cdof = el->getChild(0)->getDOF(node, n0); 
     (*drv)[cdof] = 
-      0.375*(*drv)[pdof[0]] - 0.125*(*drv)[pdof[1]]  + 0.75*(*drv)[pdof[5]];
+      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];
 
     /****************************************************************************/
     /* node in the common edge of child[0] and child[1]                         */
@@ -1275,8 +1252,8 @@ namespace AMDiS {
 
     cdof = el->getChild(0)->getDOF(node+1, n0); 
     (*drv)[cdof] =
-      -0.125*((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25*(*drv)[pdof[5]]
-      + 0.5*((*drv)[pdof[3]] + (*drv)[pdof[4]]);
+      -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
+      + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
 
     /****************************************************************************/
     /*  midpoint of edge on child[1] at the refinement edge                     */
@@ -1284,22 +1261,20 @@ namespace AMDiS {
   
     cdof = el->getChild(1)->getDOF(node+1, n0); 
     (*drv)[cdof] = 
-      -0.125*(*drv)[pdof[0]] + 0.375*(*drv)[pdof[1]]  + 0.75*(*drv)[pdof[5]];
+      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]]  + 0.75 * (*drv)[pdof[5]];
 
-    if (n > 1)
-      {
-	/****************************************************************************/
-	/*  adjust the value at the midpoint of the common edge of neigh's children */
-	/****************************************************************************/
-	el = list->getElement(1);
-	pdof = basFct->getLocalIndices(el, admin, NULL);
-
-	cdof = el->getChild(0)->getDOF(node+1, n0); 
-	(*drv)[cdof] = 
-	  -0.125*((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25*(*drv)[pdof[5]]
-	  + 0.5*((*drv)[pdof[3]] + (*drv)[pdof[4]]);
-      }
-    return;
+    if (n > 1) {
+      /****************************************************************************/
+      /*  adjust the value at the midpoint of the common edge of neigh's children */
+      /****************************************************************************/
+      el = list->getElement(1);
+      pdof = basFct->getLocalIndices(el, admin, NULL);
+      
+      cdof = el->getChild(0)->getDOF(node+1, n0); 
+      (*drv)[cdof] = 
+	-0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
+	+ 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
+    }
   }
 
   void Lagrange::refineInter2_3d(DOFIndexed<double> *drv, 

AMDiS/src/Lagrange.h

@@ -152,20 +152,6 @@ namespace AMDiS {
 					   const DOFAdmin*,