Fixed merge with parallelization branch.

AMDiS/bin/Makefile.am

@@ -18,6 +18,7 @@ if USE_PARALLEL_AMDIS
   $(PARALLEL_DIR)/ParallelProblem.h $(PARALLEL_DIR)/ParallelProblem.cc \
   $(PARALLEL_DIR)/ParMetisPartitioner.h $(PARALLEL_DIR)/ParMetisPartitioner.cc \
   $(PARALLEL_DIR)/PartitionElementData.h \
+  $(PARALLEL_DIR)/PollutionError.h $(PARALLEL_DIR)/PollutionError.cc
   PARALLEL_INCLUDES = -I$(MPI_DIR)/include -I$(PARMETIS_DIR)
   libamdis_la_CXXFLAGS += -DHAVE_PARALLEL_AMDIS=1
 else

AMDiS/bin/Makefile.in

@@ -86,14 +86,14 @@ am__libamdis_la_SOURCES_DIST = $(PARALLEL_DIR)/ParallelDomainBase.h \
 	$(PARALLEL_DIR)/ParallelProblem.cc \
 	$(PARALLEL_DIR)/ParMetisPartitioner.h \
 	$(PARALLEL_DIR)/ParMetisPartitioner.cc \
-	$(PARALLEL_DIR)/PartitionElementData.h PARALLEL_INCLUDES = \
-	-I$(MPI_DIR)/include -I$(PARMETIS_DIR) \
-	$(SOURCE_DIR)/DOFIndexed.h $(SOURCE_DIR)/DOFIndexed.cc \
-	$(SOURCE_DIR)/GNUPlotWriter.h $(SOURCE_DIR)/GNUPlotWriter.cc \
-	$(SOURCE_DIR)/VertexVector.h $(SOURCE_DIR)/VertexVector.cc \
-	$(SOURCE_DIR)/PeriodicBC.h $(SOURCE_DIR)/PeriodicBC.cc \
-	$(SOURCE_DIR)/Recovery.h $(SOURCE_DIR)/Recovery.cc \
-	$(SOURCE_DIR)/RecoveryEstimator.h \
+	$(PARALLEL_DIR)/PartitionElementData.h \
+	$(PARALLEL_DIR)/PollutionError.h \
+	$(PARALLEL_DIR)/PollutionError.cc $(SOURCE_DIR)/DOFIndexed.h \
+	$(SOURCE_DIR)/DOFIndexed.cc $(SOURCE_DIR)/GNUPlotWriter.h \
+	$(SOURCE_DIR)/GNUPlotWriter.cc $(SOURCE_DIR)/VertexVector.h \
+	$(SOURCE_DIR)/VertexVector.cc $(SOURCE_DIR)/PeriodicBC.h \
+	$(SOURCE_DIR)/PeriodicBC.cc $(SOURCE_DIR)/Recovery.h \
+	$(SOURCE_DIR)/Recovery.cc $(SOURCE_DIR)/RecoveryEstimator.h \
 	$(SOURCE_DIR)/RecoveryEstimator.cc \
 	$(SOURCE_DIR)/ResidualEstimator.h \
 	$(SOURCE_DIR)/ResidualEstimator.cc $(SOURCE_DIR)/Cholesky.h \
@@ -237,6 +237,7 @@ am__libamdis_la_SOURCES_DIST = $(PARALLEL_DIR)/ParallelDomainBase.h \
 @USE_PARALLEL_AMDIS_TRUE@	libamdis_la-ConditionalMarker.lo \
 @USE_PARALLEL_AMDIS_TRUE@	libamdis_la-ParallelProblem.lo \
 @USE_PARALLEL_AMDIS_TRUE@	libamdis_la-ParMetisPartitioner.lo \
+@USE_PARALLEL_AMDIS_TRUE@	libamdis_la-PollutionError.lo \
 @USE_PARALLEL_AMDIS_TRUE@	$(am__objects_1)
 am_libamdis_la_OBJECTS = $(am__objects_2) libamdis_la-DOFIndexed.lo \
 	libamdis_la-GNUPlotWriter.lo libamdis_la-VertexVector.lo \
