merge parallel multimesh branch to the trunk

AMDiS/src/AMDiS_fwd.h

@@ -163,6 +163,10 @@ namespace AMDiS {
   template<typename T>                                 class WorldVector;
   template<typename T>                                 class WorldMatrix;
   template<typename T>                                 class VectorOfFixVecs;
+  
+  namespace detail {
+    template <class P> class CouplingProblemStat;
+  }
 
 
   typedef mtl::dense2D<double>              ElementMatrix;

AMDiS/src/Assembler.cc

@@ -53,7 +53,6 @@ namespace AMDiS {
   Assembler::~Assembler()
   {}
 
-
   void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
 					 ElementMatrix& userMat,
 					 double factor)
@@ -64,19 +63,22 @@ namespace AMDiS {
     Element *el = elInfo->getElement();
 
     if (el != lastMatEl || !operat->isOptimized()) {
-      initElement(elInfo, elInfo);
+      initElement(elInfo);
 
       if (rememberElMat)
 	set_to_zero(elementMatrix);
 
       lastMatEl = el;
     } else {
+      // Only possible in single mesh case when one operator 
+      // is used more than twice?
       if (rememberElMat) {
-	userMat += factor * elementMatrix;
+	if (&userMat != &elementMatrix)
+	  userMat += factor * elementMatrix;
 	return;
       }
     }
-
+    
     ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
 
     if (secondOrderAssembler)
@@ -91,8 +93,7 @@ namespace AMDiS {
     if (rememberElMat && &userMat != &elementMatrix)
       userMat += factor * elementMatrix;
   }
-
-
+  
   void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
 					 const ElInfo *colElInfo,
 					 const ElInfo *smallElInfo,
@@ -101,115 +102,68 @@ namespace AMDiS {
 					 ElementMatrix& userMat,
 					 double factor)
   {
-    if (remember && (factor != 1.0 || operat->uhOld))
-      rememberElMat = true;
-  
-    Element *el = smallElInfo->getElement();   
+    Element *el = smallElInfo->getElement();
+    // initElement and calculateElementMatrix always need to be
+    // recalculated since the smaller ElInfo can be different
     lastVecEl = lastMatEl = NULL;
    
     if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
       initElement(smallElInfo, largeElInfo);
 
-    if (el != lastMatEl || !operat->isOptimized()) {
-      if (rememberElMat)
-	set_to_zero(elementMatrix);
 
+    if (el != lastMatEl || !operat->isOptimized())
       lastMatEl = el;
-    } else {
-      if (rememberElMat) {
+    else {
+      ERROR_EXIT("Never reached because lastVecEl = lastMatEl = NULL;\n");
+      if (&userMat != &elementMatrix)
 	userMat += factor * elementMatrix;
-	return;
-      }
+      return;
     }
- 
-    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
-    
-    if (secondOrderAssembler) {
-      // calculate element matrices always on smallest element
-      secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
 
-      // smallElInfo stores refinement-relation to largeElInfo
-      ElementMatrix &m = 
-	smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());  // muste be moved to next if-else block when generalized for multiple polynomial degrees
-      
-      if (!rowColFeSpaceEqual) {
-	if (smallElInfo == colElInfo)
-	  tmpMat = m * mat;	
-	else
-	  tmpMat = mat * trans(m);
-	
-	mat = tmpMat;
-      }
-    }
+    // In multimesh case, elementMatrix is always be recalculated, because each time we 
+    // need to do transformation here.
+    set_to_zero(elementMatrix);
+    
+    ElementMatrix& mat = elementMatrix;
 
-    if (firstOrderAssemblerGrdPsi) {
+    if (secondOrderAssembler)
+      secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
+    if (firstOrderAssemblerGrdPsi)
       firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
-
-      if (!rowColFeSpaceEqual) {
-	if (largeElInfo == rowElInfo) {
-	  ElementMatrix &m = 
-	    smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-	  
-	  tmpMat = m * mat;
-	} else {
-	  ElementMatrix &m = 
-	    smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-	  
-	  tmpMat = mat * trans(m);
-	}
-	
-	mat = tmpMat;
-      }
-    }
-
-    if (firstOrderAssemblerGrdPhi) {
+    if (firstOrderAssemblerGrdPhi)
       firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
-
-      if (!rowColFeSpaceEqual) {
-	if (largeElInfo == colElInfo) {
-	  ElementMatrix &m = 
-	    smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-	  
-	  tmpMat = mat * trans(m);
-	} else {
-	  ElementMatrix &m = 
-	    smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-	  
-	  tmpMat = m * mat;	
-	}
-	
-	mat = tmpMat;
-      }
-    }
-
-    if (zeroOrderAssembler) {
+    if (zeroOrderAssembler)
       zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
+    
+    
 
-      if (!rowColFeSpaceEqual) {
-	ElementMatrix &m = 
-	  smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-	
-	if (smallElInfo == colElInfo)
-	  tmpMat = m * mat;
-	else 
-	  tmpMat = mat * trans(m);
-	
-	mat = tmpMat;
-      }
-    }
+    ElementMatrix &m = 
+      smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
 
-    if (rememberElMat && &userMat != &elementMatrix)
-      userMat += factor * elementMatrix;       
+    if (!rowColFeSpaceEqual) {
+      if (smallElInfo == colElInfo) 
+	tmpMat = m * mat;
+      else
+	tmpMat = mat * trans(m);
+      
+      mat = tmpMat;
+    } else if (smallElInfo == colElInfo) {
+      tmpMat = m * mat * trans(m);
+      mat = tmpMat;
+    }
+    
+    // This is for the call in matVecAssemble
+    if (&userMat != &elementMatrix)
+      userMat += factor * elementMatrix;   
   }
 
-
   void Assembler::calculateElementVector(const ElInfo *elInfo, 
 					 ElementVector& userVec,
 					 double factor)
   {
     if (remember && factor != 1.0)
       rememberElVec = true;
-
+    
     Element *el = elInfo->getElement();
 
     if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
@@ -221,6 +175,8 @@ namespace AMDiS {
 	
       lastVecEl = el;
     } else {
+      // Only possible in single mesh case when one operator 
+      // is used more than twice at dof vector?
       if (rememberElVec) {
 	userVec += factor * elementVector;
 	return;
@@ -246,7 +202,6 @@ namespace AMDiS {
       userVec += factor * elementVector;    
   }
 
-
   void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
 					 const ElInfo *auxElInfo,
 					 const ElInfo *smallElInfo,
@@ -256,36 +211,33 @@ namespace AMDiS {
   {
     FUNCNAME("Assembler::calculateElementVector()");
 
-    if (remember && factor != 1.0)
-      rememberElVec = true;
-
     Element *el = mainElInfo->getElement();
+    // initElement and calculateElementVector always need to be
+    // recalculated since the smaller ElInfo can be different
+    lastVecEl = lastMatEl = NULL;
 
     if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
       initElement(smallElInfo, largeElInfo);
-   
-    if (el != lastVecEl || !operat->isOptimized()) {
-      if (rememberElVec)
-	set_to_zero(elementVector);
-
+    
+    if (el != lastVecEl || !operat->isOptimized())
       lastVecEl = el;
-    } else {
-      if (rememberElVec) {
-	userVec += factor * elementVector;
-	return;
-      }
+    else {
+      ERROR_EXIT("Never reached because lastVecEl = lastMatEl = NULL;\n");
+      userVec += factor * elementVector;
+      return;
     }
-    ElementVector& vec = rememberElVec ? elementVector : userVec;
+    
+    set_to_zero(elementVector);
+    ElementVector& vec = elementVector;
 
     if (operat->uhOld && remember) {
+      
       if (smallElInfo->getLevel() == largeElInfo->getLevel())
 	matVecAssemble(auxElInfo, vec);
-      else
+      else 
 	matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec);      
 
-      if (rememberElVec)
-	userVec += factor * elementVector;      
-
+      userVec += factor * elementVector;
       return;
     } 
 
@@ -306,13 +258,12 @@ namespace AMDiS {
       }      
     }
 
-    if (rememberElVec)
-      userVec += factor * elementVector;    
+    userVec += factor * elementVector;    
   }
 
-
   void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
   {
+
     Element *el = elInfo->getElement(); 
     ElementVector uhOldLoc(operat->uhOld->getFeSpace() == rowFeSpace ? 
 			   nRow : nCol);
@@ -323,7 +274,8 @@ namespace AMDiS {
       calculateElementMatrix(elInfo, elementMatrix);
     }
 
-    vec += elementMatrix*uhOldLoc;
+
+    vec += elementMatrix * uhOldLoc;
   }
 
 
@@ -331,7 +283,7 @@ namespace AMDiS {
 				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
 				 ElementVector& vec)
   {
-    FUNCNAME("Assembler::matVecAssemble()");
+     FUNCNAME("Assembler::matVecAssemble()");
 
     TEST_EXIT(rowFeSpace->getBasisFcts() == colFeSpace->getBasisFcts())
       ("Works only for equal basis functions for different components!\n");
@@ -343,7 +295,7 @@ namespace AMDiS {
       usedEl = mainElInfo->getElement();
     else
       ERROR("Mesh is incorrect.\n");
-    
+
     const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
     int nBasFcts = basFcts->getNumber();
     ElementVector uhOldLoc(nBasFcts);
@@ -353,10 +305,11 @@ namespace AMDiS {
     if (mainElInfo->getElement() != lastMatEl) {
       set_to_zero(elementMatrix);
       calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
-			     rowFeSpace == operat->uhOld->getFeSpace(), elementMatrix);  
+			     rowFeSpace == operat->uhOld->getFeSpace(), elementMatrix);    
     }
 
     vec += elementMatrix * uhOldLoc;
+
   }
 
 
@@ -472,4 +425,4 @@ namespace AMDiS {
     checkQuadratures();
   }
 
-}
+}
\ No newline at end of file

AMDiS/src/BoundaryObject.cc

@@ -56,6 +56,8 @@ namespace AMDiS {
     FUNCNAME("BoundaryObject::computeReverseMode()");
 
     bool reverseMode = false;
+    
+    TEST_EXIT_DBG(obj0.el && obj1.el)("BoundaryObject without element pointer.\n");
 
     switch (feSpace->getMesh()->getDim()) {
     case 2:

AMDiS/src/CoarseningManager3d.cc

@@ -99,7 +99,7 @@ namespace AMDiS {
 
 #if HAVE_PARALLEL_DOMAIN_AMDIS
     vector<FixRefinementPatch::EdgeInEl> refineEdges;
-    FixRefinementPatch::getOtherEl(stack, refineEdges);
+    FixRefinementPatch::getOtherEl(mesh, stack, refineEdges);
 
     // === If the refinement edge must be fixed, add also the other part of this ===
     // === edge to the refinement patch.                                         ===
@@ -110,6 +110,7 @@ namespace AMDiS {
       return;
 
       Element *otherEl = refineEdges[0].first;
+      TEST_EXIT_DBG(otherEl->getMesh() == mesh)("Something is wrong.\n");
       
       TraverseStack stack2;
       ElInfo *elInfo2 = 

AMDiS/src/CouplingIterationInterface.h

@@ -40,7 +40,7 @@ namespace AMDiS {
    * problem-list, that is filled by addProblem. Alternatively one can a access
    * each problem by an unique name.
    */
-  class CouplingIterationInterface : public ProblemIterationInterface
+  class CouplingIterationInterface : public virtual ProblemIterationInterface
   {
   public:
     virtual ~CouplingIterationInterface() {}

AMDiS/src/CouplingProblemStat.h

@@ -31,7 +31,7 @@
 #include "AMDiS_fwd.h"
 #include "ProblemStat.h"
 #include "Initfile.h"
-#include <boost/lexical_cast.hpp>
+#include "utility/to_string.hpp"
 
 namespace AMDiS {
 
@@ -43,21 +43,26 @@ namespace AMDiS {
     * computations. 
     */
     template<typename ProblemStatType>
-    class CouplingProblemStat
+    class CouplingProblemStat : public ProblemStatSeq
     {
+    protected:
+      typedef ProblemStatSeq  super;
+      
+      using super::nComponents;
+      using super::meshes;
+      using super::nMeshes;
+      using super::feSpaces;
+      using super::name;
+      using super::refinementManager;
+      using super::coarseningManager;
+      
     public:
       /// Constructor
       CouplingProblemStat(std::string name_)
-	: name(name_),
-	  nComponents(0),
-	  nMeshes(0),
-	  refinementManager(NULL),
-	  coarseningManager(NULL)
+	: super(name_),
+	  dim(-1)
       {}
 
-      /// Destructor
-      virtual ~CouplingProblemStat() {}
-
       /// add problem by number
       virtual void addProblem(ProblemStatType* prob)
       {
@@ -71,112 +76,193 @@ namespace AMDiS {
 			      Flag adoptFlag = INIT_NOTHING)
       { FUNCNAME("CouplingProblemStat::initialize()");
 
-	// create one refinement-/coarseningmanager for all problems
-	if (refinementManager != NULL && coarseningManager != NULL) { 
-	  WARNING("refinement-/coarseningmanager already created\n");
-	} else {
-	  if (!adoptProblem)
-	    createRefCoarseManager();
-	  else {
-	    refinementManager = adoptProblem->getRefinementManager();
-	    coarseningManager = adoptProblem->getCoarseningManager();
-	  }
+	super::initialize(initFlag - INIT_MESH);
+	
+	const Flag DEFAULT_INIT = (INIT_FE_SPACE | INIT_MESH | CREATE_MESH | INIT_SYSTEM | INIT_SOLVER | INIT_ESTIMATOR | INIT_MARKER | INIT_FILEWRITER);
+	for (size_t p = 0; p < problems.size(); ++p) {
+	  problems[p]->initialize(initFlag - DEFAULT_INIT);
 	}
-    
-	if (refinementManager == NULL || coarseningManager == NULL) 
-	  WARNING("no refinement-/coarseningmanager created\n");
 
-	// create Meshes and FeSpaces
+	for (size_t i = 0; i < meshes.size(); i++) {
+	  int globalRefinements = 0;
+
+	  // If AMDiS is compiled for parallel computations, the global refinements are
+	  // ignored here. Later, each rank will add the global refinements to its
+	  // private mesh.
+  #ifndef HAVE_PARALLEL_DOMAIN_AMDIS
+	  Parameters::get(meshes[i]->getName() + "->global refinements",
+			  globalRefinements);
+  #endif
 
+	  bool initMesh = initFlag.isSet(INIT_MESH);
 
+	  // Initialize the meshes if there is no serialization file.
+	  if (initMesh && meshes[i] && !(meshes[i]->isInitialized())) {
+	    meshes[i]->initialize();
+	    refinementManager->globalRefine(meshes[i], globalRefinements);
+	  }
+	}
+      }
+      
+      
+      /// Used in \ref initialize().
+      virtual void createMesh() override
+      {
 	// all problems must have the same dimension (?)
-	int dim = 0;
+	dim = 0;
 	Parameters::get(name + "->dim", dim);
 	TEST_EXIT(dim)("No problem dimension specified for \"%s->dim\"!\n",
 		      name.c_str());
 
-	std::map<std::string, Mesh*> meshByName;
+	std::map<std::pair<std::string, int>, Mesh*> meshByName;  // (name, refSet) --> Mesh*
+	typedef std::map<std::pair<std::string, int>, Mesh*>::iterator MeshIterator;
+
 	std::vector<std::set<Mesh*> > meshesForProblems(problems.size());
-	std::map<std::pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
+	
 	for (size_t i = 0; i < problems.size(); ++i) {
 	  TEST_EXIT(problems[i])("problem[%d] does not exist!\n",i);
-	  for (size_t j = 0; j < problems[i]->getNumComponents(); j++) {
-	    // mesh
-
+	  
+	  int nComponents = problems[i]->getNumComponents();
+	  
+	  int nAddComponents = 0;
+	  Parameters::get(problems[i]->getName() + "->additional components", nAddComponents);
+	  
+	  problems[i]->componentMeshes.resize(nComponents + nAddComponents);
+	
+	  for (size_t j = 0; j < nComponents + nAddComponents; j++) {
+	    // name of the mesh
 	    std::string meshName("");
 	    Parameters::get(problems[i]->getName() + "->mesh", meshName);
 	    TEST_EXIT(meshName != "")("No mesh name specified for \"%s->mesh\"!\n",
 				      problems[i]->getName().c_str());
-
-	    if (meshByName.find(meshName) == meshByName.end()) {
+				      
+	    // dimension of the mesh
+	    int mesh_dim = 0;
+	    Parameters::get(problems[i]->getName() + "->dim", mesh_dim);
+	    TEST_EXIT(dim == mesh_dim)("Mesh-dimension must be the same for all problems!\n");
+		
+	    // refinement set (optional)
+	    int refSet = 0;
+	    Parameters::get(problems[i]->getName() + "->refinement set[" + to_string(j) + "]", refSet);
+
+	    // create a new Mesh only if not already created for other problem
+	    Mesh* componentMesh;
+	    MeshIterator meshIt = meshByName.find(std::make_pair(meshName, refSet));
+	    if (meshIt == meshByName.end()) {
 	      Mesh *newMesh = new Mesh(meshName, dim);
-	      meshByName[meshName] = newMesh;
 	      meshes.push_back(newMesh);
-	      meshesForProblems[i].insert(newMesh);
+	      meshByName[std::make_pair(meshName, refSet)] = newMesh;	    
+	      componentMesh = newMesh;
 	      nMeshes++;
-	    } else
-	      meshesForProblems[i].insert(meshByName[meshName]);
+	    } else {
+	      componentMesh = meshIt->second;
+	    }
+	    problems[i]->componentMeshes[j] = componentMesh;
+	  }
+	  
+	  // copy unqiue set of meshes to problem[i]->meshes
+	  std::set<Mesh*> uniqueMeshes;
+	  for (size_t j = 0; j < problems[i]->componentMeshes.size(); ++j)
+	    uniqueMeshes.insert(problems[i]->componentMeshes[j]);
+	  problems[i]->meshes.clear();
+	  problems[i]->meshes.insert(problems[i]->meshes.begin(), uniqueMeshes.begin(), uniqueMeshes.end());
+	}
+      }
 
-	    problems[i]->setComponentMesh(j, meshByName[meshName]);
+      /// Used in \ref initialize().
+      virtual void createFeSpace(DOFAdmin *admin) override
+      {
+	std::vector<std::set<FiniteElemSpace const*> > feSpacesForProblems(problems.size());
+	std::map<std::pair<Mesh*, std::string>, FiniteElemSpace*> feSpaceMap;
+	
+	for (size_t p = 0; p < problems.size(); ++p) {
+	  TEST_EXIT(problems[p])("problem[%d] does not exist!\n",p);
+	  
+	  int nComponents = problems[p]->getNumComponents();
+	  
+	  int nAddComponents = 0;
+	  Parameters::get(problems[p]->getName() + "->additional components", nAddComponents);
+	  problems[p]->componentSpaces.resize(nComponents + nAddComponents, NULL);
+	  problems[p]->traverseInfo.resize(nComponents);
+	
+	  for (size_t i = 0; i < nComponents + nAddComponents; i++) {
+	  
+	    std::string componentString = "[" + to_string(i) + "]";
+	    
+	    std::string feSpaceName = "";
+	    std::string initFileStr = problems[p]->getName() + "->feSpace" + componentString;
+	    Parameters::get(initFileStr, feSpaceName);
+	    
+	    // synonym for "feSpace"
+	    if (feSpaceName.size() == 0) {
+	      initFileStr = problems[p]->getName() + "->finite element space" + componentString;
+	      Parameters::get(initFileStr, feSpaceName);
+	    }
+	    
+	    // for backward compatibility also accept the old syntax
+	    if (feSpaceName.size() == 0) {
+	      int degree = 1;
+	      initFileStr = problems[p]->getName() + "->polynomial degree" + componentString;
+	      Parameters::get(initFileStr, degree);
+	      TEST_EXIT(degree > 0)
+		("Poynomial degree in component %d must be larger than zero!\n", i);
+		
+	      feSpaceName = "Lagrange" + to_string(degree);
+	    } 
+	    
+	    if (feSpaceName.size() == 0)
+	      feSpaceName = "Lagrange1";
+
+	    if (feSpaceMap[std::make_pair(problems[p]->componentMeshes[i], feSpaceName)] == NULL) {      
+	      BasisFunctionCreator *basisFctCreator = 
+		dynamic_cast<BasisFunctionCreator*>(CreatorMap<BasisFunction>::getCreator(feSpaceName, initFileStr));
+	      TEST_EXIT(basisFctCreator)
+		("No valid basisfunction type found in parameter \"%s\"\n", initFileStr.c_str());
+	      basisFctCreator->setDim(dim);
 	      
-	    // feSpace
-	    int degree = 1;
-	    Parameters::get(problems[i]->getName() + "->polynomial degree[" + 
-			    boost::lexical_cast<std::string>(j) + "]", degree);
-
-	    if (feSpaceMap[std::pair<Mesh*, int>(meshByName[meshName], degree)] == NULL) {
-	      std::stringstream s;
-	      s << problems[i]->getName() << "->feSpace[" << j << "]";
-      
 	      FiniteElemSpace *newFeSpace = 
-		FiniteElemSpace::provideFeSpace(NULL, Lagrange::getLagrange(dim, degree),
-						meshByName[meshName], s.str());
-	      feSpaceMap[std::pair<Mesh*, int>(meshByName[meshName], degree)] = newFeSpace;
+		FiniteElemSpace::provideFeSpace(admin, basisFctCreator->create(),
+						problems[p]->componentMeshes[i], 
+						"FeSpace" + componentString + " (" + feSpaceName + ")");
+
+	      feSpaceMap[std::make_pair(problems[p]->componentMeshes[i], feSpaceName)] = newFeSpace;
 	      f