Work on integration of the BDDCML library.

AMDiS/src/parallel/BddcMlSolver.cc

+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+extern "C" {
+#include <bddcml_interface_c.h>
+}
+
+#include "parallel/BddcMlSolver.h"
+
+namespace AMDiS {
+
+#ifdef HAVE_BDDC_ML
+
+  void BddcMlSolver::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
+  {
+    FUNCNAME("BddcMlSolver::fillPetscMatrix()");
+  }
+
+
+  void BddcMlSolver::fillPetscRhs(SystemVector *vec)
+  {
+    FUNCNAME("BddcMlSolver::fillPetscRhs()");
+  }
+
+
+  void BddcMlSolver::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo)
+  {
+    FUNCNAME("BddcMlSolver::solvePetscMatrix()");
+    
+    int nComponents = vec.getSize();
+    const FiniteElemSüace *feSpace = vec.getFeSpace(0);
+    Mesh *mesh = feSpace->getMesh();
+    
+
+    // === First, create a continous leaf element index on each subdomain ===
+
+    std::set<int> leafElIndex;
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+    while (elInfo) {
+      leafElIndex.insert(elInfo->getElement()->getIndex());
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    map<int, int> mapElIndex;
+    int c = nLeafEls;
+    for (std::set<int>::iterator it = leafElIndex.begin();
+	 it != leafElIndex.end(); ++it)
+      mapElIndex[*it] = nLeafEls++;
+
+
+
+    int nLevel = 1;
+    int nSubdomains = meshDistributor->getMpiSize();
+    int length = 1;
+    int nSubPerProc = 1;
+    MPI_Fint c2f = MPI_Comm_c2f(MPI::COMM_WORLD);
+    int verboseLevel = 2;
+    int numbase = 0;
+
+    bddcml_init(&nLevel, &nSubdomains, &length, &nSubPerProc, 
+		&c2f, &verboseLevel, &numbase);
+
+    // global number of elements
+    int nelem = mesh->getNumberOfLeaves();
+    mpi::globalAdd(nelem);
+
+    // global number of nodes
+    int nnod = meshDistributor->getNumberOverallDofs(feSpace);
+
+    // global number of dofs
+    int ndof = nnod * nComponents;
+
+    // space dimenstion
+    int ndim = 2;
+
+    // mesh dimension
+    int meshdim = 2;
+
+    // global indes of subdomain
+    int isub = meshDistributor->getMpiRank();
+
+    // local number of elements
+    int nelems = nLeafEls;
+
+    // local number of nodes
+    int nnods = feSpace->getUsedSize();
+
+    // local number of dofs
+    int ndofs = nnods * nComponents;
+
+    // Length of array inet
+    int linet = nelems * 3;
+    
+    // Local array with indices of nodes on each element
+    int inet[linet];
+    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+    while (elInfo) {
+      int localElIndex = mapElIndex[elInfo->getElement()->getIndex()];
+      for (int i = 0; i < 3; i++)
+	inet[localElIndex * 3 + i] = elInfo->getElement()->getDof(i, 0);
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+
+    // local array of number of nodes per element
+    int nnet[nelems];
+    for (int i = 0; i < nelems; i++)
+      nnet[i] = 3;
+
+    // local array with number of DOFs per node.
+    int nndf[nnod];
+    for (int i = 0; i < nnod; i++)
+      nndf[i] = nComponents;
+
+    // array of indices of subdomain nodes in global numbering
+    int isngn[nnod];
+    for (int i = 0; i < nnod; i++)
+      nnod[i] = meshDistributor->mapLocalToGlobal(feSpace, i);
+
+    // array of indices of subdomain variables in global numbering
+    int isvgvn[ndof];
+
+    /*
+      bddcml_upload_subdomain_data(
+      &nelem,
+      &nnod,
+      &ndof,
+      &ndim,
+      &meshdim,
+      &isub,
+      &nelems,
+      &nnods,
+      &ndofs,
+      inet,
+      &linet,
+      nnet,
+      &nelems,
+      nndf,
+      &nnods,
+      isngn,
+      &nnods,
+
+      );
+    */
+
+    bddcml_finalize();
+  }
+
+
+  void BddcMlSolver::destroyMatrixData()
+  {
+    FUNCNAME("BddcMlSolver::destroyMatrixData()");
+  }
+
+#endif
+
+}

AMDiS/src/parallel/BddcMlSolver.h

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ==  http://www.amdis-fem.org                                              ==
+// ==                                                                        ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file BddcMlSolver.h */
+
+#ifndef AMDIS_BDDC_ML_SOLVER_H
+#define AMDIS_BDDC_ML_SOLVER_H
+
+#include "AMDiS_fwd.h"
+#include "parallel/PetscSolver.h"
+
+namespace AMDiS {
+
+  using namespace std;
+
+#ifdef HAVE_BDDC_ML
+
+  class BddcMlSolver : public PetscSolver
+  {
+  public:
+    BddcMlSolver()
+      : PetscSolver()
+    {}
+
+    void fillPetscMatrix(Matrix<DOFMatrix*> *mat);
+
+    void fillPetscRhs(SystemVector *vec);
+
+    void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo);    
+
+    void destroyMatrixData();
+  };
+
+#endif
+
+}
+
+#endif
+

AMDiS/src/parallel/PetscProblemStat.cc

@@ -17,6 +17,7 @@
 #include "parallel/PetscProblemStat.h"
 #include "parallel/PetscSolver.h"
 #include "parallel/MpiHelper.h"
+#include "parallel/BddcMlSolver.h"
 
 namespace AMDiS {
 
@@ -41,6 +42,12 @@ namespace AMDiS {
       petscSolver = new PetscSolverGlobalBlockMatrix();
     } else if (name == "petsc" || name == "") {
       petscSolver = new PetscSolverGlobalMatrix();
+    } else if (name == "bddcml") {
+#ifdef HAVE_BDDC_ML
+      petscSolver = new BddcMlSolver();
+#else
+      ERROR_EXIT("AMDiS was compiled without BDDC-ML support!\n");
+#endif
     } else {
       ERROR_EXIT("No parallel solver %s available!\n", name.c_str());
     }