PetscSolverNSCH is now fit for 3D

AMDiS/src/parallel/PetscSolverNSCH.cc

@@ -31,18 +31,16 @@ namespace AMDiS {
   CahnHilliardData2* data = static_cast<CahnHilliardData2*>(ctx);
   
   /// extract vectors
-  Vec b1, b2, b34, b5, x1, x2, x34, x5, b3,b4;
-  data->globalMatrixSolver->extractVectorComponent(b, 0+3, &b1);
-  data->globalMatrixSolver->extractVectorComponent(b, 1+3, &b2);
-  data->globalMatrixSolver->extractVectorComponent(b, 0, &b34, 2);  
-  data->globalMatrixSolver->extractVectorComponent(b, 2, &b5);
+  Vec b1, b2, b34, b5, x1, x2, x34, x5;
+  data->globalMatrixSolver->extractVectorComponent(b, data->dim+1, &b1);
+  data->globalMatrixSolver->extractVectorComponent(b, data->dim+2, &b2);
+  data->globalMatrixSolver->extractVectorComponent(b, 0, &b34, data->dim);  
+  data->globalMatrixSolver->extractVectorComponent(b, data->dim, &b5);
   
-  data->globalMatrixSolver->extractVectorComponent(x, 0+3, &x1);
-  data->globalMatrixSolver->extractVectorComponent(x, 1+3, &x2);
-  data->globalMatrixSolver->extractVectorComponent(x, 0, &x34, 2);  
-  data->globalMatrixSolver->extractVectorComponent(x, 2, &x5);
-  data->globalMatrixSolver->extractVectorComponent(b, 0, &b3, 1);  
-  data->globalMatrixSolver->extractVectorComponent(b, 1, &b4, 1);  
+  data->globalMatrixSolver->extractVectorComponent(x, data->dim+1, &x1);
+  data->globalMatrixSolver->extractVectorComponent(x, data->dim+2, &x2);
+  data->globalMatrixSolver->extractVectorComponent(x, 0, &x34, data->dim);  
+  data->globalMatrixSolver->extractVectorComponent(x, data->dim, &x5);
     
   
   /// solve Cahn-Hilliard Preconditioner  
@@ -190,7 +188,7 @@ namespace AMDiS {
     double wtime = MPI::Wtime();
     int dim = componentSpaces[0]->getMesh()->getDim();     
     pressureComponent=dim;
-    const FiniteElemSpace *cahnHilliardFeSpace = componentSpaces[0+3];
+    const FiniteElemSpace *cahnHilliardFeSpace = componentSpaces[dim+1];
     const FiniteElemSpace *velocityFeSpace= componentSpaces[0];
     const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
         
@@ -199,6 +197,7 @@ namespace AMDiS {
     PCShellSetContext(pc, &matShellContext);
     matShellContext.globalMatrixSolver = (this);
     matShellContext.mpiCommGlobal= &(meshDistributor->getMpiComm(0));
+    matShellContext.dim = dim;
     
     
     /// Init Cahn-Hilliard part
@@ -222,16 +221,16 @@ namespace AMDiS {
       laplaceOp2.addTerm(&sot); // -delta*K
       
       massMatrixCH.assembleOperator(massOpCH);        
-      massMatrixSolverCH = createSubSolver(0+3, "mass_");       
+      massMatrixSolverCH = createSubSolver(dim+1, "mass_");       
       massMatrixSolverCH->fillPetscMatrix(&massMatrixCH);     
       // === matrix (M + eps*sqrt(delta)*K) ===
       laplaceMatrixCH.assembleOperator(laplaceOpCH);
-      laplaceMatrixSolverCH = createSubSolver(0+3, "laplace_");
+      laplaceMatrixSolverCH = createSubSolver(dim+1, "laplace_");
       laplaceMatrixSolverCH->fillPetscMatrix(&laplaceMatrixCH);     
       
       // === matrix (-delta*K) ===
       deltaKMatrix.assembleOperator(laplaceOp2);
-      deltaKMatrixSolver = createSubSolver(0+3, "laplace2_");
+      deltaKMatrixSolver = createSubSolver(dim+1, "laplace2_");
       deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);      
     }
     
@@ -317,9 +316,9 @@ namespace AMDiS {
     conDifMatrixSolver->fillPetscMatrix(&conDifMatrix);
     matShellContext.matConDif = conDifMatrixSolver->getMatInterior();    
     
-    extractMatrixComponent(mat[0][0], 2, 1, 0, 2, &(matShellContext.matDiv));
-    extractMatrixComponent(mat[0][0], 0, 2, 2, 1, &(matShellContext.matGrad));
-    extractMatrixComponent(mat[0][0], 0, 2, 0, 2, &(matShellContext.velocityMat));
+    extractMatrixComponent(mat[0][0], dim, 1, 0, dim, &(matShellContext.matDiv));
+    extractMatrixComponent(mat[0][0], 0, dim, dim, 1, &(matShellContext.matGrad));
+    extractMatrixComponent(mat[0][0], 0, dim, 0, dim, &(matShellContext.velocityMat));
     
     ///erstelle kspVelocity
     KSPCreate((meshDistributor->getMpiComm(0)), &(matShellContext.kspVelocity));

AMDiS/src/parallel/PetscSolverNSCH.h

@@ -31,6 +31,7 @@ namespace AMDiS {
     Mat matMassCH, matMinusDeltaK, matGrad, matDiv, matConDif, matSchur, velocityMat;
     PetscSolverGlobalMatrix *globalMatrixSolver;
     double *eps, *delta;
+    int dim;
     MPI::Intracomm *mpiCommGlobal;
   };