@@ -471,10 +472,11 @@ libamdis_la_CXXFLAGS = $(am__append_1) $(am__append_3) $(am__append_5) \
 @USE_PARALLEL_AMDIS_TRUE@	$(PARALLEL_DIR)/ParMetisPartitioner.h \
 @USE_PARALLEL_AMDIS_TRUE@	$(PARALLEL_DIR)/ParMetisPartitioner.cc \
 @USE_PARALLEL_AMDIS_TRUE@	$(PARALLEL_DIR)/PartitionElementData.h \
-@USE_PARALLEL_AMDIS_TRUE@	PARALLEL_INCLUDES = \
-@USE_PARALLEL_AMDIS_TRUE@	-I$(MPI_DIR)/include \
-@USE_PARALLEL_AMDIS_TRUE@	-I$(PARMETIS_DIR) $(am__append_2)
+@USE_PARALLEL_AMDIS_TRUE@	$(PARALLEL_DIR)/PollutionError.h \
+@USE_PARALLEL_AMDIS_TRUE@	$(PARALLEL_DIR)/PollutionError.cc \
+@USE_PARALLEL_AMDIS_TRUE@	$(am__append_2)
 @USE_PARALLEL_AMDIS_FALSE@PARALLEL_INCLUDES = 
+@USE_PARALLEL_AMDIS_TRUE@PARALLEL_INCLUDES = -I$(MPI_DIR)/include -I$(PARMETIS_DIR)
 TEMPLATE_INCLUDES = -I../lib/mtl4 
 INCLUDES = $(AMDIS_INCLUDES) $(PARALLEL_INCLUDES) $(TEMPLATE_INCLUDES)
 libamdis_la_SOURCES = \
@@ -784,6 +786,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Parametric.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-PeriodicBC.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-PngWriter.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-PollutionError.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-PovrayWriter.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ProblemInstat.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ProblemInterpolScal.Plo@am__quote@
@@ -893,6 +896,13 @@ libamdis_la-ParMetisPartitioner.lo: $(PARALLEL_DIR)/ParMetisPartitioner.cc
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ParMetisPartitioner.lo `test -f '$(PARALLEL_DIR)/ParMetisPartitioner.cc' || echo '$(srcdir)/'`$(PARALLEL_DIR)/ParMetisPartitioner.cc
 
+libamdis_la-PollutionError.lo: $(PARALLEL_DIR)/PollutionError.cc
+@am__fastdepCXX_TRUE@	if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-PollutionError.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-PollutionError.Tpo" -c -o libamdis_la-PollutionError.lo `test -f '$(PARALLEL_DIR)/PollutionError.cc' || echo '$(srcdir)/'`$(PARALLEL_DIR)/PollutionError.cc; \
+@am__fastdepCXX_TRUE@	then mv -f "$(DEPDIR)/libamdis_la-PollutionError.Tpo" "$(DEPDIR)/libamdis_la-PollutionError.Plo"; else rm -f "$(DEPDIR)/libamdis_la-PollutionError.Tpo"; exit 1; fi
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='$(PARALLEL_DIR)/PollutionError.cc' object='libamdis_la-PollutionError.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-PollutionError.lo `test -f '$(PARALLEL_DIR)/PollutionError.cc' || echo '$(srcdir)/'`$(PARALLEL_DIR)/PollutionError.cc
+
 libamdis_la-DOFIndexed.lo: $(SOURCE_DIR)/DOFIndexed.cc
 @am__fastdepCXX_TRUE@	if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-DOFIndexed.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-DOFIndexed.Tpo" -c -o libamdis_la-DOFIndexed.lo `test -f '$(SOURCE_DIR)/DOFIndexed.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/DOFIndexed.cc; \
 @am__fastdepCXX_TRUE@	then mv -f "$(DEPDIR)/libamdis_la-DOFIndexed.Tpo" "$(DEPDIR)/libamdis_la-DOFIndexed.Plo"; else rm -f "$(DEPDIR)/libamdis_la-DOFIndexed.Tpo"; exit 1; fi

AMDiS/src/PollutionError.cc

