AND, YEAH, IT WORKS....

AMDiS/src/Element.cc

@@ -651,8 +651,14 @@ namespace AMDiS {
 			   bool baseDofPtr,
 			   vector<GeoIndex>* dofGeoIndex)
   {
+    FUNCNAME("Element::getAllDofs()");
+
+    int sb = dofs.size();
+
     getNodeDofs(feSpace, bound, dofs, baseDofPtr);
 
+    //    MSG("    -> getNodeDofs %d\n", dofs.size() - sb);
+
     if (dofGeoIndex != NULL) {
       // In the case dofGeoIndex should be filled, set all node DOFs to be 
       // vertex DOFs.
@@ -661,8 +667,14 @@ namespace AMDiS {
 	(*dofGeoIndex)[i] = VERTEX;
     }
 
-    if (feSpace->getBasisFcts()->getDegree() > 1)
+    if (feSpace->getBasisFcts()->getDegree() > 1) {
+
+      sb = dofs.size();
+
       getHigherOrderDofs(feSpace, bound, dofs, baseDofPtr, dofGeoIndex);
+      
+      //      MSG("     -> getHODofs %d\n", dofs.size() - sb);
+    }
 
     if (dofGeoIndex) {
       TEST_EXIT_DBG(dofs.size() == dofGeoIndex->size())

AMDiS/src/Tetrahedron.cc

@@ -249,8 +249,7 @@ namespace AMDiS {
     case 2:
     case 3:
       {
-	int n0 = 
-	  (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(FACE));
+	int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(VERTEX));
 	BoundaryObject nextBound0 = bound, nextBound1 = bound;
 	prepareNextBound(nextBound0, 0);
 	prepareNextBound(nextBound1, 1);
@@ -296,7 +295,7 @@ namespace AMDiS {
     if (!child[0])
       return;    
 
-    int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(EDGE));
+    int n0 = (baseDofPtr ? 0 : feSpace->getAdmin()->getNumberOfPreDofs(VERTEX));
     BoundaryObject nextBound0 = bound, nextBound1 = bound;
     prepareNextBound(nextBound0, 0);
     prepareNextBound(nextBound1, 1);

AMDiS/src/parallel/DofComm.cc

@@ -14,6 +14,8 @@
 #include "parallel/InteriorBoundary.h"
 #include "parallel/MeshLevelData.h"
 #include "FiniteElemSpace.h"
+#include "Debug.h"
+#include "ElementDofIterator.h"
 
 namespace AMDiS {
 
@@ -94,33 +96,10 @@ namespace AMDiS {
   {
     FUNCNAME("DofComm:serialize()");
 
-    MSG("MUSS DAS WIRKLICH SEIN????\n");
+    ERROR_EXIT("MUSS DAS WIRKLICH SEIN????\n");
   }
 
   
-  void DofComm::deserialize(istream &in, 
-			    map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > dofIndexMap)
-  {
-    FUNCNAME("DofComm::deserialize()");
-
-    MSG("MUSS DAS WIRKLICH SEIN????\n");
-  }
-
-
-  void DofComm::serialize(LevelDataType &data, ostream &out)
-  {
-    FUNCNAME("DofComm:serialize()");
-  }
-
-   
-  void DofComm::deserialize(LevelDataType &data,
-			    istream &in, 
-			    map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > dofIndexMap)
-  {
-    FUNCNAME("DofComm::deserialize()");
-  }
-
-
   int DofComm::getNumberDofs(LevelDataType &data, 
 			     int level, 
 			     const FiniteElemSpace *feSpace,

AMDiS/src/parallel/DofComm.h

@@ -76,10 +76,6 @@ namespace AMDiS {
     // Writes all data of this object to an output stream.
     void serialize(ostream &out);
 
-    // Reads the object data from an input stream.
-    void deserialize(istream &in, 
-		     map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > dofIndexMap);
-
     int getNumberDofs(LevelDataType &data, 
 		      int level, 
 		      const FiniteElemSpace *feSpace,
@@ -88,12 +84,6 @@ namespace AMDiS {
   protected:
     void createContainer(RankToBoundMap &boundary, LevelDataType &data);
 
-    void serialize(LevelDataType &data, ostream &out);
-
-    void deserialize(LevelDataType &data,
-		     istream &in, 
-		     map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > dofIndexMap);
-
   protected:
     /// This map contains for each rank the list of DOFs the current rank must 
     /// end to exchange solution DOFs at the interior boundaries.

AMDiS/src/parallel/MeshDistributor.cc

@@ -614,7 +614,6 @@ namespace AMDiS {
 
       for (int i = 0; i < vec.getSize(); i++) {
 	DOFVector<double> &dofVec = *(vec.getDOFVector(i));
-
 	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
 	  dofs.push_back(dofVec[it.getDofIndex()]);
       }

AMDiS/src/parallel/PetscSolverNavierStokes.cc

@@ -57,7 +57,6 @@ namespace AMDiS {
 		    pressureComponent);
     TEST_EXIT(pressureComponent >= 0)
       ("For using PETSc stokes solver you must define a pressure component!\n");
-    TEST_EXIT(pressureComponent == 2)("Fixed for pressure component = 2\n");
 
     Parameters::get(initFileStr + "->navierstokes->pressure null space",
 		    pressureNullSpace);
@@ -89,9 +88,13 @@ namespace AMDiS {
     TEST_EXIT(invTau)("invtau pointer not set!\n");
     TEST_EXIT(solution)("solution pointer not set!\n");
 
+    int dim = componentSpaces[pressureComponent]->getMesh()->getDim();
+
     vector<int> velocityComponents;
     velocityComponents.push_back(0);
     velocityComponents.push_back(1);
+    if (dim == 3)
+      velocityComponents.push_back(2);
 
     PCSetType(pc, PCFIELDSPLIT);
     PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
@@ -139,8 +142,6 @@ namespace AMDiS {
     massMatrixSolver->fillPetscMatrix(&massMatrix);
 
 
-    //    VtkWriter::writeFile(phase, "phase.vtu");
-
     // === Laplace matrix solver ===
 
     DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
@@ -158,9 +159,12 @@ namespace AMDiS {
 
     DOFVector<double> vx(pressureFeSpace, "vx");
     DOFVector<double> vy(pressureFeSpace, "vy");
+    DOFVector<double> vz(pressureFeSpace, "vz");
     DOFVector<double> vp(pressureFeSpace, "vp");
     vx.interpol(solution->getDOFVector(0));
     vy.interpol(solution->getDOFVector(1));
+    if (dim == 3)
+      vz.interpol(solution->getDOFVector(2));
 
     DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
     Operator conDifOp(pressureFeSpace, pressureFeSpace);
@@ -175,6 +179,11 @@ namespace AMDiS {
       conDifOp.addTerm(conDif2, GRD_PHI);
       VecAtQP_FOT *conDif3 = new VecAtQP_FOT(&vy, &idFct, 1);
       conDifOp.addTerm(conDif3, GRD_PHI);
+
+      if (dim == 3) {
+	VecAtQP_FOT *conDif4 = new VecAtQP_FOT(&vz, &idFct, 2);	
+	conDifOp.addTerm(conDif4, GRD_PHI);
+      }
     } else {
       MSG("INIT WITH PHASE!\n");
       
@@ -190,6 +199,11 @@ namespace AMDiS {
 
       Vec2AtQP_FOT *conDif3 = new Vec2AtQP_FOT(&vy, &vp, new Multiplier3(), 1);
       conDifOp.addTerm(conDif3, GRD_PHI);
+
+      if (dim == 3) {
+	Vec2AtQP_FOT *conDif4 = new Vec2AtQP_FOT(&vz, &vp, new Multiplier3(), 2);
+	conDifOp.addTerm(conDif4, GRD_PHI);
+      }
     }
 
     conDifMatrix.assembleOperator(conDifOp);