+#include "PollutionError.h"
+#include "mpi.h"
+
+namespace AMDiS
+{
+
+  PollutionErrorKonvex::BoundaryTraverseStack::BoundaryTraverseStack(Mesh *mesh, 
+								     PartitionStatus ps_)
+  {
+    TraverseStack myStack;
+    ElInfo *swap = 
+      myStack.traverseFirst(mesh,-1, 
+			    Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_NEIGH);
+    ps = ps_;
+
+    WorldVector<double>* centerWorld;
+    DimVec<double> center = DimVec<double>(2, DEFAULT_VALUE, 1.0 / 3.0);
+    while (swap != NULL) {
+      if (proofConditions(swap)) {
+	// if conditions are true, push the midpoint to the end of the list. a better 
+	// algorithm should use the boundary edge
+	centerWorld = new WorldVector<double>();
+	swap->coordToWorld(center, *centerWorld);
+	push_back(centerWorld);
+      }
+      swap = myStack.traverseNext(swap);
+    }
+  }
+	
+  bool AMDiS::PollutionErrorKonvex::BoundaryTraverseStack::proofConditions(ElInfo *elInfo) 
+  {
+    if (elInfo == NULL) 
+      return false;
+
+    Element *curEl = elInfo->getElement();
+    PartitionElementData *ped = 
+      static_cast<PartitionElementData*>(curEl->getElementData(PARTITION_ED));
+
+    TEST_EXIT_DBG(ped != NULL)("No PartitionElementData found");    
+
+    if (ped != NULL && ped->getPartitionStatus() == OVERLAP) {
+      //check neighbour
+      int nrNeigh = curEl->getGeo(NEIGH);
+      for (int i = 0; i < nrNeigh; i++) {
+	Element* curNeigh = elInfo->getNeighbour(i);
+	if (curNeigh != NULL) { 		    
+	  ped = 
+	    static_cast<PartitionElementData*>(curNeigh->getElementData(PARTITION_ED));
+	  if (ped != NULL && ped->getPartitionStatus() == ps)
+	    return true;
+	}
+      }
+    }
+    return false;
+  }
+
+  PollutionErrorKonvex::PollutionErrorKonvex(std::string name, int r)
+    : ResidualEstimator(name,r),
+      center(2,DEFAULT_VALUE, 1.0 / 3.0)
+  {
+    FUNCNAME("PollutionErrorKonvexEstimator::PollutionErrorKonvexEstimator()");
+    GET_PARAMETER(0, name + "->C1p", "%f", &C1p);
+    GET_PARAMETER(0, name + "->C2p", "%f", &C2p);   
+  }
+
+  void PollutionErrorKonvex::init(double timestep)  
+  {
+    FUNCNAME("PollutionErrorKonvex::init");
+    ResidualEstimator::init(timestep);
+
+    // set the distance (ie. width of the overlapped area)
+    // approximate the overlapped area as rectangle. then set d=max(edge1, edge2)
+    // retrieve all elements at the boundary overlapped -> out
+    
+    bts = new BoundaryTraverseStack(mesh, IN);	
+  }
+
+  void PollutionErrorKonvex::estimateElement(ElInfo *elInfo)  
+  {
+    FUNCNAME("PollutionErrorKonvex::estimateElement()");
+
+    // switch between inner and outer
+    // inner part --> use H_1Norm
+    // outer part --> use L_2Norm
+
+    PartitionElementData* ped = 
+      static_cast<PartitionElementData*>(elInfo->getElement()->getElementData(PARTITION_ED));
+    TEST_EXIT_DBG(ped != NULL)("No partition Elementdata found");
+
+    PartitionStatus ps = ped->getPartitionStatus();
+    TEST_EXIT_DBG(ps != UNDEFINED)("Elementstatus undefined");
+
+    // save old values
+    double C0 = ResidualEstimator::C0;
+    double C1 = ResidualEstimator::C1;
+
+    if (ps == OUT) {
+      //outer part
+      ResidualEstimator::norm = L2_NORM;		
+      ResidualEstimator::C0 *= C2p;
+
+      //compute distance
+      double d = getDistance(elInfo);
+      ResidualEstimator::C1 *= C2p / d;
+      
+    } else {
+      //inner part
+      ResidualEstimator::norm = H1_NORM;
+      ResidualEstimator::C0 *= C1p;
+      ResidualEstimator::C1 *= C1p;
+    }
+
+    ResidualEstimator::estimateElement(elInfo);
+    //set old values
+    ResidualEstimator::C0 = C0;
+    ResidualEstimator::C1 = C1;
+  }
+
+  double PollutionErrorKonvex::getDistance(ElInfo* element)  
+  {
+    if (bts==NULL)
+      MSG("bts is null");
+    TEST_EXIT_DBG(bts!=NULL)("no bts");
+    double distance = 0.0;
+    double curDistance = 1.0;
+    //get center in world coord
+    WorldVector<double> centerWorld;
+    element->coordToWorld(center,centerWorld);
+
+    //start traversal over all Midpoints at the boundary overlap <-> in      
+    for (BoundaryTraverseStack::iterator myIt = bts->begin(); 
+	 myIt != bts->end(); myIt++) {
+      curDistance = 0;	
+      for (int i = 0;i < dim;i++)
+	curDistance += ((*(*myIt))[i]-centerWorld[i]) * ((*(*myIt))[i] - centerWorld[i]);
+
+      distance = max(curDistance, distance);
+    }
+    return distance;
+  }
+  
+  void PollutionErrorKonvex::exit(bool output = true) 
+  {
+    ResidualEstimator::exit(output);
+    delete bts;
+  }
+}
+

AMDiS/src/PollutionError.h

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ============================================================================
+// ==                                                                        ==
+// ==  TU Dresden                                                            ==
+// ==                                                                        ==
+// ==  Institut fr Wissenschaftliches Rechnen                               ==
+// ==  Zellescher Weg 12-14                                                  ==
+// ==  01069 Dresden                                                         ==
+// ==  germany                                                               ==
+// ==                                                                        ==
+// ============================================================================
+// ==                                                                        ==
+// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
+// ==                                                                        ==
+// ============================================================================
+
+/** \file PollutionError.h */
+
+/** \defgroup Estimator Estimator module
+ * @{ <img src="estimator.png"> @}
+ */
+
+#ifndef AMDIS_POLLUTIONERROR_KONVEX_H
+#define AMDIS_POLLUTIONERROR_KONVEX_H
+
+#include <list>
+#include <time.h>
+#include "ResidualEstimator.h"
+#include "PartitionElementData.h"
+
+namespace AMDiS {
+
+  /** \ingroup Estimator
+   * \brief Estimator for parallel solutions. 
+   *
+   * estimates the error using following formulae
+   * \f$||e_h||_{1,\Omega_0}=C_1 ||u_h||_{H_1(\Omega_1)}+
+   * \frac{C_2}{d^2}||u_h||_{L_2(\Omega\setminus\Omega_1}\f$
+   * where \f$\Omega\f$ is the complete area of interest and \f$\Omega_1\f$ the area of 
+   * this process  with overlapping. \f$\Omega_0\$ is the area of this process without 
+   * overlapping,\f$d\$ contains the distance 
+   * dist(\f$\partial \Omega_0,\partial \Omega_1\f$)
+   */
+  class PollutionErrorKonvex : public ResidualEstimator 
+  {
+  public:
+    class Creator : public EstimatorCreator 
+    {
+    public:
+      Creator() : EstimatorCreator() {}
+
+      virtual ~Creator() {}
+
+      /// Returns a new ODirSolver object.
+      Estimator* create() 
+      { 
+	return new PollutionErrorKonvex(name, row);
+      }
+    };
+    
+    PollutionErrorKonvex(std::string name, int r);
+
+    //sets the overlapping width for 2d-elements
+    virtual void init(double timestep);
+
+    virtual void estimateElement(ElInfo *elInfo);
+
+    virtual void exit(bool output);
+
+  protected:
+    /// constant in front of the inner residual   
+    double C1p;
+    
+    /// constant in front of the outer residual    
+    double C2p;
+
+    ///traverse-class for edges and there normals in 2d
+    class BoundaryTraverseStack: public std::list< WorldVector<double>* >
+    {
+    public:
+      BoundaryTraverseStack(Mesh *mesh, PartitionStatus p);			
+      
+    private:
+      PartitionStatus ps;
+      bool proofConditions(ElInfo *elInfo);
+    };
+
+    /// computes the squared distance of element el to the boundary of \f$\Omega_0$\f$
+    double getDistance(ElInfo* el);
+
+    /// stack to traverse over the boundary of \f$\Omega_0$\f$
+    BoundaryTraverseStack* bts;
+
+    /// barycentric center
+    DimVec<double> center;
+  };
+}
+#endif