diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index 9e7be1c5adcc651d7a309c59b81e18355ca474c4..8cf0c72a5edafbc5f7a2d4040826bba65a40673b 100644
--- a/AMDiS/CMakeLists.txt
+++ b/AMDiS/CMakeLists.txt
@@ -75,6 +75,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/DOFIndexed.cc
${SOURCE_DIR}/CreatorMap.cc
${SOURCE_DIR}/ProblemInterpolScal.cc
${SOURCE_DIR}/ProblemInterpolVec.cc
+ ${SOURCE_DIR}/MacroInfo.cc
${SOURCE_DIR}/MacroReader.cc
${SOURCE_DIR}/ValueReader.cc
${SOURCE_DIR}/Projection.cc
@@ -165,10 +166,6 @@ if(ENABLE_PARMETIS)
CMAKE_FORCE_CXX_COMPILER(mpiCC "MPI C++ compiler")
include_directories( ${LIB_DIR}/ParMetis-3.1)
- SET(PARALLEL_AMDIS_SRC ${SOURCE_DIR}/ConditionalEstimator.cc ${SOURCE_DIR}/ConditionalMarker.cc
- ${SOURCE_DIR}/ParallelProblem.cc ${SOURCE_DIR}/ParMetisPartitioner.cc
- ${SOURCE_DIR}/PollutionError.cc)
- SET(COMPILEFLAGS "${COMPILEFLAGS} -DHAVE_PARALLEL_AMDIS=1")
INSTALL(FILES ${LIB_DIR}/ParMetis-3.1/parmetis.h
${LIB_DIR}/ParMetis-3.1/libparmetis.a
${LIB_DIR}/ParMetis-3.1/libmetis.a
@@ -179,6 +176,7 @@ endif(ENABLE_PARMETIS)
if(ENABLE_PARALLEL_DOMAIN)
include_directories(${PETSC_DIR}/include ${PETSC_DIR}/${PETSC_ARCH}/include)
SET(PARALLEL_DOMAIN_AMDIS_SRC
+ ${SOURCE_DIR}/parallel/ParMetisPartitioner.cc
${SOURCE_DIR}/parallel/MeshDistributor.cc
${SOURCE_DIR}/parallel/StdMpi.cc
${SOURCE_DIR}/parallel/ParallelDebug.cc
@@ -247,7 +245,7 @@ SET(COMPOSITE_FEM_SRC ${COMPOSITE_SOURCE_DIR}/CFE_Integration.cc ${COMPOSITE_SOU
include_directories(${MTL_DIR})
include_directories(${SOURCE_DIR})
-add_library(amdis SHARED ${AMDIS_SRC} ${PARALLEL_AMDIS_SRC} ${PARALLEL_DOMAIN_AMDIS_SRC})
+add_library(amdis SHARED ${AMDIS_SRC} ${PARALLEL_DOMAIN_AMDIS_SRC})
add_library(compositeFEM SHARED ${COMPOSITE_FEM_SRC})
target_link_libraries(compositeFEM amdis)
diff --git a/AMDiS/bin/Makefile.am b/AMDiS/bin/Makefile.am
index 245467093bc5ce494b9fe698569cc04bffa318ce..d4427386da8aa5c2bac8867eade65d86a554dad4 100644
--- a/AMDiS/bin/Makefile.am
+++ b/AMDiS/bin/Makefile.am
@@ -10,25 +10,13 @@ AMDIS_INCLUDES = -I$(SOURCE_DIR)
libamdis_la_CXXFLAGS =
-if USE_PARALLEL_AMDIS
- PARALLEL_AMDIS_SOURCES = \
- $(PARALLEL_DIR)/ConditionalEstimator.h $(PARALLEL_DIR)/ConditionalEstimator.cc \
- $(PARALLEL_DIR)/ConditionalMarker.h \ $(PARALLEL_DIR)/ConditionalMarker.cc \
- $(PARALLEL_DIR)/ParallelError.h $(PARALLEL_DIR)/ParallelError.hh \
- $(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
- PARALLEL_AMDIS_SOURCES =
- PARALLEL_INCLUDES =
-endif
+PARALLEL_AMDIS_SOURCES =
+PARALLEL_INCLUDES =
if USE_PARALLEL_DOMAIN_AMDIS
PARALLEL_AMDIS_SOURCES += \
$(SOURCE_DIR)/parallel/StdMpi.h $(SOURCE_DIR)/parallel/StdMpi.cc \
+ $(SOURCE_DIR)/parallel/ParMetisPartitioner.h $(SOURCE_DIR)/parallel/ParMetisPartitioner.cc \
$(SOURCE_DIR)/parallel/MeshDistributor.h $(SOURCE_DIR)/parallel/MeshDistributor.cc \
$(SOURCE_DIR)/parallel/ParallelDebug.h $(SOURCE_DIR)/parallel/ParallelDebug.cc \
$(SOURCE_DIR)/parallel/ParallelProblemStatBase.h \
@@ -98,6 +86,7 @@ $(SOURCE_DIR)/ProblemInterpolVec.h $(SOURCE_DIR)/ProblemInterpolVec.cc \
$(SOURCE_DIR)/Serializable.h \
$(SOURCE_DIR)/BallProject.h \
$(SOURCE_DIR)/CylinderProject.h \
+$(SOURCE_DIR)/MacroInfo.h $(SOURCE_DIR)/MacroInfo.cc \
$(SOURCE_DIR)/MacroReader.h $(SOURCE_DIR)/MacroReader.cc \
$(SOURCE_DIR)/ValueReader.h $(SOURCE_DIR)/ValueReader.cc \
$(SOURCE_DIR)/Projection.h $(SOURCE_DIR)/Projection.cc \
diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in
index 3c44b066eee8093d6244166fcafee2fc6479b698..3f02633fdd0f983b274e505a6787dc564ce4ab9d 100644
--- a/AMDiS/bin/Makefile.in
+++ b/AMDiS/bin/Makefile.in
@@ -34,9 +34,9 @@ PRE_UNINSTALL = :
POST_UNINSTALL = :
build_triplet = @build@
host_triplet = @host@
-@USE_PARALLEL_AMDIS_TRUE@am__append_1 = -DHAVE_PARALLEL_AMDIS=1
-@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__append_2 = \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__append_1 = \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/StdMpi.h $(SOURCE_DIR)/parallel/StdMpi.cc \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/ParMetisPartitioner.h $(SOURCE_DIR)/parallel/ParMetisPartitioner.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/MeshDistributor.h $(SOURCE_DIR)/parallel/MeshDistributor.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/ParallelDebug.h $(SOURCE_DIR)/parallel/ParallelDebug.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/ParallelProblemStatBase.h \
@@ -44,19 +44,19 @@ host_triplet = @host@
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/MpiHelper.h $(SOURCE_DIR)/parallel/MpiHelper.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ $(SOURCE_DIR)/parallel/ElementObjectData.h $(SOURCE_DIR)/parallel/ElementObjectData.cc
-@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__append_3 = -DHAVE_PARALLEL_DOMAIN_AMDIS=1
-@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__append_4 = -I$(PETSC_DIR)/include -I$(PETSC_DIR)/$(PETSC_ARCH)/include
-@ENABLE_UMFPACK_TRUE@am__append_5 = -DHAVE_UMFPACK=1 -DMTL_HAS_UMFPACK
-@ENABLE_UMFPACK_TRUE@am__append_6 = -I$(LIB_DIR)/UFconfig \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__append_2 = -DHAVE_PARALLEL_DOMAIN_AMDIS=1
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__append_3 = -I$(PETSC_DIR)/include -I$(PETSC_DIR)/$(PETSC_ARCH)/include
+@ENABLE_UMFPACK_TRUE@am__append_4 = -DHAVE_UMFPACK=1 -DMTL_HAS_UMFPACK
+@ENABLE_UMFPACK_TRUE@am__append_5 = -I$(LIB_DIR)/UFconfig \
@ENABLE_UMFPACK_TRUE@ -I$(LIB_DIR)/AMD/Include \
@ENABLE_UMFPACK_TRUE@ -I$(LIB_DIR)/UMFPACK/Include
-@ENABLE_MKL_TRUE@am__append_7 = -DHAVE_MKL=1 -I${MKL_INC}
-@ENABLE_DUNE_TRUE@am__append_8 = -DHAVE_DUNE=1
-@ENABLE_DUNE_TRUE@am__append_9 = -I$(DUNE_DIR)
-@ENABLE_BOOST_TRUE@am__append_10 = -DHAVE_BOOST=1
-@AMDIS_DEBUG_TRUE@am__append_11 = -g -O0 -Wall -DDEBUG=1 $(OPENMP_FLAG) $(INCLUDES) #-pedantic
-@AMDIS_DEBUG_FALSE@am__append_12 = -O2 -Wall -DDEBUG=0 -DNDEBUG $(OPENMP_FLAG) -ftemplate-depth-100 $(INCLUDES) #-pedantic
+@ENABLE_MKL_TRUE@am__append_6 = -DHAVE_MKL=1 -I${MKL_INC}
+@ENABLE_DUNE_TRUE@am__append_7 = -DHAVE_DUNE=1
+@ENABLE_DUNE_TRUE@am__append_8 = -I$(DUNE_DIR)
+@ENABLE_BOOST_TRUE@am__append_9 = -DHAVE_BOOST=1
+@AMDIS_DEBUG_TRUE@am__append_10 = -g -O0 -Wall -DDEBUG=1 $(OPENMP_FLAG) $(INCLUDES) #-pedantic
+@AMDIS_DEBUG_FALSE@am__append_11 = -O2 -Wall -DDEBUG=0 -DNDEBUG $(OPENMP_FLAG) -ftemplate-depth-100 $(INCLUDES) #-pedantic
subdir = bin
DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
@@ -92,6 +92,8 @@ LTLIBRARIES = $(lib_LTLIBRARIES)
libamdis_la_LIBADD =
am__libamdis_la_SOURCES_DIST = $(SOURCE_DIR)/parallel/StdMpi.h \
$(SOURCE_DIR)/parallel/StdMpi.cc \
+ $(SOURCE_DIR)/parallel/ParMetisPartitioner.h \
+ $(SOURCE_DIR)/parallel/ParMetisPartitioner.cc \
$(SOURCE_DIR)/parallel/MeshDistributor.h \
$(SOURCE_DIR)/parallel/MeshDistributor.cc \
$(SOURCE_DIR)/parallel/ParallelDebug.h \
@@ -103,24 +105,12 @@ am__libamdis_la_SOURCES_DIST = $(SOURCE_DIR)/parallel/StdMpi.h \
$(SOURCE_DIR)/parallel/MpiHelper.cc \
$(SOURCE_DIR)/parallel/ElementObjectData.h \
$(SOURCE_DIR)/parallel/ElementObjectData.cc \
- $(PARALLEL_DIR)/ConditionalEstimator.h \
- $(PARALLEL_DIR)/ConditionalEstimator.cc \
- $(PARALLEL_DIR)/ConditionalMarker.h \ \
- $(PARALLEL_DIR)/ConditionalMarker.cc \
- $(PARALLEL_DIR)/ParallelError.h \
- $(PARALLEL_DIR)/ParallelError.hh \
- $(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 $(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)/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 \
@@ -146,7 +136,8 @@ am__libamdis_la_SOURCES_DIST = $(SOURCE_DIR)/parallel/StdMpi.h \
$(SOURCE_DIR)/ProblemInterpolVec.h \
$(SOURCE_DIR)/ProblemInterpolVec.cc \
$(SOURCE_DIR)/Serializable.h $(SOURCE_DIR)/BallProject.h \
- $(SOURCE_DIR)/CylinderProject.h $(SOURCE_DIR)/MacroReader.h \
+ $(SOURCE_DIR)/CylinderProject.h $(SOURCE_DIR)/MacroInfo.h \
+ $(SOURCE_DIR)/MacroInfo.cc $(SOURCE_DIR)/MacroReader.h \
$(SOURCE_DIR)/MacroReader.cc $(SOURCE_DIR)/ValueReader.h \
$(SOURCE_DIR)/ValueReader.cc $(SOURCE_DIR)/Projection.h \
$(SOURCE_DIR)/Projection.cc $(SOURCE_DIR)/SubAssembler.h \
@@ -199,7 +190,7 @@ am__libamdis_la_SOURCES_DIST = $(SOURCE_DIR)/parallel/StdMpi.h \
$(SOURCE_DIR)/Estimator.cc $(SOURCE_DIR)/FiniteElemSpace.h \
$(SOURCE_DIR)/FixVec.h $(SOURCE_DIR)/FixVec.hh \
$(SOURCE_DIR)/FixVecConvert.h $(SOURCE_DIR)/Flag.h \
- $(SOURCE_DIR)/Global.h $(SOURCE_DIR)/UmfPackSolver.h \
+ $(SOURCE_DIR)/Global.h $(SOURCE_DIR)/UmfPackSolver.h \ \
$(SOURCE_DIR)/UmfPackSolver.hh $(SOURCE_DIR)/Lagrange.h \
$(SOURCE_DIR)/Line.h $(SOURCE_DIR)/MacroElement.h \
$(SOURCE_DIR)/MacroWriter.h $(SOURCE_DIR)/Markings.h \
@@ -267,19 +258,13 @@ am__libamdis_la_SOURCES_DIST = $(SOURCE_DIR)/parallel/StdMpi.h \
$(SOURCE_DIR)/time/RosenbrockMethod.h \
$(SOURCE_DIR)/time/RosenbrockMethod.cc
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__objects_1 = libamdis_la-StdMpi.lo \
+@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-ParMetisPartitioner.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-MeshDistributor.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-ParallelDebug.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-PetscSolver.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-MpiHelper.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-ElementObjectData.lo
-@USE_PARALLEL_AMDIS_FALSE@am__objects_2 = $(am__objects_1)
-@USE_PARALLEL_AMDIS_TRUE@am__objects_2 = \
-@USE_PARALLEL_AMDIS_TRUE@ libamdis_la-ConditionalEstimator.lo \
-@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__objects_2 = $(am__objects_1)
am_libamdis_la_OBJECTS = $(am__objects_2) libamdis_la-DOFIndexed.lo \
libamdis_la-GNUPlotWriter.lo libamdis_la-VertexVector.lo \
libamdis_la-PeriodicBC.lo libamdis_la-Recovery.lo \
@@ -292,9 +277,10 @@ am_libamdis_la_OBJECTS = $(am__objects_2) libamdis_la-DOFIndexed.lo \
libamdis_la-ElementData.lo \
libamdis_la-ComponentTraverseInfo.lo libamdis_la-CreatorMap.lo \
libamdis_la-ProblemInterpolScal.lo \
- libamdis_la-ProblemInterpolVec.lo libamdis_la-MacroReader.lo \
- libamdis_la-ValueReader.lo libamdis_la-Projection.lo \
- libamdis_la-SubAssembler.lo libamdis_la-ZeroOrderAssembler.lo \
+ libamdis_la-ProblemInterpolVec.lo libamdis_la-MacroInfo.lo \
+ libamdis_la-MacroReader.lo libamdis_la-ValueReader.lo \
+ libamdis_la-Projection.lo libamdis_la-SubAssembler.lo \
+ libamdis_la-ZeroOrderAssembler.lo \
libamdis_la-FirstOrderAssembler.lo \
libamdis_la-SecondOrderAssembler.lo libamdis_la-Assembler.lo \
libamdis_la-AdaptInfo.lo libamdis_la-Marker.lo \
@@ -511,28 +497,13 @@ SOURCE_DIR = ../src
LIB_DIR = ../lib
PARALLEL_DIR = $(SOURCE_DIR)
PARMETIS_DIR = ../lib/ParMetis-3.1
-AMDIS_INCLUDES = -I$(SOURCE_DIR) $(am__append_4) $(am__append_6) \
- $(am__append_9)
-libamdis_la_CXXFLAGS = $(am__append_1) $(am__append_3) $(am__append_5) \
- $(am__append_7) $(am__append_8) $(am__append_10) \
- $(am__append_11) $(am__append_12)
-@USE_PARALLEL_AMDIS_FALSE@PARALLEL_AMDIS_SOURCES = $(am__append_2)
-@USE_PARALLEL_AMDIS_TRUE@PARALLEL_AMDIS_SOURCES = $(PARALLEL_DIR)/ConditionalEstimator.h \
-@USE_PARALLEL_AMDIS_TRUE@ $(PARALLEL_DIR)/ConditionalEstimator.cc \
-@USE_PARALLEL_AMDIS_TRUE@ $(PARALLEL_DIR)/ConditionalMarker.h \ \
-@USE_PARALLEL_AMDIS_TRUE@ $(PARALLEL_DIR)/ConditionalMarker.cc \
-@USE_PARALLEL_AMDIS_TRUE@ $(PARALLEL_DIR)/ParallelError.h \
-@USE_PARALLEL_AMDIS_TRUE@ $(PARALLEL_DIR)/ParallelError.hh \
-@USE_PARALLEL_AMDIS_TRUE@ $(PARALLEL_DIR)/ParallelProblem.h \
-@USE_PARALLEL_AMDIS_TRUE@ $(PARALLEL_DIR)/ParallelProblem.cc \
-@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_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)
+AMDIS_INCLUDES = -I$(SOURCE_DIR) $(am__append_3) $(am__append_5) \
+ $(am__append_8)
+libamdis_la_CXXFLAGS = $(am__append_2) $(am__append_4) $(am__append_6) \
+ $(am__append_7) $(am__append_9) $(am__append_10) \
+ $(am__append_11)
+PARALLEL_AMDIS_SOURCES = $(am__append_1)
+PARALLEL_INCLUDES =
TEMPLATE_INCLUDES = -I../lib/mtl4
INCLUDES = $(AMDIS_INCLUDES) $(PARALLEL_INCLUDES) $(TEMPLATE_INCLUDES)
libamdis_la_SOURCES = \
@@ -563,6 +534,7 @@ $(SOURCE_DIR)/ProblemInterpolVec.h $(SOURCE_DIR)/ProblemInterpolVec.cc \
$(SOURCE_DIR)/Serializable.h \
$(SOURCE_DIR)/BallProject.h \
$(SOURCE_DIR)/CylinderProject.h \
+$(SOURCE_DIR)/MacroInfo.h $(SOURCE_DIR)/MacroInfo.cc \
$(SOURCE_DIR)/MacroReader.h $(SOURCE_DIR)/MacroReader.cc \
$(SOURCE_DIR)/ValueReader.h $(SOURCE_DIR)/ValueReader.cc \
$(SOURCE_DIR)/Projection.h $(SOURCE_DIR)/Projection.cc \
@@ -810,8 +782,6 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-CoarseningManager2d.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-CoarseningManager3d.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ComponentTraverseInfo.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ConditionalEstimator.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ConditionalMarker.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-CreatorMap.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-DOFAdmin.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-DOFIndexed.Plo@am__quote@
@@ -844,6 +814,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-LeafData.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Line.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-MacroElement.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-MacroInfo.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-MacroReader.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-MacroWriter.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Marker.Plo@am__quote@
@@ -856,13 +827,11 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-OperatorTerm.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ParMetisPartitioner.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ParallelDebug.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ParallelProblem.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Parameters.Plo@am__quote@
@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-PetscSolver.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@
@@ -943,6 +912,13 @@ libamdis_la-StdMpi.lo: $(SOURCE_DIR)/parallel/StdMpi.cc
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-StdMpi.lo `test -f '$(SOURCE_DIR)/parallel/StdMpi.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/StdMpi.cc
+libamdis_la-ParMetisPartitioner.lo: $(SOURCE_DIR)/parallel/ParMetisPartitioner.cc
+@am__fastdepCXX_TRUE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-ParMetisPartitioner.lo -MD -MP -MF $(DEPDIR)/libamdis_la-ParMetisPartitioner.Tpo -c -o libamdis_la-ParMetisPartitioner.lo `test -f '$(SOURCE_DIR)/parallel/ParMetisPartitioner.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/ParMetisPartitioner.cc
+@am__fastdepCXX_TRUE@ $(am__mv) $(DEPDIR)/libamdis_la-ParMetisPartitioner.Tpo $(DEPDIR)/libamdis_la-ParMetisPartitioner.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/parallel/ParMetisPartitioner.cc' object='libamdis_la-ParMetisPartitioner.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ParMetisPartitioner.lo `test -f '$(SOURCE_DIR)/parallel/ParMetisPartitioner.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/ParMetisPartitioner.cc
+
libamdis_la-MeshDistributor.lo: $(SOURCE_DIR)/parallel/MeshDistributor.cc
@am__fastdepCXX_TRUE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-MeshDistributor.lo -MD -MP -MF $(DEPDIR)/libamdis_la-MeshDistributor.Tpo -c -o libamdis_la-MeshDistributor.lo `test -f '$(SOURCE_DIR)/parallel/MeshDistributor.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/MeshDistributor.cc
@am__fastdepCXX_TRUE@ $(am__mv) $(DEPDIR)/libamdis_la-MeshDistributor.Tpo $(DEPDIR)/libamdis_la-MeshDistributor.Plo
@@ -978,41 +954,6 @@ libamdis_la-ElementObjectData.lo: $(SOURCE_DIR)/parallel/ElementObjectData.cc
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ElementObjectData.lo `test -f '$(SOURCE_DIR)/parallel/ElementObjectData.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/parallel/ElementObjectData.cc
-libamdis_la-ConditionalEstimator.lo: $(PARALLEL_DIR)/ConditionalEstimator.cc
-@am__fastdepCXX_TRUE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-ConditionalEstimator.lo -MD -MP -MF $(DEPDIR)/libamdis_la-ConditionalEstimator.Tpo -c -o libamdis_la-ConditionalEstimator.lo `test -f '$(PARALLEL_DIR)/ConditionalEstimator.cc' || echo '$(srcdir)/'`$(PARALLEL_DIR)/ConditionalEstimator.cc
-@am__fastdepCXX_TRUE@ $(am__mv) $(DEPDIR)/libamdis_la-ConditionalEstimator.Tpo $(DEPDIR)/libamdis_la-ConditionalEstimator.Plo
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(PARALLEL_DIR)/ConditionalEstimator.cc' object='libamdis_la-ConditionalEstimator.lo' libtool=yes @AMDEPBACKSLASH@
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ConditionalEstimator.lo `test -f '$(PARALLEL_DIR)/ConditionalEstimator.cc' || echo '$(srcdir)/'`$(PARALLEL_DIR)/ConditionalEstimator.cc
-
-libamdis_la-ConditionalMarker.lo: $(PARALLEL_DIR)/ConditionalMarker.cc
-@am__fastdepCXX_TRUE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-ConditionalMarker.lo -MD -MP -MF $(DEPDIR)/libamdis_la-ConditionalMarker.Tpo -c -o libamdis_la-ConditionalMarker.lo `test -f '$(PARALLEL_DIR)/ConditionalMarker.cc' || echo '$(srcdir)/'`$(PARALLEL_DIR)/ConditionalMarker.cc
-@am__fastdepCXX_TRUE@ $(am__mv) $(DEPDIR)/libamdis_la-ConditionalMarker.Tpo $(DEPDIR)/libamdis_la-ConditionalMarker.Plo
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(PARALLEL_DIR)/ConditionalMarker.cc' object='libamdis_la-ConditionalMarker.lo' libtool=yes @AMDEPBACKSLASH@
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ConditionalMarker.lo `test -f '$(PARALLEL_DIR)/ConditionalMarker.cc' || echo '$(srcdir)/'`$(PARALLEL_DIR)/ConditionalMarker.cc
-
-libamdis_la-ParallelProblem.lo: $(PARALLEL_DIR)/ParallelProblem.cc
-@am__fastdepCXX_TRUE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-ParallelProblem.lo -MD -MP -MF $(DEPDIR)/libamdis_la-ParallelProblem.Tpo -c -o libamdis_la-ParallelProblem.lo `test -f '$(PARALLEL_DIR)/ParallelProblem.cc' || echo '$(srcdir)/'`$(PARALLEL_DIR)/ParallelProblem.cc
-@am__fastdepCXX_TRUE@ $(am__mv) $(DEPDIR)/libamdis_la-ParallelProblem.Tpo $(DEPDIR)/libamdis_la-ParallelProblem.Plo
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(PARALLEL_DIR)/ParallelProblem.cc' object='libamdis_la-ParallelProblem.lo' libtool=yes @AMDEPBACKSLASH@
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ParallelProblem.lo `test -f '$(PARALLEL_DIR)/ParallelProblem.cc' || echo '$(srcdir)/'`$(PARALLEL_DIR)/ParallelProblem.cc
-
-libamdis_la-ParMetisPartitioner.lo: $(PARALLEL_DIR)/ParMetisPartitioner.cc
-@am__fastdepCXX_TRUE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-ParMetisPartitioner.lo -MD -MP -MF $(DEPDIR)/libamdis_la-ParMetisPartitioner.Tpo -c -o libamdis_la-ParMetisPartitioner.lo `test -f '$(PARALLEL_DIR)/ParMetisPartitioner.cc' || echo '$(srcdir)/'`$(PARALLEL_DIR)/ParMetisPartitioner.cc
-@am__fastdepCXX_TRUE@ $(am__mv) $(DEPDIR)/libamdis_la-ParMetisPartitioner.Tpo $(DEPDIR)/libamdis_la-ParMetisPartitioner.Plo
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(PARALLEL_DIR)/ParMetisPartitioner.cc' object='libamdis_la-ParMetisPartitioner.lo' libtool=yes @AMDEPBACKSLASH@
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --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@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --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@ $(am__mv) $(DEPDIR)/libamdis_la-PollutionError.Tpo $(DEPDIR)/libamdis_la-PollutionError.Plo
-@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 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --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@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --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@ $(am__mv) $(DEPDIR)/libamdis_la-DOFIndexed.Tpo $(DEPDIR)/libamdis_la-DOFIndexed.Plo
@@ -1146,6 +1087,13 @@ libamdis_la-ProblemInterpolVec.lo: $(SOURCE_DIR)/ProblemInterpolVec.cc
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ProblemInterpolVec.lo `test -f '$(SOURCE_DIR)/ProblemInterpolVec.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/ProblemInterpolVec.cc
+libamdis_la-MacroInfo.lo: $(SOURCE_DIR)/MacroInfo.cc
+@am__fastdepCXX_TRUE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-MacroInfo.lo -MD -MP -MF $(DEPDIR)/libamdis_la-MacroInfo.Tpo -c -o libamdis_la-MacroInfo.lo `test -f '$(SOURCE_DIR)/MacroInfo.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/MacroInfo.cc
+@am__fastdepCXX_TRUE@ $(am__mv) $(DEPDIR)/libamdis_la-MacroInfo.Tpo $(DEPDIR)/libamdis_la-MacroInfo.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/MacroInfo.cc' object='libamdis_la-MacroInfo.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-MacroInfo.lo `test -f '$(SOURCE_DIR)/MacroInfo.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/MacroInfo.cc
+
libamdis_la-MacroReader.lo: $(SOURCE_DIR)/MacroReader.cc
@am__fastdepCXX_TRUE@ $(LIBTOOL) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-MacroReader.lo -MD -MP -MF $(DEPDIR)/libamdis_la-MacroReader.Tpo -c -o libamdis_la-MacroReader.lo `test -f '$(SOURCE_DIR)/MacroReader.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/MacroReader.cc
@am__fastdepCXX_TRUE@ $(am__mv) $(DEPDIR)/libamdis_la-MacroReader.Tpo $(DEPDIR)/libamdis_la-MacroReader.Plo
diff --git a/AMDiS/src/CoarseningManager1d.cc b/AMDiS/src/CoarseningManager1d.cc
index 719f682d83f587efe623037b3cf5b0803111d290..9b19b9c5e2b0ce05057e9456c03861949e0456d2 100644
--- a/AMDiS/src/CoarseningManager1d.cc
+++ b/AMDiS/src/CoarseningManager1d.cc
@@ -75,7 +75,7 @@ namespace AMDiS {
}
if (mesh->getNumberOfDOFs(CENTER) && !mesh->queryCoarseDOFs()) {
- parent->setDOF(mesh->getNode(CENTER), mesh->getDOF(CENTER));
+ parent->setDOF(mesh->getNode(CENTER), mesh->getDof(CENTER));
}
/*--------------------------------------------------------------------------*/
diff --git a/AMDiS/src/CoarseningManager2d.cc b/AMDiS/src/CoarseningManager2d.cc
index 1a4d6b3bc43d454dc95a4b4079fc518087dcdab8..24454798df6c91bbddfb63a438623cf9ac01f6d6 100644
--- a/AMDiS/src/CoarseningManager2d.cc
+++ b/AMDiS/src/CoarseningManager2d.cc
@@ -84,7 +84,7 @@ namespace AMDiS {
// get new dof on el at the midpoint of the coarsening edge
if (!el->getDOF(node + 2)) {
- el->setDOF(node + 2, mesh->getDOF(EDGE));
+ el->setDOF(node + 2, mesh->getDof(EDGE));
if (neigh)
neigh->setDOF(node + 2, const_cast<int*>(el->getDOF(node + 2)));
}
diff --git a/AMDiS/src/CoarseningManager3d.cc b/AMDiS/src/CoarseningManager3d.cc
index df4cb88bac45d1a38b16b345fc7900ac83d518ef..8115f126cd139b9ef2a3b39714f9a84246b55194 100644
--- a/AMDiS/src/CoarseningManager3d.cc
+++ b/AMDiS/src/CoarseningManager3d.cc
@@ -324,7 +324,7 @@ namespace AMDiS {
/****************************************************************************/
node = mesh->getNode(EDGE);
if (!(dof = const_cast<int*>( el->getDOF(node))))
- dof = mesh->getDOF(EDGE);
+ dof = mesh->getDof(EDGE);
} else {
dof = NULL;
}
diff --git a/AMDiS/src/ConditionalEstimator.cc b/AMDiS/src/ConditionalEstimator.cc
deleted file mode 100644
index 3ab73099913414a8fb61fa85a1f39026b059d24a..0000000000000000000000000000000000000000
--- a/AMDiS/src/ConditionalEstimator.cc
+++ /dev/null
@@ -1,120 +0,0 @@
-#include "ConditionalEstimator.h"
-#include "ElInfo.h"
-#include "PartitionElementData.h"
-#include "Element.h"
-#include "Traverse.h"
-#include "mpi.h"
-
-namespace AMDiS {
-
- ConditionalEstimator::ConditionalEstimator(Estimator* decorated)
- : decoratedEstimator_(decorated),
- elementCount_(0),
- row_(decorated ? decorated->getRow() : 0)
- {
- FUNCNAME("ConditionalEstimator::ConditionalEstimator(Estimator* decorated)");
-
- if (decorated!=NULL) {
- name = decorated->getName();
- int swap;
- int retVal =
- Parameters::getGlobalParameter(0, name + "->estimate outer", "%d", &swap);
-
- if (retVal > 0)
- estimateOut = (bool)swap;
- else
- estimateOut = false;
- }
- }
-
- double ConditionalEstimator::estimate(double ts)
- {
- if (decoratedEstimator_) {
-
- double partition_sum = 0.0;
- elementCount_ = 0;
- decoratedEstimator_->init(ts);
-
- mesh = decoratedEstimator_->getMesh();
- traverseFlag = decoratedEstimator_->getTraverseFlag();
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
- while (elInfo) {
- PartitionElementData *elData = dynamic_cast<PartitionElementData*>
- (elInfo->getElement()->getElementData(PARTITION_ED));
-
- TEST_EXIT_DBG(elInfo->getElement()->isLeaf())("element not leaf\n");
- TEST_EXIT_DBG(elData)("no partition data on leaf element %d (rank %d)\n",
- elInfo->getElement()->getIndex(),
- MPI::COMM_WORLD.Get_rank());
-
- PartitionStatus status = elData->getPartitionStatus();
-
- if (status == IN || status == OVERLAP || (estimateOut && status==OUT))
- decoratedEstimator_->estimateElement(elInfo);
- else
- elInfo->getElement()->setMark(0);
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
- while (elInfo) {
- PartitionElementData *elData = dynamic_cast<PartitionElementData*>
- (elInfo->getElement()->getElementData(PARTITION_ED));
- PartitionStatus status = elData->getPartitionStatus();
- if (status == IN || status==OVERLAP || (estimateOut && status==OUT)) {
- elementCount_++;
- partition_sum += elInfo->getElement()->getEstimation(row_);
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- // !!! test !!!
- double total_sum = 0.0;
- MPI::COMM_WORLD.Allreduce(&partition_sum,
- &total_sum,
- 1,
- MPI_DOUBLE,
- MPI_SUM);
-
- decoratedEstimator_->setErrorSum(total_sum);
-
- double total_max = 0.0;
- est_max = decoratedEstimator_->getErrorMax();
- MPI::COMM_WORLD.Allreduce(&est_max,
- &total_max,
- 1,
- MPI_DOUBLE,
- MPI_MAX);
-
- decoratedEstimator_->setErrorMax(total_max);
-
- decoratedEstimator_->exit();
-
- est_sum = decoratedEstimator_->getErrorSum();
- est_max = decoratedEstimator_->getErrorMax();
-
- // !!! test !!!
-#if 0
- decoratedEstimator_->exit();
-
- est_sum = sqrt(partition_sum); //decoratedEstimator_->getErrorSum();
- est_t_sum = decoratedEstimator_->getTimeEst();
- est_max = decoratedEstimator_->getErrorMax();
-#endif
-
- // MSG("rank %d , estimate %e (total %e) elements %d (total %d)\n",
- // MPI::COMM_WORLD.Get_rank(), est_sum, total_sum,
- // elementCount_,
- // mesh->getNumberOfLeaves());
-
- return est_sum;
- } else {
- return 0.0;
- }
- }
-
-}
diff --git a/AMDiS/src/ConditionalEstimator.h b/AMDiS/src/ConditionalEstimator.h
deleted file mode 100644
index 94940d7fc5cce1b4b19bb322edadc6b18d45c184..0000000000000000000000000000000000000000
--- a/AMDiS/src/ConditionalEstimator.h
+++ /dev/null
@@ -1,72 +0,0 @@
-// ============================================================================
-// == ==
-// == AMDiS - Adaptive multidimensional simulations ==
-// == ==
-// ============================================================================
-// == ==
-// == crystal growth group ==
-// == ==
-// == Stiftung caesar ==
-// == Ludwig-Erhard-Allee 2 ==
-// == 53175 Bonn ==
-// == germany ==
-// == ==
-// ============================================================================
-// == ==
-// == http://www.caesar.de/cg/AMDiS ==
-// == ==
-// ============================================================================
-
-/** \file ConditionalEstimator.h */
-
-/** \defgroup Estimator Estimator module
- * @{ <img src="estimator.png"> @}
- */
-
-#ifndef AMDIS_CONDITIONALESTIMATOR_H
-#define AMDIS_CONDITIONALESTIMATOR_H
-
-#include "Estimator.h"
-
-namespace AMDiS {
-
- /**
- * \ingroup Estimator
- *
- * \brief
- * Estimator for scalar problems.
- */
- class ConditionalEstimator : public Estimator
- {
- public:
- /**
- * Reads the estimateOuter parameter. default value is false!
- * (i.e.: the ouer part is not estimated)
- */
- ConditionalEstimator(Estimator *decorated);
-
- double estimate(double timestep);
-
- inline int getElementCount()
- {
- return elementCount_;
- }
-
- inline Estimator *getDecoratedEstimator()
- {
- return decoratedEstimator_;
- }
-
- protected:
- Estimator *decoratedEstimator_;
-
- int elementCount_;
-
- int row_;
-
- bool estimateOut;
- };
-
-}
-
-#endif // AMDIS_ESTIMATOR_H
diff --git a/AMDiS/src/ConditionalMarker.cc b/AMDiS/src/ConditionalMarker.cc
deleted file mode 100644
index 9aeeb1493460bfe63c5b8cdb242d5e04fab72b91..0000000000000000000000000000000000000000
--- a/AMDiS/src/ConditionalMarker.cc
+++ /dev/null
@@ -1,85 +0,0 @@
-#include "ConditionalMarker.h"
-
-namespace AMDiS
-{
- ConditionalMarker::ConditionalMarker(const std::string name_,
- int row,
- Marker *decoratedMarker,
- int globalCoarseGridLevel,
- int localCoarseGridLevel)
- : Marker(name_, row),
- decoratedMarker_(decoratedMarker),
- globalCoarseGridLevel_(globalCoarseGridLevel),
- localCoarseGridLevel_(localCoarseGridLevel)
- {
- if (decoratedMarker != NULL)
- name = decoratedMarker->getName();
- else
- std::cout << "no decorated Marker, mark outer path is" << name << "->mark outer" << std::endl;
-
- int swap;
- int retVal = Parameters::getGlobalParameter(0, name + "->mark outer", "%d", &swap);
-
- if (retVal>0)
- markOuter = (bool)swap;
- else
- markOuter = false;
-
- MSG("markOuter:" + markOuter);
- }
-
- void ConditionalMarker::initMarking(AdaptInfo *adaptInfo, Mesh *mesh)
- {
- if (decoratedMarker_)
- decoratedMarker_->initMarking(adaptInfo, mesh);
- }
-
- void ConditionalMarker::finishMarking(AdaptInfo *adaptInfo)
- {
- if (decoratedMarker_) {
- int tmp = decoratedMarker_->getElMarkRefine();
- MPI::COMM_WORLD.Allreduce(&tmp,
- &elMarkRefine,
- 1,
- MPI_INT,
- MPI_SUM);
- tmp = decoratedMarker_->getElMarkCoarsen();
- MPI::COMM_WORLD.Allreduce(&tmp,
- &elMarkCoarsen,
- 1,
- MPI_INT,
- MPI_SUM);
- decoratedMarker_->finishMarking(adaptInfo);
- }
- }
-
-
- void ConditionalMarker::markElement(AdaptInfo *adaptInfo, ElInfo *elInfo)
- {
- FUNCNAME("ConditionalMarker::markElement()");
-
- if (decoratedMarker_) {
- PartitionElementData *elData =
- dynamic_cast<PartitionElementData*>(elInfo->getElement()->
- getElementData(PARTITION_ED));
-
- TEST_EXIT_DBG(elData)("no partition data\n");
-
- decoratedMarker_->markElement(adaptInfo, elInfo);
-
- if (!markOuter && elData->getPartitionStatus() == OUT) {
- Element *element = elInfo->getElement();
- if (element->getMark() > 0)
- element->setMark(0);
- }
-
- int minLevel =
- elData->getPartitionStatus() != OUT ?
- localCoarseGridLevel_ :
- globalCoarseGridLevel_;
-
- if (elData->getLevel() + elInfo->getElement()->getMark() < minLevel)
- elInfo->getElement()->setMark(-(elData->getLevel() - minLevel));
- }
- }
-}
diff --git a/AMDiS/src/ConditionalMarker.h b/AMDiS/src/ConditionalMarker.h
deleted file mode 100644
index 8328d6572966288f26da052718a1058b3d0047d9..0000000000000000000000000000000000000000
--- a/AMDiS/src/ConditionalMarker.h
+++ /dev/null
@@ -1,67 +0,0 @@
-// ============================================================================
-// == ==
-// == AMDiS - Adaptive multidimensional simulations ==
-// == ==
-// ============================================================================
-// == ==
-// == crystal growth group ==
-// == ==
-// == Stiftung caesar ==
-// == Ludwig-Erhard-Allee 2 ==
-// == 53175 Bonn ==
-// == germany ==
-// == ==
-// ============================================================================
-// == ==
-// == http://www.caesar.de/cg/AMDiS ==
-// == ==
-// ============================================================================
-
-/** \file ConditionalMarker.h */
-
-#ifndef AMDIS_CONDITIONALMARKER_H
-#define AMDIS_CONDITIONALMARKER_H
-
-#include "mpi.h"
-#include "Marker.h"
-#include "Marker.h"
-#include "PartitionElementData.h"
-#include "Traverse.h"
-#include "ElInfo.h"
-
-namespace AMDiS {
-
- /**
- * \ingroup Adaption
- *
- * \brief
- *
- */
- class ConditionalMarker : public Marker
- {
- public:
- /// Constructor.
- ConditionalMarker(const std::string name,
- int row,
- Marker *decoratedMarker,
- int globalCoarseGridLevel,
- int localCoarseGridLevel);
-
-
- /// Implementation of MarkScal::initMarking().
- virtual void initMarking(AdaptInfo *adaptInfo, Mesh *mesh);
-
- virtual void finishMarking(AdaptInfo *adaptInfo);
-
- void markElement(AdaptInfo *adaptInfo, ElInfo *elInfo);
-
- protected:
- Marker *decoratedMarker_;
- int globalCoarseGridLevel_;
- int localCoarseGridLevel_;
- bool markOuter;
- };
-
-}
-
-#endif
diff --git a/AMDiS/src/FixVecConvert.h b/AMDiS/src/FixVecConvert.h
index ad9d710e0ffa46f023cb3f94be8b77c30c5dbd22..4fddd19be1767162b6f32fc10e640e2095a7a832 100644
--- a/AMDiS/src/FixVecConvert.h
+++ b/AMDiS/src/FixVecConvert.h
@@ -32,7 +32,7 @@ namespace AMDiS {
public:
static FixVec<T,d1>& convertVec(FixVec<T,d2>& rhs, Mesh* mesh)
{
- TEST_EXIT(mesh->getGeo(d1)==mesh->getGeo(d2))("Incompatible dimensions.\n");
+ TEST_EXIT(mesh->getGeo(d1) == mesh->getGeo(d2))("Incompatible dimensions!\n");
return reinterpret_cast<FixVec<T,d1>&>(rhs);
}
};
diff --git a/AMDiS/src/MacroInfo.cc b/AMDiS/src/MacroInfo.cc
new file mode 100644
index 0000000000000000000000000000000000000000..486d3c366713660668b42273a01bff7e70da9489
--- /dev/null
+++ b/AMDiS/src/MacroInfo.cc
@@ -0,0 +1,575 @@
+#include "MacroInfo.h"
+#include "Mesh.h"
+#include "MacroReader.h"
+#include "FixVecConvert.h"
+#include "SurfaceRegion_ED.h"
+#include "ElementRegion_ED.h"
+
+namespace AMDiS {
+
+ void MacroInfo::fill(Mesh *pmesh, int nElements, int nVertices)
+ {
+ FUNCNAME("MacroInfo::fill()");
+
+ TEST_EXIT(pmesh)("no mesh\n");
+
+ int dim = pmesh->getDim();
+ mesh = pmesh;
+
+ mesh->setNumberOfElements(nElements);
+ mesh->setNumberOfLeaves(nElements);
+ mesh->setNumberOfVertices(nVertices);
+
+ for (int i = 0; i < nElements; i++) {
+ MacroElement *newMacro = new MacroElement(mesh->getDim());
+ mel.push_back(newMacro);
+ mesh->addMacroElement(mel[i]);
+ }
+
+ dof = new DegreeOfFreedom*[nVertices];
+ coords = new WorldVector<double>[nVertices];
+ mel_vertex = new int*[nElements];
+
+ for (int i = 0; i < nElements; i++)
+ mel_vertex[i] = new int[mesh->getGeo(VERTEX)];
+
+ for (int i = 0; i < nVertices; i++)
+ dof[i] = mesh->getDof(VERTEX);
+
+ for (int i = 0; i < nElements; i++) {
+ mel[i]->element = mesh->createNewElement();
+ mel[i]->index = i;
+
+ if (dim == 3)
+ mel[i]->elType = 0;
+ }
+ neigh_set = false;
+ bound_set = false;
+ }
+
+
+ void MacroInfo::clear()
+ {
+ for (int i = 0; i < mesh->getNumberOfMacros(); i++)
+ delete [] mel_vertex[i];
+
+ delete [] mel_vertex;
+ delete [] coords;
+ coords = NULL;
+ delete [] dof;
+ dof = NULL;
+
+ mesh = NULL;
+ neigh_set = false;
+ }
+
+
+ /****************************************************************************/
+ /* read_indices() reads dim+1 indices from file into id[0-dim], */
+ /* returns true if dim+1 inputs arguments could be read successfully by */
+ /* fscanf(), else false */
+ /****************************************************************************/
+
+ int MacroInfo::read_indices(FILE *file, DimVec<int> &id)
+ {
+ int dim = mesh->getDim();
+
+ for (int i = 0; i <= dim; i++)
+ if (fscanf(file, "%d", &id[i]) != 1)
+ return(false);
+
+ return(true);
+ }
+
+#define N_KEYS 14
+#define N_MIN_KEYS 7
+ static const char *keys[N_KEYS] = {
+ "DIM", // 0
+ "DIM_OF_WORLD", // 1
+ "number of vertices", // 2
+ "number of elements", // 3
+ "vertex coordinates", // 4
+ "element vertices", // 5
+ "element boundaries", // 6
+ "element neighbours", // 7
+ "element type", // 8
+ "projections", // 9
+ "element region", // 10
+ "surface region", // 11
+ "mesh name", // 12
+ "time" // 13
+ };
+
+
+ static int get_key_no(const char *key)
+ {
+ for (int i = 0; i < N_KEYS; i++)
+ if (!strcmp(keys[i], key))
+ return i;
+
+ return -1;
+ }
+
+#include <ctype.h>
+
+ static const char *read_key(const char *line)
+ {
+ static char key[100];
+ char *k = key;
+
+ while (isspace(*line))
+ line++;
+ while ((*k++ = *line++) != ':');
+ *--k = '\0';
+
+ return const_cast<const char *>(key);
+ }
+
+
+ void MacroInfo::readAMDiSMacro(std::string filename, Mesh* mesh)
+ {
+ FUNCNAME("MacroInfo::readAMDiSMacro()");
+
+ int dim, dow;
+ int nElements, nVertices;
+ int j, k;
+ double dbl;
+ char line[256];
+ int line_no, n_keys, sort_key[N_KEYS], nv_key, ne_key;
+ int key_def[N_KEYS] = {0,0,0,0,0,0,0,0,0,0,0,0};
+ const char *key;
+ DimVec<int> *ind = NULL;
+
+ TEST_EXIT(filename != "")("No filename specified!\n");
+ TEST_EXIT(filename.length() < static_cast<unsigned int>(127))
+ ("can only handle filenames up to 127 characters\n");
+
+ FILE *file = fopen(filename.c_str(), "r");
+ TEST_EXIT(file)("cannot open file %s\n", filename.c_str());
+
+
+ // === Looking for all keys in the macro file. ===
+
+ line_no = n_keys = 0;
+ while (fgets(line, 255, file)) {
+ line_no++;
+ if (!strchr(line, ':'))
+ continue;
+ key = read_key(line);
+ int i_key = get_key_no(key);
+ TEST_EXIT(i_key >= 0)
+ ("macro file %s must not contain key %s on line %d\n",
+ filename.c_str(), key, line_no);
+ TEST_EXIT(!key_def[i_key])("key %s defined second time on line %d in file %s\n");
+
+ sort_key[n_keys++] = i_key;
+ key_def[i_key] = true;
+ }
+ fclose(file);
+
+
+ // === Test, if there is data for every key and if all is defined in ===
+ // === right order. ===
+
+ for (int i_key = 0; i_key < N_MIN_KEYS; i_key++) {
+ for (j = 0; j < n_keys; j++)
+ if (sort_key[j] == i_key)
+ break;
+
+ TEST_EXIT(j < n_keys)("You do not have specified data for %s in %s\n",
+ keys[i_key], filename.c_str());
+
+ for (j = 0; j < n_keys; j++)
+ if (sort_key[j] == 2) break;
+ nv_key = j;
+ for (j = 0; j < n_keys; j++)
+ if (sort_key[j] == 3) break;
+ ne_key = j;
+
+ switch (i_key) {
+ case 0:
+ case 1:
+ TEST_EXIT(sort_key[i_key] < 2)
+ ("You have to specify DIM or mesh->getGeo(WORLD) before all other data\n");
+ break;
+ case 4:
+ TEST_EXIT(nv_key < i_key)
+ ("Before reading data for %s, you have to specify the %s in file\n",
+ keys[4], keys[2], filename.c_str());
+ break;
+ case 5:
+ TEST_EXIT(nv_key < i_key && ne_key < i_key)
+ ("Before reading data for %s, you have to specify the %s and %s in file %s\n",
+ keys[5], keys[3], keys[2], filename.c_str());
+ case 6:
+ case 7:
+ case 8:
+ TEST_EXIT(ne_key < i_key)
+ ("Before reading data for %s, you have to specify the %s in file %s\n",
+ keys[i_key], keys[3], filename.c_str());
+ }
+ }
+
+ for (int i_key = 0; i_key < N_KEYS; i_key++)
+ key_def[i_key] = false;
+
+
+ // === And now, reading data. ===
+
+ file = fopen(filename.c_str(), "r");
+ TEST_EXIT(file)("cannot open file %s\n", filename.c_str());
+
+ int result;
+
+ for (int i_key = 0; i_key < n_keys; i_key++) {
+
+ switch (sort_key[i_key]) {
+
+ case 0:
+ // line "DIM"
+ result = fscanf(file, "%*s %d", &dim);
+ TEST_EXIT(result == 1)("cannot read DIM correctly in file %s\n", filename.c_str());
+
+ ind = new DimVec<int>(dim, NO_INIT);
+
+ key_def[0] = true;
+ break;
+
+ case 1:
+ // line "DIM_OF_WORLD"
+ result = fscanf(file, "%*s %d", &dow);
+ TEST_EXIT(result == 1)
+ ("cannot read Global::getGeo(WORLD) correctly in file %s\n", filename.c_str());
+ TEST_EXIT(dow == Global::getGeo(WORLD))
+ ("dimension of world = %d != Global::getGeo(WORLD) = %d\n",
+ dow, Global::getGeo(WORLD));
+
+ key_def[1] = true;
+ break;
+
+ case 2:
+ // line "number of vertices"
+ result = fscanf(file, "%*s %*s %*s %d", &nVertices);
+ TEST_EXIT(result == 1)
+ ("cannot read number of vertices correctly in file %s\n", filename.c_str());
+ TEST_EXIT(nVertices > 0)
+ ("number of vertices = %d must be bigger than 0\n", nVertices);
+
+ key_def[2] = true;
+ if (key_def[3])
+ fill(mesh, nElements, nVertices);
+ break;
+
+ case 3:
+ // line "number of elements"
+ result = fscanf(file, "%*s %*s %*s %d", &nElements);
+ TEST_EXIT(result == 1)
+ ("cannot read number of elements correctly in file %s\n", filename.c_str());
+ TEST_EXIT(nElements > 0)
+ ("number of elements = %d must be bigger than 0\n", nElements);
+
+ key_def[3] = true;
+ if (key_def[2])
+ fill(mesh, nElements, nVertices);
+ break;
+
+ case 4:
+ // block "vertex coordinates"
+ fscanf(file, "%*s %*s");
+ for (int i = 0; i < nVertices; i++) {
+ for (j = 0; j <Global::getGeo(WORLD) ; j++) {
+ result = fscanf(file, "%lf", &dbl);
+ TEST_EXIT(result == 1)
+ ("error while reading coordinates, check file %s\n", filename.c_str());
+ coords[i][j] = dbl;
+ }
+ }
+ key_def[4] = true;
+ break;
+
+ case 5:
+ // block "element vertices"
+ fscanf(file, "%*s %*s");
+
+ // === Global number of vertices on a single element. ===
+
+ for (int i = 0; i < nElements; i++) {
+ result = read_indices(file, *ind);
+ TEST_EXIT(result)
+ ("cannot read vertex indices of element %d in file %s\n", i, filename.c_str());
+
+ for (k = 0; k < mesh->getGeo(VERTEX); k++)
+ mel_vertex[i][k] = (*ind)[k];
+ }
+
+ key_def[5] = true;
+ break;
+
+ case 6:
+ // block "element boundaries"
+ fscanf(file, "%*s %*s");
+
+
+ // === MEL boundary pointers. ===
+
+ for (int i = 0; i < nElements; i++) {
+ // boundary information of ith element
+
+ result = read_indices(file, *ind);
+ TEST_EXIT(result)
+ ("cannot read boundary type of element %d in file %s\n", i, filename.c_str());
+
+ // fill boundary of macro-element
+ MacroReader::fillMelBoundary(mesh,
+ mel[i],
+ VecConv<int,NEIGH,PARTS>::convertVec((*ind), mesh));
+ }
+
+ this->fillBoundaryInfo(mesh);
+
+ bound_set = true;
+ key_def[6] = true;
+ break;
+
+ case 7:
+ // block "element neighbours"
+ fscanf(file, "%*s %*s");
+
+ // === Fill MEL neighbour pointers: ===
+ // === if they are specified in the file: read them from file, ===
+ // === else init them by a call of fill_mel_neighbour() ===
+
+ neigh_set = true;
+ for (int i = 0; i < nElements; i++) {
+ // neighbour information about ith element
+
+ if (read_indices(file, *ind)) {
+ MacroReader::fillMelNeigh(mel[i], mel,
+ VecConv<int,NEIGH,PARTS>::convertVec((*ind),
+ mesh));
+ } else {
+ // setting of neighbours fails
+
+ neigh_set = false;
+ break;
+ }
+ }
+
+ key_def[7] = true;
+ break;
+
+ case 8:
+ // block "element type"
+ fscanf(file, "%*s %*s");
+
+ // === MEL elType ===
+
+ if (dim == 2 || dim == 1)
+ ERROR("there is no element type in 2d and 2d; ignoring data for elType\n");
+
+ for (int i = 0; i < nElements; i++) {
+ result = fscanf(file, "%d", &j);
+ TEST_EXIT(result == 1)
+ ("cannot read elType of element %d in file %s\n", i, filename.c_str());
+ if (dim == 3)
+ (mel)[i]->elType = j;
+ }
+
+ key_def[8] = true;
+ break;
+
+ case 9:
+ // block "projections"
+ {
+ fscanf(file, "%*s");
+
+ int nFaces = mesh->getGeo(FACE);
+ int nEdgesAtBoundary = 0;
+
+ for (k = 1; k < dim; k++)
+ nEdgesAtBoundary += k;
+
+ for (int i = 0; i < nElements; i++) {
+ result = read_indices(file, *ind);
+ TEST_EXIT(result)
+ ("cannot read boundary projector of element %d in file %s\n", i, filename.c_str());
+
+ Projection *projector = Projection::getProjection((*ind)[0]);
+
+ if (projector && projector->getType() == VOLUME_PROJECTION) {
+ mel[i]->setProjection(0, projector);
+ } else { // boundary projection
+ for (j = 0; j < mesh->getGeo(NEIGH); j++) {
+ projector = Projection::getProjection((*ind)[j]);
+ if (projector) {
+ mel[i]->setProjection(j, projector);
+ if (dim > 2) {
+ for (k = 0; k < nEdgesAtBoundary; k++) {
+ int edgeNr = Global::getReferenceElement(dim)->getEdgeOfFace(j, k);
+ mel[i]->setProjection(nFaces + edgeNr, projector);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ key_def[9] = true;
+ break;
+
+ case 10:
+ // block "element region"
+ fscanf(file, "%*s %*s");
+
+ // === MEL regions. ===
+
+ for (int i = 0; i < nElements; i++) {
+ result = fscanf(file, "%d", &j);
+ TEST_EXIT(result == 1)
+ ("cannot read region of element %d in file %s\n", i, filename.c_str());
+ if (j >= 0) {
+ Element *el = mel[i]->getElement();
+ ElementRegion_ED *elementRegion =
+ new ElementRegion_ED(el->getElementData());
+ elementRegion->setRegion(j);
+ el->setElementData(elementRegion);
+ }
+ }
+ key_def[10] = true;
+ break;
+
+ case 11:
+ // block "surface region"
+ fscanf(file, "%*s %*s");
+ for (int i = 0; i < nElements; i++) {
+ result = read_indices(file, *ind);
+ TEST_EXIT(result)
+ ("cannot read surface regions of element %d in file %s\n", i, filename.c_str());
+
+ Element *el = mel[i]->getElement();
+
+ for (j = 0; j < mesh->getGeo(NEIGH); j++) {
+ if ((*ind)[j] >= 0) {
+ SurfaceRegion_ED *surfaceRegion =
+ new SurfaceRegion_ED(el->getElementData());
+ surfaceRegion->setSide(j);
+ surfaceRegion->setRegion((*ind)[j]);
+ el->setElementData(surfaceRegion);
+ }
+ }
+ }
+ key_def[11] = true;
+ break;
+
+ case 12:
+ // line "mesh name"
+ fscanf(file, "%*s %*s %*s");
+ break;
+
+ case 13:
+ // line "time"
+ fscanf(file, "%*s %*s");
+ break;
+ }
+ }
+
+ if (ind)
+ delete ind;
+
+ fclose(file);
+ }
+
+
+ void MacroInfo::dirichletBoundary()
+ {
+ for (int i = 0; i < static_cast<int>(mel.size()); i++)
+ for (int k = 0; k < mesh->getGeo(NEIGH); k++)
+ if (mel[i]->neighbour[k])
+ mel[i]->boundary[k] = INTERIOR;
+ else
+ mel[i]->boundary[k] = DIRICHLET;
+ }
+
+
+ void MacroInfo::fillBoundaryInfo(Mesh *mesh)
+ {
+ FUNCNAME("MacroInfo::fillBoundaryInfo()");
+
+ int i;
+ int nVertices = mesh->getNumberOfVertices();
+ std::deque<MacroElement*>::iterator melIt;
+ std::vector<BoundaryType> bound(nVertices);
+
+ switch (mesh->getDim()) {
+ case 1:
+ break;
+ case 2:
+ for (i = 0; i < nVertices; i++)
+ bound[i] = INTERIOR;
+
+ for (i = 0, melIt = mesh->firstMacroElement();
+ melIt != mesh->endOfMacroElements();
+ ++melIt, i++) {
+ for (int j = 0; j < mesh->getGeo(NEIGH); j++) {
+ if ((*melIt)->getBoundary(j) != INTERIOR) {
+ if ((*melIt)->getBoundary(j) >= DIRICHLET) {
+ int j1 = mel_vertex[i][(j + 1) % 3];
+ int j2 = mel_vertex[i][(j + 2) % 3];
+
+ bound[j1] = max(bound[j1], (*melIt)->getBoundary(j));
+ bound[j2] = max(bound[j2], (*melIt)->getBoundary(j));
+ } else if ((*melIt)->getBoundary(j) <= NEUMANN) {
+ int j1 = mel_vertex[i][(j + 1) % 3];
+ int j2 = mel_vertex[i][(j + 2) % 3];
+
+ if (bound[j1] != INTERIOR)
+ bound[j1] = max(bound[j1], (*melIt)->getBoundary(j));
+ else
+ bound[j1] = (*melIt)->getBoundary(j);
+
+ if (bound[j2] != INTERIOR)
+ bound[j2] = max(bound[j2], (*melIt)->getBoundary(j));
+ else
+ bound[j2] = (*melIt)->getBoundary(j);
+ }
+ }
+ }
+ }
+
+ for (i = 0, melIt = mesh->firstMacroElement();
+ melIt != mesh->endOfMacroElements();
+ ++melIt, i++)
+ for (int j = 0; j < mesh->getGeo(VERTEX); j++)
+ (*melIt)->setBoundary(3 + j, bound[mel_vertex[i][j]]);
+
+ break;
+ case 3:
+ for (i = 0; i < nVertices; i++)
+ bound[i] = INTERIOR;
+
+ for (i = 0, melIt = mesh->firstMacroElement();
+ melIt != mesh->endOfMacroElements();
+ ++melIt, i++) {
+ for (int j = 0; j < mesh->getGeo(NEIGH); j++) {
+ for (int k = 1; k < 4; k++)
+ bound[mel_vertex[i][(j + k) % 4]] =
+ ((*melIt)->getBoundary(j) != INTERIOR) ?
+ newBound((*melIt)->getBoundary(j),
+ bound[mel_vertex[i][(j + k) % 4]]) :
+ bound[mel_vertex[i][(j + k) % 4]];
+ }
+ }
+
+ for (i = 0, melIt = mesh->firstMacroElement();
+ melIt != mesh->endOfMacroElements();
+ ++melIt, i++)
+ for (int j = 0; j < mesh->getGeo(VERTEX); j++)
+ (*melIt)->setBoundary(10 + j, bound[mel_vertex[i][j]]);
+
+ break;
+ default:
+ ERROR_EXIT("invalid dim\n");
+ }
+ }
+
+}
diff --git a/AMDiS/src/MacroInfo.h b/AMDiS/src/MacroInfo.h
new file mode 100644
index 0000000000000000000000000000000000000000..a9cbadcaa835d1898f5da45a91e0be7ead8c6988
--- /dev/null
+++ b/AMDiS/src/MacroInfo.h
@@ -0,0 +1,88 @@
+// ============================================================================
+// == ==
+// == AMDiS - Adaptive multidimensional simulations ==
+// == ==
+// ============================================================================
+// == ==
+// == TU Dresden ==
+// == ==
+// == Institut für Wissenschaftliches Rechnen ==
+// == Zellescher Weg 12-14 ==
+// == 01069 Dresden ==
+// == germany ==
+// == ==
+// ============================================================================
+// == ==
+// == https://gforge.zih.tu-dresden.de/projects/amdis/ ==
+// == ==
+// ============================================================================
+
+/** \file MacroInfo.h */
+
+#ifndef AMDIS_MACROINFO_H
+#define AMDIS_MACROINFO_H
+
+#include <deque>
+#include "AMDiS_fwd.h"
+#include "Global.h"
+
+namespace AMDiS {
+
+ /** \ingroup Input
+ * \brief
+ * Used for reading a macro triangulation
+ */
+ class MacroInfo
+ {
+ public:
+ /// Pointer to the Mesh
+ Mesh *mesh;
+
+ /// list of macro elements
+ std::deque<MacroElement*> mel;
+
+ /// vector of all vertex dofs
+ DegreeOfFreedom **dof;
+
+ /// coords[j][k]: kth coordinate of global vertex j
+ WorldVector<double> *coords;
+
+ /// mel_vertex[i][k]: global index of kth vertex of element i
+ int **mel_vertex;
+
+ /// true, if neighbour information is in macro file
+ bool neigh_set;
+
+ /// true, if boundary information is in macro file
+ bool bound_set;
+
+ public:
+ /** \brief
+ * Reads macro triangulation from ascii file in AMDiS format.
+ * Fills MacroInfo structure.
+ * Called by Mesh::readMacro(), fills missing information
+ */
+ void readAMDiSMacro(std::string filename, Mesh* mesh);
+
+ /// Fills MacroInfo structure and some pointers in mesh
+ void fill(Mesh *mesh, int nElements, int nVertices);
+
+ /// Frees memory of MacroInfo
+ void clear();
+
+ /** \brief
+ * Sets the boundary of all edges/faces with no neigbour to a straight
+ * line/face with dirichlet boundary type
+ */
+ void dirichletBoundary();
+
+ void fillBoundaryInfo(Mesh *mesh);
+
+ protected:
+ /// Reads indices from macro file
+ int read_indices(FILE *file, DimVec<int> &id);
+ };
+
+}
+
+#endif
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index c315629eb71031c0452c59592ca3a1a00c2a3a71..7b2d43313e824c8f29ba1ca4e8b6bc417d3de8b9 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -1,18 +1,16 @@
#include "MacroReader.h"
#include "MacroWriter.h"
#include "MacroElement.h"
+#include "MacroInfo.h"
#include "Boundary.h"
#include "FiniteElemSpace.h"
#include "Mesh.h"
#include <string.h>
#include "FixVec.h"
-#include "FixVecConvert.h"
#include "PeriodicMap.h"
#include "ElInfo.h"
#include "Parameters.h"
#include "DOFIterator.h"
-#include "SurfaceRegion_ED.h"
-#include "ElementRegion_ED.h"
#include "LeafData.h"
#include "VertexVector.h"
#include <map>
@@ -281,7 +279,7 @@ namespace AMDiS {
if (mesh->getNumberOfDOFs(CENTER)) {
for (int i = 0; i < mesh->getNumberOfMacros(); i++)
const_cast<Element*>(mel[i]->getElement())->
- setDOF(mesh->getNode(CENTER), mesh->getDOF(CENTER));
+ setDOF(mesh->getNode(CENTER), mesh->getDof(CENTER));
}
/****************************************************************************/
@@ -321,591 +319,6 @@ namespace AMDiS {
return macroInfo;
}
- /****************************************************************************/
- /* fill macro info structure and some pointers in mesh ... */
- /****************************************************************************/
-
- void MacroInfo::fill(Mesh *pmesh, int ne, int nv)
- {
- FUNCNAME("MacroInfo::fill()");
-
- TEST_EXIT(pmesh)("no mesh\n");
-
- int dim = pmesh->getDim();
- mesh = pmesh;
-
- mesh->setNumberOfElements(ne);
- mesh->setNumberOfLeaves(ne);
- mesh->setNumberOfVertices(nv);
-
- for (int i = 0; i < ne; i++) {
- MacroElement *newMacro = new MacroElement(mesh->getDim());
- mel.push_back(newMacro);
- mesh->addMacroElement(mel[i]);
- }
-
- dof = new DegreeOfFreedom*[nv];
- coords = new WorldVector<double>[nv];
- mel_vertex = new int*[ne];
-
- for (int i = 0; i < ne; i++)
- mel_vertex[i] = new int[mesh->getGeo(VERTEX)];
-
- for (int i = 0; i < nv; i++)
- dof[i] = mesh->getDOF(VERTEX);
-
- for (int i = 0; i < ne; i++) {
- mel[i]->element = mesh->createNewElement();
- mel[i]->index = i;
-
- if (dim == 3)
- mel[i]->elType = 0;
- }
- neigh_set = false;
- bound_set = false;
- }
-
- void MacroInfo::clear()
- {
- for (int i = 0; i < mesh->getNumberOfMacros(); i++)
- delete [] mel_vertex[i];
-
- delete [] mel_vertex;
- delete [] coords;
- coords = NULL;
- delete [] dof;
- dof = NULL;
-
- mesh = NULL;
- neigh_set = false;
- }
-
- /****************************************************************************/
- /****************************************************************************/
- /* tool for reading macro triangulations in ALBERT-format */
- /****************************************************************************/
- /****************************************************************************/
-
- /****************************************************************************/
- /* read_indices() reads dim+1 indices from file into id[0-dim], */
- /* returns true if dim+1 inputs arguments could be read successfully by */
- /* fscanf(), else false */
- /****************************************************************************/
-
- int MacroInfo::read_indices(FILE *file, DimVec<int> &id)
- {
- int dim = mesh->getDim();
-
- for (int i = 0; i <= dim; i++)
- if (fscanf(file, "%d", &id[i]) != 1)
- return(false);
-
- return(true);
- }
-
-#define N_KEYS 14
-#define N_MIN_KEYS 7
- static const char *keys[N_KEYS] = {
- "DIM", // 0
- "DIM_OF_WORLD", // 1
- "number of vertices", // 2
- "number of elements", // 3
- "vertex coordinates", // 4
- "element vertices", // 5
- "element boundaries", // 6
- "element neighbours", // 7
- "element type", // 8
- "projections", // 9
- "element region", // 10
- "surface region", // 11
- "mesh name", // 12
- "time" // 13
- };
-
- static int get_key_no(const char *key)
- {
- for (int i = 0; i < N_KEYS; i++)
- if (!strcmp(keys[i], key))
- return i;
-
- return -1;
- }
-
-#include <ctype.h>
-
- static const char *read_key(const char *line)
- {
- static char key[100];
- char *k = key;
-
- while (isspace(*line))
- line++;
- while ((*k++ = *line++) != ':');
- *--k = '\0';
-
- return(const_cast<const char *>( key));
- }
-
- /****************************************************************************/
- /* read_albert_macro(): */
- /* read macro triangulation from ascii file in ALBERT format */
- /* fills macro_info structure */
- /* called by read_macro(), fills missing information */
- /****************************************************************************/
-
-
- void MacroInfo::readAMDiSMacro(std::string filename, Mesh* mesh)
- {
- FUNCNAME("MacroInfo::readAMDiSMacro()");
-
- int dim;
- int dow, nv, ne, j, k;
- double dbl;
- char line[256];
- int line_no, n_keys, sort_key[N_KEYS], nv_key, ne_key;
- int key_def[N_KEYS] = {0,0,0,0,0,0,0,0,0,0,0,0};
- const char *key;
- DimVec<int> *ind = NULL;
-
- TEST_EXIT(filename != "")("No filename specified!\n");
- TEST_EXIT(filename.length() < static_cast<unsigned int>(127))
- ("can only handle filenames up to 127 characters\n");
-
- FILE *file = fopen(filename.c_str(), "r");
- TEST_EXIT(file)("cannot open file %s\n", filename.c_str());
-
- /****************************************************************************/
- /* looking for all keys in the macro file ... */
- /****************************************************************************/
-
- line_no = n_keys = 0;
- while (fgets(line, 255, file)) {
- line_no++;
- if (!strchr(line, ':'))
- continue;
- key = read_key(line);
- int i_key = get_key_no(key);
- TEST_EXIT(i_key >= 0)
- ("macro file %s must not contain key %s on line %d\n",
- filename.c_str(), key, line_no);
- TEST_EXIT(!key_def[i_key])("key %s defined second time on line %d in file %s\n");
-
- sort_key[n_keys++] = i_key;
- key_def[i_key] = true;
- }
- fclose(file);
-
-
- /*******************************************************************************/
- /* Test, if there is data for every key and if all is defined in right order. */
- /*******************************************************************************/
-
- for (int i_key = 0; i_key < N_MIN_KEYS; i_key++) {
- for (j = 0; j < n_keys; j++)
- if (sort_key[j] == i_key) break;
- TEST_EXIT(j < n_keys)("You do not have specified data for %s in %s\n",
- keys[i_key], filename.c_str());
-
- for (j = 0; j < n_keys; j++)
- if (sort_key[j] == 2) break;
- nv_key = j;
- for (j = 0; j < n_keys; j++)
- if (sort_key[j] == 3) break;
- ne_key = j;
-
- switch (i_key) {
- case 0:
- case 1:
- TEST_EXIT(sort_key[i_key] < 2)
- ("You have to specify DIM or mesh->getGeo(WORLD) before all other data\n");
- break;
- case 4:
- TEST_EXIT(nv_key < i_key)
- ("Before reading data for %s, you have to specify the %s in file\n",
- keys[4], keys[2], filename.c_str());
- break;
- case 5:
- TEST_EXIT(nv_key < i_key && ne_key < i_key)
- ("Before reading data for %s, you have to specify the %s and %s in file %s\n",
- keys[5], keys[3], keys[2], filename.c_str());
- case 6:
- case 7:
- case 8:
- TEST_EXIT(ne_key < i_key)
- ("Before reading data for %s, you have to specify the %s in file %s\n",
- keys[i_key], keys[3], filename.c_str());
- }
- }
-
- for (int i_key = 0; i_key < N_KEYS; i_key++)
- key_def[i_key] = false;
-
- /****************************************************************************/
- /* and now, reading data ... */
- /****************************************************************************/
-
- file = fopen(filename.c_str(), "r");
- TEST_EXIT(file)("cannot open file %s\n", filename.c_str());
-
- int result;
-
- for (int i_key = 0; i_key < n_keys; i_key++) {
-
- switch (sort_key[i_key]) {
-
- case 0:
- // line "DIM"
- result = fscanf(file, "%*s %d", &dim);
- TEST_EXIT(result == 1)("cannot read DIM correctly in file %s\n", filename.c_str());
-
- ind = new DimVec<int>(dim, NO_INIT);
-
- key_def[0] = true;
- break;
-
- case 1:
- // line "DIM_OF_WORLD"
- result = fscanf(file, "%*s %d", &dow);
- TEST_EXIT(result == 1)
- ("cannot read Global::getGeo(WORLD) correctly in file %s\n", filename.c_str());
- TEST_EXIT(dow == Global::getGeo(WORLD))
- ("dimension of world = %d != Global::getGeo(WORLD) = %d\n",
- dow, Global::getGeo(WORLD));
-
- key_def[1] = true;
- break;
-
- case 2:
- // line "number of vertices"
- result = fscanf(file, "%*s %*s %*s %d", &nv);
- TEST_EXIT(result == 1)
- ("cannot read number of vertices correctly in file %s\n", filename.c_str());
- TEST_EXIT(nv > 0)
- ("number of vertices = %d must be bigger than 0\n", nv);
-
- key_def[2] = true;
- if (key_def[3])
- fill(mesh, ne, nv);
- break;
-
- case 3:
- // line "number of elements"
- result = fscanf(file, "%*s %*s %*s %d", &ne);
- TEST_EXIT(result == 1)
- ("cannot read number of elements correctly in file %s\n", filename.c_str());
- TEST_EXIT(ne > 0)
- ("number of elements = %d must be bigger than 0\n", ne);
-
- key_def[3] = true;
- if (key_def[2])
- fill(mesh, ne, nv);
- break;
-
- case 4:
- // block "vertex coordinates"
- fscanf(file, "%*s %*s");
- for (int i = 0; i < nv; i++) {
- for (j = 0; j <Global::getGeo(WORLD) ; j++) {
- result = fscanf(file, "%lf", &dbl);
- TEST_EXIT(result == 1)
- ("error while reading coordinates, check file %s\n", filename.c_str());
- coords[i][j] = dbl;
- }
- }
- key_def[4] = true;
- break;
-
- case 5:
- // block "element vertices"
- fscanf(file, "%*s %*s");
-
- /****************************************************************************/
- /* global number of vertices on a single element */
- /****************************************************************************/
-
- for (int i = 0; i < ne; i++) {
- result = read_indices(file, *ind);
- TEST_EXIT(result)
- ("cannot read vertex indices of element %d in file %s\n", i, filename.c_str());
-
- for (k = 0; k < mesh->getGeo(VERTEX); k++)
- mel_vertex[i][k] = (*ind)[k];
- }
-
- key_def[5] = true;
- break;
-
- case 6:
- // block "element boundaries"
- fscanf(file, "%*s %*s");
-
- /****************************************************************************/
- /* MEL boundary pointers */
- /****************************************************************************/
- for (int i = 0; i < ne; i++) {
- // boundary information of ith element
-
- result = read_indices(file, *ind);
- TEST_EXIT(result)
- ("cannot read boundary type of element %d in file %s\n", i, filename.c_str());
-
- // fill boundary of macro-element
- MacroReader::fillMelBoundary(mesh,
- mel[i],
- VecConv<int,NEIGH,PARTS>::convertVec((*ind), mesh));
- }
-
- this->fillBoundaryInfo(mesh);
-
- bound_set = true;
- key_def[6] = true;
- break;
-
- case 7:
- // block "element neighbours"
- fscanf(file, "%*s %*s");
-
- /****************************************************************************/
- /* fill MEL neighbour pointers: */
- /* if they are specified in the file: read them from file, */
- /* else init them by a call of fill_mel_neighbour() */
- /****************************************************************************/
- neigh_set = true;
- for (int i = 0; i < ne; i++) {
- // neighbour information about ith element
-
- if (read_indices(file, *ind)) {
- MacroReader::fillMelNeigh(mel[i], mel,
- VecConv<int,NEIGH,PARTS>::convertVec((*ind),
- mesh));
- } else {
- neigh_set = false; /* setting of neighbours fails :-( */
- break;
- }
- }
-
- key_def[7] = true;
- break;
-
- case 8:
- // block "element type"
- fscanf(file, "%*s %*s");
- /****************************************************************************/
- /* MEL elType */
- /****************************************************************************/
-
- if (dim == 2 || dim == 1)
- ERROR("there is no element type in 2d and 2d; ignoring data for elType\n");
-
- for (int i = 0; i < ne; i++) {
- result = fscanf(file, "%d", &j);
- TEST_EXIT(result == 1)
- ("cannot read elType of element %d in file %s\n", i, filename.c_str());
- if (dim == 3)
- (mel)[i]->elType = j;
- }
-
- key_def[8] = true;
- break;
-
- case 9:
- // block "projections"
- {
- fscanf(file, "%*s");
-
- int numFaces = mesh->getGeo(FACE);
- int numEdgesAtBoundary = 0;
-
- for (k = 1; k < dim; k++)
- numEdgesAtBoundary += k;
-
- for (int i = 0; i < ne; i++) {
- result = read_indices(file, *ind);
- TEST_EXIT(result)
- ("cannot read boundary projector of element %d in file %s\n", i, filename.c_str());
-
- Projection *projector = Projection::getProjection((*ind)[0]);
- if (projector && projector->getType() == VOLUME_PROJECTION) {
- mel[i]->setProjection(0, projector);
- } else { // boundary projection
- for (j = 0; j < mesh->getGeo(NEIGH); j++) {
- projector = Projection::getProjection((*ind)[j]);
- if (projector) {
- mel[i]->setProjection(j, projector);
- if (dim > 2) {
- for (k = 0; k < numEdgesAtBoundary; k++) {
- int edgeNr = Global::getReferenceElement(dim)->getEdgeOfFace(j, k);
- mel[i]->setProjection(numFaces + edgeNr, projector);
- }
- }
- }
- }
- }
- }
- }
- key_def[9] = true;
- break;
-
- case 10:
- // block "element region"
- fscanf(file, "%*s %*s");
- /****************************************************************************/
- /* MEL regions */
- /****************************************************************************/
-
- for (int i = 0; i < ne; i++) {
- result = fscanf(file, "%d", &j);
- TEST_EXIT(result == 1)
- ("cannot read region of element %d in file %s\n", i, filename.c_str());
- if (j >= 0) {
- Element *el = mel[i]->getElement();
- ElementRegion_ED *elementRegion =
- new ElementRegion_ED(el->getElementData());
- elementRegion->setRegion(j);
- el->setElementData(elementRegion);
- }
- }
- key_def[10] = true;
- break;
-
- case 11:
- // block "surface region"
- fscanf(file, "%*s %*s");
- for (int i = 0; i < ne; i++) {
- result = read_indices(file, *ind);
- TEST_EXIT(result)
- ("cannot read surface regions of element %d in file %s\n", i, filename.c_str());
-
- Element *el = mel[i]->getElement();
-
- for (j = 0; j < mesh->getGeo(NEIGH); j++) {
- if ((*ind)[j] >= 0) {
- SurfaceRegion_ED *surfaceRegion =
- new SurfaceRegion_ED(el->getElementData());
- surfaceRegion->setSide(j);
- surfaceRegion->setRegion((*ind)[j]);
- el->setElementData(surfaceRegion);
- }
- }
- }
- key_def[11] = true;
- break;
-
- case 12:
- // line "mesh name"
- fscanf(file, "%*s %*s %*s");
- break;
-
- case 13:
- // line "time"
- fscanf(file, "%*s %*s");
- break;
- }
- }
-
- if (ind)
- delete ind;
-
- fclose(file);
- }
-
- /****************************************************************************/
- /* sets the boundary of all edges/faces with no neigbour to a straight */
- /* line/face with Dirichlet boundary type */
- /****************************************************************************/
-
- void MacroInfo::dirichletBoundary()
- {
- for (int i = 0; i < static_cast<int>(mel.size()); i++) {
- for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
- if (mel[i]->neighbour[k])
- mel[i]->boundary[k] = INTERIOR;
- else
- mel[i]->boundary[k] = DIRICHLET;
- }
- }
- }
-
-
- void MacroInfo::fillBoundaryInfo(Mesh *mesh)
- {
- FUNCNAME("MacroInfo::fillBoundaryInfo()");
-
- int i, nv = mesh->getNumberOfVertices();
- std::deque<MacroElement*>::iterator melIt;
- std::vector<BoundaryType> bound(nv);
-
- switch (mesh->getDim()) {
- case 1:
- break;
- case 2:
- for (i = 0; i < nv; i++)
- bound[i] = INTERIOR;
-
- for (i = 0, melIt = mesh->firstMacroElement();
- melIt != mesh->endOfMacroElements();
- ++melIt, i++) {
- for (int j = 0; j < mesh->getGeo(NEIGH); j++) {
- if ((*melIt)->getBoundary(j) != INTERIOR) {
- if ((*melIt)->getBoundary(j) >= DIRICHLET) {
- int j1 = mel_vertex[i][(j + 1) % 3];
- int j2 = mel_vertex[i][(j + 2) % 3];
-
- bound[j1] = max(bound[j1], (*melIt)->getBoundary(j));
- bound[j2] = max(bound[j2], (*melIt)->getBoundary(j));
- } else if ((*melIt)->getBoundary(j) <= NEUMANN) {
- int j1 = mel_vertex[i][(j + 1) % 3];
- int j2 = mel_vertex[i][(j + 2) % 3];
-
- if (bound[j1] != INTERIOR)
- bound[j1] = max(bound[j1], (*melIt)->getBoundary(j));
- else
- bound[j1] = (*melIt)->getBoundary(j);
-
- if (bound[j2] != INTERIOR)
- bound[j2] = max(bound[j2], (*melIt)->getBoundary(j));
- else
- bound[j2] = (*melIt)->getBoundary(j);
- }
- }
- }
- }
-
- for (i = 0, melIt = mesh->firstMacroElement();
- melIt != mesh->endOfMacroElements();
- ++melIt, i++)
- for (int j = 0; j < mesh->getGeo(VERTEX); j++)
- (*melIt)->setBoundary(3 + j, bound[mel_vertex[i][j]]);
-
- break;
- case 3:
- for (i = 0; i < nv; i++)
- bound[i] = INTERIOR;
-
- for (i = 0, melIt = mesh->firstMacroElement();
- melIt != mesh->endOfMacroElements();
- ++melIt, i++) {
- for (int j = 0; j < mesh->getGeo(NEIGH); j++) {
- for (int k = 1; k < 4; k++)
- bound[mel_vertex[i][(j + k) % 4]] =
- ((*melIt)->getBoundary(j) != INTERIOR) ?
- newBound((*melIt)->getBoundary(j),
- bound[mel_vertex[i][(j + k) % 4]]) :
- bound[mel_vertex[i][(j + k) % 4]];
- }
- }
-
- for (i = 0, melIt = mesh->firstMacroElement();
- melIt != mesh->endOfMacroElements();
- ++melIt, i++)
- for (int j = 0; j < mesh->getGeo(VERTEX); j++)
- (*melIt)->setBoundary(10 + j, bound[mel_vertex[i][j]]);
-
- break;
- default:
- ERROR_EXIT("invalid dim\n");
- }
- }
void MacroReader::computeNeighbours(Mesh *mesh)
{
@@ -1024,14 +437,14 @@ namespace AMDiS {
mesh->incrementNumberOfEdges(1);
if (mesh->getNumberOfDOFs(EDGE)) {
- dof = el->setDOF(lnode + i, mesh->getDOF(EDGE));
+ dof = el->setDOF(lnode + i, mesh->getDof(EDGE));
if ((*mel)->getNeighbour(i)) {
Element *neigh =
const_cast<Element*>((*mel)->getNeighbour(i)->getElement());
if (periodic[i])
- neigh->setDOF(lnode + (*mel)->getOppVertex(i), mesh->getDOF(EDGE));
+ neigh->setDOF(lnode + (*mel)->getOppVertex(i), mesh->getDof(EDGE));
else
neigh->setDOF(lnode + (*mel)->getOppVertex(i), dof);
}
@@ -1089,7 +502,7 @@ namespace AMDiS {
lnode + k, (*(mel + i))->getIndex());
const_cast<Element*>((*(mel + i))->getElement())->setDOF(lnode + k,
- mesh->getDOF(FACE));
+ mesh->getDof(FACE));
if (neigh) {
ov = (*(mel + i))->getOppVertex(k);
@@ -1101,7 +514,7 @@ namespace AMDiS {
const_cast<Element*>((*(mel + i))->getNeighbour(k)->getElement());
if (periodic[k])
- neighEl->setDOF(lnode+ov, mesh->getDOF(FACE));
+ neighEl->setDOF(lnode+ov, mesh->getDof(FACE));
else
neighEl->setDOF(lnode+ov, const_cast<int*>((*(mel + i))->getElement()->
getDOF(lnode + k)));
@@ -1283,7 +696,7 @@ namespace AMDiS {
/****************************************************************************/
return false;
} else {
- edge_dof = el->setDOF(node+edge_no,mesh->getDOF(EDGE));
+ edge_dof = el->setDOF(node+edge_no,mesh->getDof(EDGE));
}
}
diff --git a/AMDiS/src/MacroReader.h b/AMDiS/src/MacroReader.h
index 10acb48ad21b36b3f4d188a1bba2a8662bd73cd3..ce6127f05e611e570e01a6beb71360e5ad7f6e5f 100644
--- a/AMDiS/src/MacroReader.h
+++ b/AMDiS/src/MacroReader.h
@@ -93,60 +93,6 @@ namespace AMDiS {
friend class MacroInfo;
};
- /** \ingroup Input
- * \brief
- * Used for reading a macro triangulation
- */
- class MacroInfo
- {
- public:
- /// Pointer to the Mesh
- Mesh *mesh;
-
- /// list of macro elements
- std::deque<MacroElement*> mel;
-
- /// vector of all vertex dofs
- DegreeOfFreedom **dof;
-
- /// coords[j][k]: kth coordinate of global vertex j
- WorldVector<double> *coords;
-
- /// mel_vertex[i][k]: global index of kth vertex of element i
- int **mel_vertex;
-
- /// true, if neighbour information is in macro file
- bool neigh_set;
-
- /// true, if boundary information is in macro file
- bool bound_set;
-
- public:
- /// Fills MacroInfo structure and some pointers in mesh
- void fill(Mesh *mesh, int ne, int nv);
-
- /** \brief
- * Reads macro triangulation from ascii file in AMDiS format.
- * Fills MacroInfo structure.
- * Called by Mesh::readMacro(), fills missing information
- */
- void readAMDiSMacro(std::string filename, Mesh* mesh);
-
- /// Frees memory of MacroInfo
- void clear();
-
- /** \brief
- * Sets the boundary of all edges/faces with no neigbour to a straight
- * line/face with dirichlet boundary type
- */
- void dirichletBoundary();
-
- void fillBoundaryInfo(Mesh *mesh);
-
- protected:
- /// Reads indices from macro file
- int read_indices(FILE *file, DimVec<int> &id);
- };
}
#endif
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index dd691e3d5b9281ff955d9ec7c1f2c2ce59a57a99..6ee53d80ed330eb1a0e15b8fd69efa928e73a256 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -11,6 +11,7 @@
#include "ElementDofIterator.h"
#include "MacroElement.h"
#include "MacroReader.h"
+#include "MacroInfo.h"
#include "Mesh.h"
#include "Traverse.h"
#include "Parameters.h"
@@ -470,9 +471,9 @@ namespace AMDiS {
}
- DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
+ DegreeOfFreedom *Mesh::getDof(GeoIndex position)
{
- FUNCNAME("Mesh::getDOF()");
+ FUNCNAME("Mesh::getDof()");
TEST_EXIT_DBG(position >= CENTER && position <= FACE)
("unknown position %d\n", position);
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index eb265b9c2cf7d63542f7199ed36c0dca1f69867d..b7445eabb4d280a9da4c204f319d790812592da8 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -183,7 +183,7 @@ namespace AMDiS {
* DOFs of all DOFAdmin objects belonging to this mesh.
* The return value is a pointer to the first allocated DOF.
*/
- DegreeOfFreedom *getDOF(GeoIndex position);
+ DegreeOfFreedom *getDof(GeoIndex position);
/// Returns *(\ref admin[i]) of the mesh
inline const DOFAdmin& getDOFAdmin(int i) const
diff --git a/AMDiS/src/ParallelError.h b/AMDiS/src/ParallelError.h
deleted file mode 100644
index df92ed5a1edc37572e53f75dcd22d42f07cd16f7..0000000000000000000000000000000000000000
--- a/AMDiS/src/ParallelError.h
+++ /dev/null
@@ -1,149 +0,0 @@
-// ============================================================================
-// == ==
-// == AMDiS - Adaptive multidimensional simulations ==
-// == ==
-// ============================================================================
-// == ==
-// == TU Dresden ==
-// == ==
-// == Institut f�r Wissenschaftliches Rechnen ==
-// == Zellescher Weg 12-14 ==
-// == 01069 Dresden ==
-// == germany ==
-// == ==
-// ============================================================================
-// == ==
-// == https://gforge.zih.tu-dresden.de/projects/amdis/ ==
-// == ==
-// ============================================================================
-
-/** \file ParallelError.h */
-
-#ifndef AMDIS_PARALLELERROR_H
-#define AMDIS_PARALLELERROR_H
-
-#include "AMDiS_fwd.h"
-#include "Global.h"
-#include "Mesh.h"
-
-namespace AMDiS {
-
- /** \ingroup Common
- * \brief
- * True error calculator
- */
- template<typename T>
- class ParallelError
- {
- public:
- static void setWriteLeafData() { writeInLeafData = true; }
- static void unsetWriteLeafData() { writeInLeafData = false; }
- static bool writeLeafData() { return writeInLeafData; }
-
- static double maxErrAtQp(const AbstractFunction<T, WorldVector<double> >& u,
- const DOFVector<T>& uh,
- const Quadrature* q);
- static double
- H1Err(const AbstractFunction<WorldVector<T>, WorldVector<double> >& grdU,
- const DOFVector<T>& uh,
- //const Quadrature* quad,
- int relErr,
- double* max,
- bool writeLeafData = false,
- int comp = 0);
-
- static double L2Err(const AbstractFunction<T, WorldVector<double> >& u,
- const DOFVector<T>& uh,
- //Quadrature* quad,
- int relErr,
- double* max,
- bool writeLeafData = false,
- int comp = 0);
-
- // methods for traversal
- static int maxErrAtQpFct(ElInfo* elinfo);
- static int H1ErrFct(ElInfo* elinfo);
- static int L2ErrFct(ElInfo* elinfo);
- static int relFct(ElInfo* elinfo);
-
- public:
- static const T& errUFct(const DimVec<double>& lambda);
-
- static const WorldVector<T>& grdErrUFct(const DimVec<double>& lambda);
-
- /** \brief
- *
- */
- class AbstrFctErrU : public AbstractFunction<T, DimVec<double> >
- {
- public:
- AbstrFctErrU() : AbstractFunction<T, DimVec<double> >(0) {};
-
- inline const T& operator()(const DimVec<double> & x) const {
- return ParallelError<T>::errUFct(x);
- };
- };
-
- static AbstrFctErrU errU;
-
- /** \brief
- *
- */
- class AbstrFctGrdErrU : public AbstractFunction<WorldVector<T>, DimVec<double> >
- {
- public:
- AbstrFctGrdErrU() : AbstractFunction<WorldVector<T>, DimVec<double> >(0) {};
-
- inline const WorldVector<T>& operator()(const DimVec<double> & x) const {
- return ParallelError<T> :: grdErrUFct (x);
- };
- };
-
- static AbstrFctGrdErrU grdErrU;
-
- private:
- static ElInfo* elinfo;
- /* static const Parametric* el_parametric; */
- static const FastQuadrature* quadFast;
- static const AbstractFunction<T, WorldVector<double> >* pU;
- static const AbstractFunction<WorldVector<T>, WorldVector<double> >* pGrdU;
- static const BasisFunction* basFct;
- static const DOFVector<T>* errUh;
- static double maxErr;
- static double l2Err2;
- static double l2Norm2;
- static int relative;
- static double relNorm2;
- static double h1Err2;
- static double h1Norm2;
- static bool writeInLeafData;
- static int component;
- };
-
- template<typename T> ElInfo* ParallelError<T>::elinfo = NULL;
- /* template<typename T> const Parametric* ParallelError<T>::el_parametric = NULL; */
- template<typename T> const FastQuadrature* ParallelError<T>::quadFast = NULL;
- template<typename T> const AbstractFunction<T, WorldVector<double> >* ParallelError<T>::pU = NULL;
- template<typename T> const AbstractFunction<WorldVector<T>, WorldVector<double> >* ParallelError<T>::pGrdU = NULL;
- template<typename T> const BasisFunction* ParallelError<T>::basFct = NULL;
- template<typename T> const DOFVector<T>* ParallelError<T>::errUh = NULL;
- template<typename T> double ParallelError<T>::maxErr = 0.0;
- template<typename T> double ParallelError<T>::l2Err2 = 0.0;
- template<typename T> double ParallelError<T>::l2Norm2 = 0.0;
- template<typename T> int ParallelError<T>::relative = 0;
- template<typename T> double ParallelError<T>::relNorm2 = 0.0;
- template<typename T> double ParallelError<T>::h1Err2 = 0.0;
- template<typename T> double ParallelError<T>::h1Norm2 = 0.0;
- template<typename T> typename ParallelError<T>::AbstrFctErrU ParallelError<T>::errU;
- template<typename T> typename ParallelError<T>::AbstrFctGrdErrU ParallelError<T>::grdErrU;
- template<typename T> bool ParallelError<T>::writeInLeafData = false;
- template<typename T> int ParallelError<T>::component = 0;
-
-}
-
-#include "ParallelError.hh"
-
-#endif // AMDIS_ERROR_H
-
-
-
diff --git a/AMDiS/src/ParallelError.hh b/AMDiS/src/ParallelError.hh
deleted file mode 100644
index 6de8f88d2bc8fc3cae0ae110782d9a87b17d1a8a..0000000000000000000000000000000000000000
--- a/AMDiS/src/ParallelError.hh
+++ /dev/null
@@ -1,468 +0,0 @@
-#include "Mesh.h"
-#include "Parametric.h"
-#include "Quadrature.h"
-#include "PartitionElementData.h"
-
-namespace AMDiS {
-
- template<typename T>
- const T& ParallelError<T>::errUFct(const DimVec<double>& lambda)
- {
- WorldVector<double> x;
- // if (el_parametric) {
- // ERROR_EXIT("not yet\n");
- // } else {
- elinfo->coordToWorld(lambda, &x);
- // }
- return((*pU)(x));
- }
-
- template<typename T>
- const WorldVector<T>& ParallelError<T>::grdErrUFct(const DimVec<double>& lambda)
- {
- WorldVector<double> x;
- // if (el_parametric) {
- // ERROR_EXIT("not yet\n");
- // } else {
- elinfo->coordToWorld(lambda, &x);
- // }
- return((*pGrdU)(x));
- }
-
-
- template<typename T>
- int ParallelError<T>::maxErrAtQpFct(ElInfo* el_info)
- {
- int i;
- double err;
- const double *u_vec;
- // Parametric *parametric = el_info->getMesh()->getParametric();
-
- elinfo = el_info;
- // if (parametric && parametric->initElement(el_info)) {
- // el_parametric = parametric;
- // } else {
- // el_parametric = NULL;
- // }
-
- u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL);
-
- //uh_el = errUh->getLocalVector(el_info->getElement(), NULL);
- //uh_vec = quadFast->uhAtQp(uh_el, NULL);
-
- mtl::dense_vector<double> uh_vec;
- errUh->getVecAtQPs(el_info, NULL, quadFast, uh_vec);
-
- int numPoints = quadFast->getNumPoints();
-
- for (i = 0; i < numPoints; i++)
- {
- err = u_vec[i] > uh_vec[i] ? u_vec[i] - uh_vec[i] : uh_vec[i] - u_vec[i];
- maxErr = max(maxErr, err);
- }
-
- return;
- }
-
- template<typename T>
- double ParallelError<T>::maxErrAtQp(const AbstractFunction<T, WorldVector<double> >& u,
- const DOFVector<T>& uh,
- const Quadrature* q)
- {
- FUNCNAME("ParallelError<T>::maxErrAtQp");
- const FiniteElemSpace *fe_space;
-
- if (!(pU = &u))
- {
- ERROR("no function u specified; doing nothing\n");
- return(-1.0);
- }
- if (!(errUh = &uh) || !(fe_space = uh->getFeSpace()))
- {
- ERROR("no discrete function or no fe_space for it; doing nothing\n");
- return(-1.0);
- }
-
- if (!(basFct = fe_space->getBasisFcts()))
- {
- ERROR("no basis functions at discrete solution ; doing nothing\n");
- return(-1.0);
- }
-
- int dim = fe_space->getMesh()->getDim();
-
- if (!q)
- q = Quadrature::provideQuadrature(dim,
- 2*fe_space->getBasisFcts()->getDegree()
- - 2);
- quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI);
-
- maxErr = 0.0;
-
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
- Mesh::FILL_COORDS |
- Mesh::CALL_LEAF_EL);
- while(elInfo) {
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (elInfo->getElement()->getElementData(PARTITION_ED));
- TEST_EXIT(partitionData)("no partition data\n");
- if(partitionData && partitionData->getPartitionStatus() == IN) {
- maxErrAtQpFct(elInfo);
- }
- elInfo = stack.traverseNext(elInfo);
- }
-
- // fe_space->getMesh()->traverse(-1,
- // Mesh::FILL_COORDS |
- // Mesh::CALL_LEAF_EL,
- // maxErrAtQpFct);
-
- return(maxErr);
- }
-
- template<typename T>
- int ParallelError<T>::H1ErrFct(ElInfo* el_info)
- {
- int i, j;
- double err, err_2, h1_err_el, norm_el, norm2, det, exact;
- //int dim = el_info->getMesh()->getDim();
- //const DimVec<WorldVector<double> > &Lambda = el_info->getGrdLambda();
- //const T *uh_el;
- const WorldVector<double> *grdu_vec, *grduh_vec;
- // const Parametric *parametric = el_info->getMesh()->getParametric();
-
- elinfo = el_info;
- // if (parametric && parametric->initElement(el_info)) {
- // el_parametric = parametric;
- // }
- // else {
- // el_parametric = NULL;
- // }
-
- grdu_vec = quadFast->getQuadrature()->grdFAtQp(grdErrU, NULL);
-
- //uh_el = errUh->getLocalVector(el_info->getElement(), NULL);
-
- // if (el_parametric) {
- // el_parametric->grd_lambda(el_info, NULL, 1, NULL, &Lambda, &det);
- // }
- // else {
- det = el_info->getDet();
- // }
-
- //grduh_vec = quadFast->grdUhAtQp(Lambda, uh_el, NULL);
- grduh_vec = errUh->getGrdAtQPs(elinfo, NULL, quadFast, NULL);
-
- int numPoints = quadFast->getNumPoints();
- int dow = Global::getGeo(WORLD);
-
- for (h1_err_el = i = 0; i < numPoints; i++)
- {
- for (err_2 = j = 0; j < dow; j++)
- {
- err = grdu_vec[i][j] - grduh_vec[i][j];
- err_2 += sqr(err);
- }
- h1_err_el += quadFast->getWeight(i)*err_2;
- }
-
- exact = det*h1_err_el;
- h1Err2 += exact;
- maxErr = max(maxErr, exact);
-
- if(writeInLeafData)
- el_info->getElement()->setEstimation(exact, component);
-
- if (relative)
- {
- for (norm_el = i = 0; i < numPoints; i++)
- {
- for (norm2 = j = 0; j < dow; j++)
- norm2 += sqr(grdu_vec[i][j]);
- norm_el += quadFast->getWeight(i)*norm2;
- }
- h1Norm2 += det*norm_el;
- }
-
- return 0;
- }
-
- template<typename T>
- double ParallelError<T>::H1Err(
- const AbstractFunction<WorldVector<T>, WorldVector<double> >& grdU,
- const DOFVector<T>& uh,
- //const Quadrature* quad,
- int relErr,
- double* max,
- bool writeLeafData,
- int comp)
- {
- FUNCNAME("ParallelError<T>::H1Err");
- const FiniteElemSpace *fe_space;
-
- writeInLeafData = writeLeafData;
-
- component = comp;
-
- Quadrature *q = NULL;
-
- pGrdU = &grdU;
-
- errUh = &uh;
-
- if(!(fe_space = uh.getFeSpace()))
- {
- ERROR("no fe_space for uh; doing nothing\n");
- return(0.0);
- }
-
- if (!(basFct = fe_space->getBasisFcts()))
- {
- ERROR("no basis functions at discrete solution ; doing nothing\n");
- return(0.0);
- }
-
- int dim = fe_space->getMesh()->getDim();
- int deg = grdU.getDegree();
- int degree = deg ? deg : 2 * fe_space->getBasisFcts()->getDegree() - 2;
-
- //if (!quad)
- q = Quadrature::provideQuadrature(dim, degree);
- quadFast = FastQuadrature::provideFastQuadrature(basFct,
- *q,
- INIT_GRD_PHI);
-
- relative = relErr;
-
- maxErr = h1Err2 = h1Norm2 = 0.0;
-
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
- Mesh::FILL_COORDS |
- Mesh::CALL_LEAF_EL |
- Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA);
- while(elInfo) {
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (elInfo->getElement()->getElementData(PARTITION_ED));
- TEST_EXIT(partitionData)("no partition data\n");
- if(partitionData && partitionData->getPartitionStatus() == IN) {
- H1ErrFct(elInfo);
- }
- elInfo = stack.traverseNext(elInfo);
- }
-
- // fe_space->getMesh()->traverse(-1,
- // Mesh::FILL_COORDS |
- // Mesh::CALL_LEAF_EL |
- // Mesh::FILL_DET |
- // Mesh::FILL_GRD_LAMBDA,
- // H1ErrFct);
-
- if (relative) {
- relNorm2 = h1Norm2+1.e-15;
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
- Mesh::CALL_LEAF_EL);
- while(elInfo) {
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (elInfo->getElement()->getElementData(PARTITION_ED));
- TEST_EXIT(partitionData)("no partition data\n");
- if(partitionData && partitionData->getPartitionStatus() == IN) {
- relFct(elInfo);
- }
- elInfo = stack.traverseNext(elInfo);
- }
-
- // fe_space->getMesh()->traverse(-1,
- // Mesh::CALL_LEAF_EL,
- // relFct);
-
- h1Err2 /= relNorm2;
- maxErr /= relNorm2;
- }
-
- if (max) {
- double totalMax;
- MPI::COMM_WORLD.Allreduce(&maxErr, &totalMax, 1, MPI_DOUBLE, MPI_MAX);
- *max = totalMax;
- }
-
- double totalH1Err2;
- MPI::COMM_WORLD.Allreduce(&h1Err2, &totalH1Err2, 1, MPI_DOUBLE, MPI_SUM);
-
- return(sqrt(totalH1Err2));
- }
-
- template<typename T>
- int ParallelError<T>::L2ErrFct(ElInfo* el_info)
- {
- int i;
- double err, det, l2_err_el, norm_el, exact;
- const double *u_vec;
- mtl::dense_vector<double> uh_vec;
-
- // Parametric *parametric = const_cast<Parametric*>(el_info->getMesh()->getParametric());
-
- elinfo = el_info;
-
- // if (parametric && parametric->initElement(el_info)) {
- // el_parametric = parametric;
- // }
- // else {
- // el_parametric = NULL;
- // }
-
- u_vec = quadFast->getQuadrature()->fAtQp(errU, NULL);
-
- //uh_el = errUh->getLocalVector(el_info->getElement(), NULL);
- //uh_vec = quadFast->uhAtQp(uh_el, NULL);
-
- errUh->getVecAtQPs(el_info, NULL, quadFast, uh_vec);
-
- // if (el_parametric) {
- // el_parametric->det(el_info, NULL, 1, NULL, &det);
- // }
- // else {
- // det = el_info->calcDet();
- // }
-
- det = el_info->getDet();
-
- int numPoints = quadFast->getNumPoints();
-
- for (l2_err_el = i = 0; i < numPoints; i++)
- {
- err = u_vec[i] - uh_vec[i];
- l2_err_el += quadFast->getWeight(i)*sqr(err);
- }
-
- exact = det*l2_err_el;
- l2Err2 += exact;
- maxErr = max(maxErr, exact);
-
- if(writeInLeafData) {
- el_info->getElement()->setEstimation(exact, component);
- }
-
- if (relative)
- {
- for (norm_el = i = 0; i < numPoints; i++)
- norm_el += quadFast->getWeight(i)*sqr(u_vec[i]);
- l2Norm2 += det*norm_el;
- }
-
- return 0;
- }
-
- template<typename T>
- double ParallelError<T>::L2Err(const AbstractFunction<T, WorldVector<double> >& u,
- const DOFVector<T>& uh,
- //Quadrature* quad,
- int relErr,
- double* max,
- bool writeLeafData,
- int comp)
- {
- FUNCNAME("ParallelError<T>::L2Err()");
-
- const FiniteElemSpace *fe_space;
- Quadrature *q = NULL;
- writeInLeafData = writeLeafData;
- component = comp;
-
- if (!(pU = &u)) {
- ERROR("no function u specified; doing nothing\n");
- return(0.0);
- }
-
- if (!(errUh = &uh) || !(fe_space = uh.getFeSpace())) {
- ERROR("no discrete function or no fe_space for it; doing nothing\n");
- return(0.0);
- }
-
- if (!(basFct = fe_space->getBasisFcts())) {
- ERROR("no basis functions at discrete solution ; doing nothing\n");
- return(0.0);
- }
-
- int dim = fe_space->getMesh()->getDim();
- int deg = u.getDegree();
- int degree = deg ? deg : 2 * fe_space->getBasisFcts()->getDegree() - 2;
-
- //if (!quad)
- q = Quadrature::provideQuadrature(dim, degree);
- quadFast = FastQuadrature::provideFastQuadrature(basFct, *q, INIT_PHI);
-
- relative = relErr;
-
- maxErr = l2Err2 = l2Norm2 = 0.0;
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
- Mesh::FILL_COORDS |
- Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA |
- Mesh::CALL_LEAF_EL);
- while(elInfo) {
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (elInfo->getElement()->getElementData(PARTITION_ED));
- TEST_EXIT(partitionData)("no partition data\n");
- if(partitionData && partitionData->getPartitionStatus() == IN) {
- L2ErrFct(elInfo);
- }
- elInfo = stack.traverseNext(elInfo);
- }
-
- // fe_space->getMesh()->traverse(-1,
- // Mesh::FILL_COORDS |
- // Mesh::FILL_DET |
- // Mesh::FILL_GRD_LAMBDA |
- // Mesh::CALL_LEAF_EL,
- // L2ErrFct);
-
- if (relative) {
- relNorm2 = l2Norm2+1.e-15;
- elInfo = stack.traverseFirst(fe_space->getMesh(), -1,
- Mesh::CALL_LEAF_EL);
- while(elInfo) {
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (elInfo->getElement()->getElementData(PARTITION_ED));
- TEST_EXIT(partitionData)("no partition data\n");
- if(partitionData && partitionData->getPartitionStatus() == IN) {
- relFct(elInfo);
- }
- elInfo = stack.traverseNext(elInfo);
- }
-
- // fe_space->getMesh()->traverse(-1, Mesh::CALL_LEAF_EL, relFct);
- l2Err2 /= relNorm2;
- }
-
- if (max) {
- double totalMax;
- MPI::COMM_WORLD.Allreduce(&maxErr, &totalMax, 1, MPI_DOUBLE, MPI_MAX);
- *max = totalMax;
- }
-
- double totalL2Err2;
- MPI::COMM_WORLD.Allreduce(&l2Err2, &totalL2Err2, 1, MPI_DOUBLE, MPI_SUM);
-
- return(sqrt(totalL2Err2));
- }
-
- template<typename T>
- int ParallelError<T>::relFct(ElInfo* elinfo)
- {
- double exact = elinfo->getElement()->getEstimation(component);
- exact /= relNorm2;
- if(writeInLeafData)
- elinfo->getElement()->setEstimation(exact, component);
- return 0;
- }
-
-}
diff --git a/AMDiS/src/ParallelProblem.cc b/AMDiS/src/ParallelProblem.cc
deleted file mode 100644
index 219f33e38d5d333894b3a0ac94fabda2d8769846..0000000000000000000000000000000000000000
--- a/AMDiS/src/ParallelProblem.cc
+++ /dev/null
@@ -1,1895 +0,0 @@
-#include <queue>
-#include <time.h>
-#include "boost/lexical_cast.hpp"
-#include "ParallelProblem.h"
-#include "ProblemScal.h"
-#include "ProblemVec.h"
-#include "ProblemInstat.h"
-#include "AdaptInfo.h"
-#include "AdaptStationary.h"
-#include "ConditionalEstimator.h"
-#include "ConditionalMarker.h"
-#include "Traverse.h"
-#include "ElInfo.h"
-#include "Element.h"
-#include "MacroElement.h"
-#include "PartitionElementData.h"
-#include "ParMetisPartitioner.h"
-#include "Mesh.h"
-#include "DualTraverse.h"
-#include "MeshStructure.h"
-#include "DOFVector.h"
-#include "FiniteElemSpace.h"
-#include "RefinementManager.h"
-#include "CoarseningManager.h"
-#include "Lagrange.h"
-#include "ElementFileWriter.h"
-#include "MacroWriter.h"
-#include "ValueWriter.h"
-#include "SystemVector.h"
-#include "VtkWriter.h"
-#include "mpi.h"
-
-namespace AMDiS {
-
- using boost::lexical_cast;
-
- bool elementInPartition(ElInfo *elInfo)
- {
- PartitionElementData *elementData = dynamic_cast<PartitionElementData*>
- (elInfo->getElement()->getElementData(PARTITION_ED));
- if (elementData && elementData->getPartitionStatus() == IN) {
- return true;
- } else {
- return false;
- }
- }
-
- class MyDualTraverse : public DualTraverse
- {
- public:
- MyDualTraverse(int coarseLevel)
- : coarseLevel_(coarseLevel)
- {}
-
- bool skipEl1(ElInfo *elInfo)
- {
- PartitionElementData *elementData = dynamic_cast<PartitionElementData*>
- (elInfo->getElement()->getElementData(PARTITION_ED));
- if (elementData) {
- if (elInfo->getElement()->isLeaf() &&
- elementData->getLevel() < coarseLevel_)
- return false;
- if(elementData->getLevel() == coarseLevel_)
- return false;
- }
- return true;
- }
-
- /*bool skipEl2(ElInfo *elInfo)
- {
- return skipEl1(elInfo);
- }*/
- private:
- int coarseLevel_;
- };
-
- // =========================================================================
- // ===== class ParallelProblemBase =========================================
- // =========================================================================
-
- ParallelProblemBase::ParallelProblemBase(std::string name,
- ProblemIterationInterface *iterationIF,
- ProblemTimeInterface *timeIF)
- : iterationIF_(iterationIF),
- timeIF_(timeIF),
- debugMode(0),
- debugServerProcess(false)
- {
- mpiRank = MPI::COMM_WORLD.Get_rank();
- mpiSize = MPI::COMM_WORLD.Get_size();
-
- GET_PARAMETER(0, name + "->debug mode", "%d", &debugMode);
-
- groupIsSet=false;
- if (debugMode) {
- MPI::Group group = MPI::COMM_WORLD.Get_group();
-
- int rankSize = mpiSize - 1;
- int *ranks = new int[rankSize];
- for (int i = 0; i < rankSize; i++)
- ranks[i] = i + 1;
-
- std::cout<<"set amdis-group"<<std::endl;
- amdisGroup = group.Incl(rankSize, ranks);
- groupIsSet=true;
- mpiComm = MPI::COMM_WORLD.Create(amdisGroup);
-
- if (mpiComm != MPI::COMM_NULL) {
- mpiRank = mpiComm.Get_rank();
- mpiSize = mpiComm.Get_size();
- debugServerProcess = false;
- } else {
- debugServerProcess = true;
- }
-
- delete [] ranks;
- } else {
- mpiComm = MPI::COMM_WORLD;
- }
- }
-
- void ParallelProblemBase::exitParallelization(AdaptInfo *adaptInfo)
- {
- if (!timeIF_)
- closeTimestep(adaptInfo);
-#ifdef DEBUG
- std::cout<<"free amdisGroup"<<std::endl;
-#endif
- if(groupIsSet)
- amdisGroup.Free();
-#ifdef DEBUG
- std::cout<<"freeD amdisGroup"<<std::endl;
-#endif
- }
-
- void ParallelProblemBase::closeTimestep(AdaptInfo *adaptInfo)
- {
- if (mpiSize > 1 && doBuildGlobalSolution(adaptInfo)) {
- if (debugMode && mpiRank == 0) {
- // Send adaptInfo inforamtion to debug server
- double sendTime = adaptInfo->getTime();
- double sendTimestep = adaptInfo->getTimestep();
- MPI::COMM_WORLD.Send(&sendTime, 1, MPI_DOUBLE, 0, 100);
- MPI::COMM_WORLD.Send(&sendTimestep, 1, MPI_DOUBLE, 0, 100);
- }
- synchronizeMeshes(adaptInfo);
- exchangeRankSolutions(adaptInfo);
- buildGlobalSolution(adaptInfo);
- }
-
- if (timeIF_)
- timeIF_->closeTimestep(adaptInfo);
- }
-
- Flag ParallelProblemBase::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
- {
- Flag flag;
-
- if (mpiSize > 1 && toDo.isSet(MARK | ADAPT)) {
- flag = iterationIF_->oneIteration(adaptInfo, MARK | ADAPT);
-
- double localWeightSum = setElemWeights(adaptInfo);
- if (doPartitioning(adaptInfo, localWeightSum)) {
- clock_t partitioningStart = clock();
-
- synchronizeMeshes(adaptInfo);
- partitionMesh(adaptInfo);
- refineOverlap(adaptInfo);
- createOverlap(adaptInfo);
- synchronizeMeshes(adaptInfo);
- exchangeDOFVectors(adaptInfo);
- coarsenOutOfPartition(adaptInfo);
-
- clock_t partitioningEnd = clock();
- partitioningTime = TIME_USED(partitioningStart,
- partitioningEnd);
- computationStart = partitioningEnd;
- }
-
- flag |= iterationIF_->oneIteration(adaptInfo, toDo & ~(MARK | ADAPT));
- } else {
- flag = iterationIF_->oneIteration(adaptInfo, toDo);
- }
-
- // synchronize adaption flag
- unsigned long *flagBuffer = new unsigned long[mpiSize];
-
- unsigned long localFlag = flag.getFlags();
-
- mpiComm.Allgather(&localFlag, 1, MPI_UNSIGNED_LONG,
- flagBuffer, 1, MPI_UNSIGNED_LONG);
- for (int i = 0; i < mpiSize; i++)
- flag.setFlag(flagBuffer[i]);
-
- delete [] flagBuffer;
-
- return flag;
- }
-
- // =========================================================================
- // ===== class ParallelProblem =============================================
- // =========================================================================
-
- ParallelProblem::ParallelProblem(std::string name,
- ProblemIterationInterface *iterationIF,
- ProblemTimeInterface *timeIF,
- std::vector<DOFVector<double>*> vectors,
- Mesh *mesh_,
- RefinementManager *rm,
- CoarseningManager *cm)
- : ParallelProblemBase(name, iterationIF, timeIF),
- name_(name),
- mesh(mesh_),
- refinementManager(rm),
- coarseningManager(cm),
- repartitionSteps_(1),
- puEveryTimestep_(false),
- dofVectors_(vectors),
- upperPartThreshold_(1.5),
- lowerPartThreshold_(2.0/3.0),
- globalCoarseGridLevel_(0),
- localCoarseGridLevel_(0),
- globalRefinements_(0),
- adaptiveThresholds_(0),
- thresholdIncFactor_(2.0),
- thresholdDecFactor_(0.5),
- repartTimeFactor_(10.0)
- {
- GET_PARAMETER(0, name_ + "->upper part threshold", "%f",
- &upperPartThreshold_);
- GET_PARAMETER(0, name_ + "->lower part threshold", "%f",
- &lowerPartThreshold_);
- GET_PARAMETER(0, name_ + "->global coarse grid level", "%d",
- &globalCoarseGridLevel_);
- GET_PARAMETER(0, name_ + "->local coarse grid level", "%d",
- &localCoarseGridLevel_);
- GET_PARAMETER(0, name_ + "->global refinements", "%d",
- &globalRefinements_);
-
-
- TEST_EXIT(localCoarseGridLevel_ >= globalCoarseGridLevel_)
- ("local coarse grid level < global coarse grid level\n");
-
- partitioner = new ParMetisPartitioner(mesh, &mpiComm);
-
- GET_PARAMETER(0, name_ + "->adaptive thresholds", "%d",
- &adaptiveThresholds_);
- GET_PARAMETER(0, name_ + "->threshold inc factor", "%f",
- &thresholdIncFactor_);
- GET_PARAMETER(0, name_ + "->threshold dec factor", "%f",
- &thresholdDecFactor_);
- GET_PARAMETER(0, name_ + "->repart time factor", "%f",
- &repartTimeFactor_);
-
-
- TEST_EXIT(lowerPartThreshold_ <= 1.0)("invalid lower part threshold\n");
- TEST_EXIT(upperPartThreshold_ >= 1.0)("invalid upper part threshold\n");
-
- if (adaptiveThresholds_) {
- TEST_EXIT(thresholdDecFactor_ <= 1.0)("invalid threshold dec factor\n");
- TEST_EXIT(thresholdIncFactor_ >= 1.0)("invalid threshold inc factor\n");
- }
- minUpperTH_ = upperPartThreshold_;
- maxLowerTH_ = lowerPartThreshold_;
- }
-
- ParallelProblem::~ParallelProblem()
- {
- delete partitioner;
- }
-
- bool ParallelProblem::doPartitioning(AdaptInfo *adaptInfo, double localWeightSum)
- {
- FUNCNAME("ParallelProblem::doPartitioning()");
-
- double *weightSum = new double[mpiSize];
- int *partArray = new int[mpiSize];
- int part = 0;
-
- mpiComm.Gather(&localWeightSum, 1, MPI_DOUBLE, weightSum, 1, MPI_DOUBLE, 0);
-
- if (mpiRank == 0) {
-
- double average = 0.0;
- for (int i = 0; i < mpiSize; i++)
- average += weightSum[i];
-
- average /= mpiSize;
-
- for (int i = 0; i < mpiSize; i++) {
- if ((weightSum[i] / average) > upperPartThreshold_) {
- part = 1;
- break;
- }
- if ((weightSum[i] / average) < lowerPartThreshold_) {
- part = 1;
- break;
- }
- }
-
- double computationTime = TIME_USED(computationStart, clock());
- if (adaptiveThresholds_) {
-
- bool timeOver = ((computationTime / partitioningTime) >= repartTimeFactor_);
-
- if (part == 1 && !timeOver) {
- // inc thresholds
- upperPartThreshold_ *= thresholdIncFactor_;
- lowerPartThreshold_ /= thresholdIncFactor_;
-
- // avoid repartitioning
- part = 0;
- }
-
- if (part == 0 && timeOver) {
- // dec thresholds
- upperPartThreshold_ *= thresholdDecFactor_;
- lowerPartThreshold_ /= thresholdDecFactor_;
-
- upperPartThreshold_ = max(minUpperTH_, upperPartThreshold_);
- lowerPartThreshold_ = min(maxLowerTH_, lowerPartThreshold_);
- }
- }
-
- for (int i = 0; i < mpiSize; i++) {
- partArray[i] = part;
- }
- }
-
- mpiComm.Scatter(partArray, 1, MPI_INT,
- &part, 1, MPI_INT, 0);
-
- delete [] weightSum;
- delete [] partArray;
-
- return (part == 1);
- }
-
- bool ParallelProblem::doBuildGlobalSolution(AdaptInfo *adaptInfo)
- {
- return true;
- }
-
- void ParallelProblem::partitionMesh(AdaptInfo *adaptInfo)
- {
- static bool initial = true;
-
- if (initial) {
- initial = false;
- partitioner->fillCoarsePartitionVec(&oldPartitionVec);
- partitioner->partition(&elemWeights, INITIAL);
- } else {
- oldPartitionVec = partitionVec;
- partitioner->partition(&elemWeights, ADAPTIVE_REPART, 100.0 /*0.000001*/);
- }
-
- partitioner->fillCoarsePartitionVec(&partitionVec);
- }
-
- void ParallelProblem::refineOverlap(AdaptInfo *adaptInfo)
- {
- int dim = mesh->getDim();
- bool finished = (localCoarseGridLevel_ == 0);
-
- while (!finished) {
- std::map<DegreeOfFreedom, int> inOut; // 1: in, 2: out, 3: border dof
-
- // mark in/out/border dofs
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- while (elInfo) {
- Element *element = elInfo->getElement();
- PartitionElementData *partitionData =
- dynamic_cast<PartitionElementData*>(elInfo->getElement()->
- getElementData(PARTITION_ED));
-
- const DegreeOfFreedom **dofs = element->getDOF();
-
- if (partitionData->getPartitionStatus() == IN) {
- for (int i = 0; i < dim + 1; i++) {
- DegreeOfFreedom dof = dofs[i][0];
- if (inOut[dof] == 2)
- inOut[dof] = 3;
- if (inOut[dof] == 0)
- inOut[dof] = 1;
- }
- } else {
- for (int i = 0; i < dim + 1; i++) {
- DegreeOfFreedom dof = dofs[i][0];
- if (inOut[dof] == 1)
- inOut[dof] = 3;
- if (inOut[dof] == 0)
- inOut[dof] = 2;
- }
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- // refine overlap-border and inner elements
- finished = true;
- bool marked = false;
- elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- while (elInfo) {
- Element *element = elInfo->getElement();
- PartitionElementData *partitionData =
- dynamic_cast<PartitionElementData*>(elInfo->getElement()->getElementData(PARTITION_ED));
-
- int level = partitionData->getLevel();
-
- if (level < localCoarseGridLevel_) {
- if (partitionData->getPartitionStatus() != IN) {
- const DegreeOfFreedom **dofs = element->getDOF();
- for (int i = 0; i < dim + 1; i++) {
- DegreeOfFreedom dof = dofs[i][0];
- if (inOut[dof] == 3) {
- element->setMark(1);
- marked = true;
- if ((level + 1) < localCoarseGridLevel_)
- finished = false;
- break;
- }
- }
- } else {
- element->setMark(1);
- marked = true;
- if ((level + 1) < localCoarseGridLevel_)
- finished = false;
- }
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
- if (marked)
- refinementManager->refineMesh(mesh);
- }
- }
-
- void ParallelProblem::globalRefineOutOfPartition(AdaptInfo *adaptInfo)
- {
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- while (elInfo) {
- PartitionElementData *partitionData =
- dynamic_cast<PartitionElementData*>(elInfo->getElement()->
- getElementData(PARTITION_ED));
- int refinements = globalCoarseGridLevel_ - partitionData->getLevel();
- elInfo->getElement()->setMark(max(0, refinements));
- elInfo = stack.traverseNext(elInfo);
- }
-
- refinementManager->refineMesh(mesh);
- }
-
- void ParallelProblem::coarsenOutOfPartition(AdaptInfo *adaptInfo)
- {
- Flag meshCoarsened = 1;
- while (meshCoarsened.getFlags() != 0) {
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- while (elInfo) {
- Element *element = elInfo->getElement();
- PartitionElementData *partitionData =
- dynamic_cast<PartitionElementData*>
- (element->getElementData(PARTITION_ED));
- if (partitionData->getPartitionStatus() == OUT) {
- int mark = min(0, -partitionData->getLevel() + globalCoarseGridLevel_);
- element->setMark(mark);
- }
- elInfo = stack.traverseNext(elInfo);
- }
- meshCoarsened = coarseningManager->coarsenMesh(mesh);
- }
- mpiComm.Barrier();
- }
-
- void ParallelProblem::exchangeMeshStructureCodes(MeshStructure *structures)
- {
- // every process creates a mesh structure code from its mesh.
- structures[mpiRank].init(mesh);
- const std::vector<unsigned long int>& myCode = structures[mpiRank].getCode();
-
- // broadcast code sizes
- int *codeSize = new int[mpiSize];
- int tmp = static_cast<int>(myCode.size());
- mpiComm.Allgather(&tmp, 1, MPI_INT, codeSize, 1, MPI_INT);
- if (debugMode) {
- // send code sizes also to debug server
- MPI::COMM_WORLD.Gather(&tmp, 1, MPI_INT, NULL, 1, MPI_INT, 0);
- }
-
- // broadcast number of elements
- int *elements = new int[mpiSize];
- tmp = structures[mpiRank].getNumElements();
- mpiComm.Allgather(&tmp, 1, MPI_INT, elements, 1, MPI_INT);
- if (debugMode) {
- // send number of elements also to debug server
- MPI::COMM_WORLD.Gather(&tmp, 1, MPI_INT, NULL, 1, MPI_INT, 0);
- }
-
- // broadcast codes
- int *codeOffset = new int[mpiSize];
- int codeSizeSum = 0;
- for (int rank = 0; rank < mpiSize; rank++) {
- codeOffset[rank] = codeSizeSum;
- codeSizeSum += codeSize[rank];
- }
-
- unsigned long int *code = new unsigned long int[codeSizeSum];
- unsigned long int *localCode = new unsigned long int[codeSize[mpiRank]];
- unsigned long int *ptr;
- std::vector<unsigned long int>::const_iterator it, end = myCode.end();
-
- for (ptr = localCode, it = myCode.begin(); it != end; ++it, ++ptr)
- *ptr = *it;
-
- mpiComm.Allgatherv(localCode, codeSize[mpiRank],
- MPI_UNSIGNED_LONG,
- code, codeSize, codeOffset,
- MPI_UNSIGNED_LONG);
- if (debugMode) {
- // send codes also to debug server
- MPI::COMM_WORLD.Send(localCode, codeSize[mpiRank],
- MPI_UNSIGNED_LONG, 0, 100);
- }
-
- for (int rank = 0; rank < mpiSize; rank++) {
- if (rank != mpiRank) {
- std::vector<unsigned long int> remoteCode;
- unsigned long int *ptr;
- unsigned long int *begin = code + codeOffset[rank];
- unsigned long int *end = begin + codeSize[rank];
- for (ptr = begin; ptr != end; ++ptr) {
- remoteCode.push_back(*ptr);
- }
- structures[rank].init(remoteCode, elements[rank]);
- }
- }
-
- delete [] elements;
- delete [] code;
- delete [] localCode;
- delete [] codeOffset;
- delete [] codeSize;
- }
-
- void ParallelProblem::synchronizeMeshes(AdaptInfo *adaptInfo)
- {
- FUNCNAME("ParallelProblem::synchronizeMeshes()");
-
- MeshStructure *structures = new MeshStructure[mpiSize];
-
- // build composite mesh structure
- exchangeMeshStructureCodes(structures);
-
- // merge codes
- for (int rank = 0; rank < mpiSize; rank++)
- if (rank != mpiRank)
- structures[mpiRank].merge(&structures[rank]);
-
- // build finest mesh on the rank partition
- structures[mpiRank].fitMeshToStructure(mesh,
- refinementManager,
- true);
-
- delete [] structures;
- }
-
-
- bool ParallelProblem::writeElement(ElInfo *elInfo)
- {
- Element *element = elInfo->getElement();
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (element->getElementData(PARTITION_ED));
- TEST_EXIT_DBG(partitionData)("no partition data\n");
- PartitionStatus status = partitionData->getPartitionStatus();
- if (status == IN)
- return true;
- else
- return false;
- }
-
- void ParallelProblem::exchangeRankSolutions(AdaptInfo *adaptInfo,
- Mesh *workMesh,
- std::vector<DOFVector<double>*> rankSolutions)
- {
- FUNCNAME("ParallelProblem::exchangeRankSolutions()");
-
- ParallelProblem::fillVertexPartitions(localCoarseGridLevel_, 1, true, overlapDistance_);
- overlapDistance_.clear();
-
- const FiniteElemSpace *feSpace = rankSolutions[0]->getFeSpace();
- int dim = workMesh->getDim();
- const BasisFunction *basFcts = feSpace->getBasisFcts();
- int numFcts = basFcts->getNumber();
- DegreeOfFreedom *coarseDOFs = new DegreeOfFreedom[numFcts];
- DegreeOfFreedom *fineDOFs = new DegreeOfFreedom[numFcts];
- DOFAdmin *admin = feSpace->getAdmin();
-
- std::vector<std::vector<DegreeOfFreedom> > sendOrder(mpiSize);
- std::vector<std::vector<DegreeOfFreedom> > recvOrder(mpiSize);
-
- elementPartitions_.clear();
-
- int elementPartition = -1;
- Element *coarseElement = NULL;
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(workMesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- Element *element = elInfo->getElement();
-
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (element->getElementData(PARTITION_ED));
-
- if (partitionData) {
- if (partitionData->getLevel() == 0) {
- elementPartition = partitionVec[element->getIndex()];
- }
-
- PartitionStatus status = partitionData->getPartitionStatus();
-
- if (status != OUT) {
- if (partitionData->getLevel() == localCoarseGridLevel_) {
- basFcts->getLocalIndices(element, admin, coarseDOFs);
-
- // collect other partitions element belongs to
- for (int i = 0; i < dim + 1; i++) {
- std::set<int>::iterator setBegin = vertexPartitions[coarseDOFs[i]].begin();
- std::set<int>::iterator setEnd = vertexPartitions[coarseDOFs[i]].end();
- for (std::set<int>::iterator setIt = setBegin; setIt != setEnd; ++setIt) {
- elementPartitions_[element].insert(*setIt);
- }
- }
-
- coarseElement = element;
- }
-
-
- if (element->isLeaf()) {
- basFcts->getLocalIndices(element, admin, fineDOFs);
-
- for (int i = 0; i < numFcts; i++) {
- if (status == OVERLAP) {
- // send dofs
- sendOrder[elementPartition].push_back(fineDOFs[i]);
- }
- if (status == IN) {
- // recv dofs
- TEST_EXIT(elementPartition == mpiRank)("???\n");
- std::set<int>::iterator setBegin = elementPartitions_[coarseElement].begin();
- std::set<int>::iterator setEnd = elementPartitions_[coarseElement].end();
- for (std::set<int>::iterator setIt = setBegin; setIt != setEnd; ++setIt) {
- if (*setIt != mpiRank) {
- recvOrder[*setIt].push_back(fineDOFs[i]);
- }
- }
- }
- }
- }
- }
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- // create send and recv buffers and fill send buffers
- DOFVector<double> *solution = rankSolutions[mpiRank];
-
- std::map<int, double*> sendBuffer;
- std::map<int, double*> recvBuffer;
- std::map<int, int> sendBufferSize;
- std::map<int, int> recvBufferSize;
-
- for (int partition = 0; partition < mpiSize; partition++) {
- if (partition != mpiRank) {
- int sendSize = static_cast<int>(sendOrder[partition].size());
- int recvSize = static_cast<int>(recvOrder[partition].size());
-
- sendBufferSize[partition] = sendSize;
- recvBufferSize[partition] = recvSize;
- if (sendSize > 0) {
- sendBuffer[partition] = new double[sendSize];
- std::vector<DegreeOfFreedom>::iterator dofIt;
- dofIt = sendOrder[partition].begin();
- double *bufferIt, *bufferBegin, *bufferEnd;
- bufferBegin = sendBuffer[partition];
- bufferEnd = bufferBegin + sendSize;
- for (bufferIt = bufferBegin; bufferIt < bufferEnd; ++bufferIt, ++dofIt)
- *bufferIt = (*solution)[*dofIt];
- }
- if (recvSize > 0)
- recvBuffer[partition] = new double[recvSize];
- }
- }
-
- // non-blocking sends
- for (int partition = 0; partition < mpiSize; partition++) {
- if (partition != mpiRank) {
- if (sendBufferSize[partition] > 0) {
- mpiComm.Isend(sendBuffer[partition],
- sendBufferSize[partition],
- MPI_DOUBLE,
- partition,
- 0);
- }
- }
- }
-
- // blocking recieves
- for (int partition = 0; partition < mpiSize; partition++) {
- if (partition != mpiRank) {
- if (recvBufferSize[partition] > 0) {
- mpiComm.Recv(recvBuffer[partition],
- recvBufferSize[partition],
- MPI_DOUBLE,
- partition,
- 0);
- }
- }
- }
-
- // wait for end of communication
- mpiComm.Barrier();
-
- // copy values into rank solutions
- for (int partition = 0; partition < mpiSize; partition++) {
- if (partition != mpiRank) {
- std::vector<DegreeOfFreedom>::iterator dofIt = recvOrder[partition].begin();
- for (int i = 0; i < recvBufferSize[partition]; i++) {
- (*(rankSolutions[partition]))[*dofIt] = recvBuffer[partition][i];
- ++dofIt;
- }
- }
- }
-
- // free send and recv buffers
- for (int partition = 0; partition < mpiSize; partition++) {
- if (partition != mpiRank) {
- if (sendBufferSize[partition] > 0)
- delete [] sendBuffer[partition];
- if (recvBufferSize[partition] > 0)
- delete [] recvBuffer[partition];
- }
- }
-
- delete [] coarseDOFs;
- delete [] fineDOFs;
- }
-
- void ParallelProblem::exchangeDOFVector(AdaptInfo *adaptInfo,
- DOFVector<double> *values)
- {
- partitioner->fillLeafPartitionVec(&oldPartitionVec, &oldPartitionVec);
- partitioner->fillLeafPartitionVec(&partitionVec, &partitionVec);
-
- // === get send and recieve orders ===
- std::vector<std::vector<DegreeOfFreedom> > sendOrder;
- std::vector<std::vector<DegreeOfFreedom> > recvOrder;
- sendOrder.resize(mpiSize);
- recvOrder.resize(mpiSize);
-
- const FiniteElemSpace *feSpace = values->getFeSpace();
- const BasisFunction *basFcts = feSpace->getBasisFcts();
- int numFcts = basFcts->getNumber();
- DegreeOfFreedom *dofs = new DegreeOfFreedom[numFcts];
- DOFAdmin *admin = feSpace->getAdmin();
-
- Mesh *mesh = feSpace->getMesh();
- TraverseStack stack;
- ElInfo *elInfo;
-
- elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- while(elInfo) {
- Element *element = elInfo->getElement();
- int index = element->getIndex();
- int oldPartition = oldPartitionVec[index];
- int newPartition = partitionVec[index];
-
- if (oldPartition != newPartition) {
- // get dof indices
- basFcts->getLocalIndices(element, admin, dofs);
-
- if (oldPartition == mpiRank) {
- // send element values to new partition
- for (int i = 0; i < numFcts; i++)
- sendOrder[newPartition].push_back(dofs[i]);
- }
- if (newPartition == mpiRank) {
- // recv element values from old partition
- for (int i = 0; i < numFcts; i++)
- recvOrder[oldPartition].push_back(dofs[i]);
- }
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- delete [] dofs;
-
- // === create send and recv buffers and fill send buffers ===
- std::map<int, double*> sendBuffer;
- std::map<int, double*> recvBuffer;
- std::map<int, int> sendBufferSize;
- std::map<int, int> recvBufferSize;
-
- int partition;
- for (partition = 0; partition < mpiSize; partition++) {
- if (partition != mpiRank) {
- int sendSize = static_cast<int>(sendOrder[partition].size());
- int recvSize = static_cast<int>(recvOrder[partition].size());
-
- sendBufferSize[partition] = sendSize;
- recvBufferSize[partition] = recvSize;
- if (sendSize > 0) {
- sendBuffer[partition] = new double[sendSize];
- std::vector<DegreeOfFreedom>::iterator dofIt;
- dofIt = sendOrder[partition].begin();
- double *bufferIt, *bufferBegin, *bufferEnd;
- bufferBegin = sendBuffer[partition];
- bufferEnd = bufferBegin + sendSize;
- for (bufferIt = bufferBegin; bufferIt < bufferEnd; ++bufferIt, ++dofIt)
- *bufferIt = (*values)[*dofIt];
- }
- if (recvSize > 0)
- recvBuffer[partition] = new double[recvSize];
- }
- }
-
- // === non-blocking sends ===
- for (partition = 0; partition < mpiSize; partition++) {
- if (partition != mpiRank) {
- if (sendBufferSize[partition] > 0) {
- mpiComm.Isend(sendBuffer[partition],
- sendBufferSize[partition],
- MPI_DOUBLE,
- partition,
- 0);
- }
- }
- }
-
- // === blocking receives ===
- for (partition = 0; partition < mpiSize; partition++) {
- if (partition != mpiRank) {
- if (recvBufferSize[partition] > 0) {
- mpiComm.Recv(recvBuffer[partition],
- recvBufferSize[partition],
- MPI_DOUBLE,
- partition,
- 0);
- }
- }
- }
-
- // === wait for end of MPI communication ===
- mpiComm.Barrier();
-
- // === copy received values into DOFVector ===
- for (partition = 0; partition < mpiSize; partition++) {
- if (partition != mpiRank) {
- std::vector<DegreeOfFreedom>::iterator dofIt = recvOrder[partition].begin();
- for (int i = 0; i < recvBufferSize[partition]; i++) {
- (*values)[*dofIt] = recvBuffer[partition][i];
- ++dofIt;
- }
- }
- }
-
- // === free send and receive buffers ===
- for (partition = 0; partition < mpiSize; partition++) {
- if (partition != mpiRank) {
- if (sendBufferSize[partition] > 0)
- delete [] sendBuffer[partition];
- if (recvBufferSize[partition] > 0)
- delete [] recvBuffer[partition];
- }
- }
- }
-
- void ParallelProblem::buildGlobalSolution(AdaptInfo *adaptInfo,
- std::vector<DOFVector<double>*> rankSolutions,
- DOFVector<double> *globalSolution)
- {
- FUNCNAME("ParallelProblem::buildGlobalSolution()");
-
- const FiniteElemSpace *feSpace = globalSolution->getFeSpace();
- int dim = mesh->getDim();
- const BasisFunction *basFcts = feSpace->getBasisFcts();
- int numFcts = basFcts->getNumber();
- DegreeOfFreedom *coarseDOFs = new DegreeOfFreedom[numFcts];
- DegreeOfFreedom *fineDOFs = new DegreeOfFreedom[numFcts];
- DOFAdmin *admin = feSpace->getAdmin();
- Lagrange *linearFunctions = Lagrange::getLagrange(dim, 1);
-
- MSG("Building global solution\n");
-
- // compute w[DOF][partition]->value
- std::map<DegreeOfFreedom, std::map<int, double> > w;
- std::map<DegreeOfFreedom, std::map<int, double> >::iterator wIt, wBegin, wEnd;
-
- std::map<DegreeOfFreedom, double> sumW;
-
- Element *lastCoarseElement = NULL;
-
- WorldVector<double> worldCoord;
- DimVec<double> baryCoord(dim, NO_INIT);
-
- std::set<int>::iterator partIt, partBegin, partEnd;
-
- std::map<DegreeOfFreedom, bool> visited;
-
- MyDualTraverse dualTraverse(localCoarseGridLevel_);
- ElInfo *elInfo1, *elInfo2;
- ElInfo *large, *small;
-
- bool cont = dualTraverse.traverseFirst(mesh, mesh,
- -1, -1,
- Mesh::CALL_EVERY_EL_PREORDER |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET,
- Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS,
- &elInfo1, &elInfo2,
- &small, &large);
-
- while (cont) {
- Element *element1 = elInfo1->getElement();
- Element *element2 = elInfo2->getElement();
- PartitionElementData *partitionData =
- dynamic_cast<PartitionElementData*>
- (element1->getElementData(PARTITION_ED));
-#ifdef DEBUG
- if(element1->getElementData()==NULL)
- std::cout<<"even no elementData, id:"<<element1->getIndex()<<std::endl;
- if(partitionData==NULL) std::cout<<"elem_id:"<<element1->getIndex()<<std::endl;
-#endif
- TEST_EXIT_DBG(partitionData!=NULL)("PARTITION_ED data not found\n");
- if (partitionData->getPartitionStatus() == IN) {
-
- // get coarse dofs
- if (element1 != lastCoarseElement) {
- basFcts->getLocalIndices(element1, admin, coarseDOFs);
- lastCoarseElement = element1;
- }
-
- if (elementPartitions_[element1].size() > 1) {
- // get fine dofs
- basFcts->getLocalIndices(element2, admin, fineDOFs);
-
- // for all fine DOFs
- for (int i = 0; i < numFcts; i++) {
- if (!visited[fineDOFs[i]]) {
- visited[fineDOFs[i]] = true;
-
- elInfo2->coordToWorld(*(basFcts->getCoords(i)), worldCoord);
- elInfo1->worldToCoord(worldCoord, &baryCoord);
-
- // for all coarse vertex DOFs
- for (int j = 0; j < dim + 1; j++) {
- partBegin = vertexPartitions[coarseDOFs[j]].begin();
- partEnd = vertexPartitions[coarseDOFs[j]].end();
- for (partIt = partBegin; partIt != partEnd; ++partIt) {
- int partition = *partIt/* - 1*/;
- double val = (*(linearFunctions->getPhi(j)))(baryCoord);
- w[fineDOFs[i]][partition] += val;
-
- sumW[fineDOFs[i]] += val;
- }
- }
- }
- }
- }
- }
- cont = dualTraverse.traverseNext(&elInfo1, &elInfo2,
- &small, &large);
- }
-
- delete [] coarseDOFs;
- delete [] fineDOFs;
-
- MSG("PU ...\n");
-
- wBegin = w.begin();
- wEnd = w.end();
-
- for (wIt = wBegin; wIt != wEnd; ++wIt) {
- DegreeOfFreedom dof = wIt->first;
- (*globalSolution)[dof] = 0.0;
- }
-
- for (wIt = wBegin; wIt != wEnd; ++wIt) {
- DegreeOfFreedom dof = wIt->first;
- std::map<int, double>::iterator partIt, partBegin, partEnd;
- partBegin = wIt->second.begin();
- partEnd = wIt->second.end();
-
- for (partIt = partBegin; partIt != partEnd; ++partIt) {
- int partition = partIt->first;
- double wDOF = partIt->second;
- (*globalSolution)[dof] += wDOF / sumW[dof] * (*(rankSolutions[partition]))[dof];
- }
- }
- }
-
- double ParallelProblem::errors2map(std::map<int, double> &errVec,
- int comp,
- bool add)
- {
- int elNum = -1;
- double totalErr, error;
-
- if (!add)
- errVec.clear();
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh,
- -1,
- Mesh::CALL_EVERY_EL_PREORDER);
- totalErr = 0.0;
- while(elInfo) {
- Element *element = elInfo->getElement();
-
- // get partition data
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (element->getElementData(PARTITION_ED));
-
- if(partitionData && partitionData->getPartitionStatus() == IN) {
- if(partitionData->getLevel() == 0) {
- elNum = element->getIndex();
- }
- TEST_EXIT(elNum != -1)("invalid element number\n");
- if(element->isLeaf()) {
- error = element->getEstimation(comp);
- errVec[elNum] += error;
- totalErr += error;
- }
- }
- elInfo = stack.traverseNext(elInfo);
- }
-
- return totalErr;
- }
-
- double ParallelProblem::setElemWeights(AdaptInfo *adaptInfo)
- {
- double localWeightSum = 0.0;
- int elNum = -1;
-
- elemWeights.clear();
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1,
- Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- Element *element = elInfo->getElement();
-
- // get partition data
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (element->getElementData(PARTITION_ED));
-
- if (partitionData && partitionData->getPartitionStatus() == IN) {
- if (partitionData->getLevel() == 0) {
- elNum = element->getIndex();
- }
- TEST_EXIT(elNum != -1)("invalid element number\n");
- if (element->isLeaf()) {
- elemWeights[elNum] += 1.0;
- localWeightSum += 1.0;
- }
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- return localWeightSum;
- }
-
-
- void ParallelProblem::globalRefinements()
- {
- if (globalRefinements_ <= 0)
- return;
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
- while (elInfo) {
- Element *element = elInfo->getElement();
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (element->getElementData(PARTITION_ED));
-
- if (partitionData && partitionData->getPartitionStatus() == IN) {
- element->setMark(globalRefinements_);
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- refinementManager->refineMesh(mesh);
- }
-
-
- void ParallelProblem::createOverlap(int level, int overlap, bool openOverlap,
- std::map<Element*, int> &overlapDistance)
- {
- int dim = mesh->getDim();
-
- // === create dual graph (one common node) and prepare breadth-first search ===
- std::map<DegreeOfFreedom, std::vector<Element*> > vertexElements;
- std::queue<Element*> overlapQueue;
-
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- Element *element = elInfo->getElement();
- PartitionElementData *partitionData =
- dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
-
- if (partitionData) {
- if (partitionData->getLevel() == level) {
- for (int i = 0; i < dim + 1; i++) {
- vertexElements[element->getDOF(i, 0)].push_back(element);
- }
-
- if(partitionData->getPartitionStatus() == IN) {
- overlapDistance[element] = 0;
- overlapQueue.push(element);
- } else {
- overlapDistance[element] = -1; // still unknown
- }
- }
- }
- elInfo = stack.traverseNext(elInfo);
- }
-
- // === breadth-first search on dual graph ===
- std::vector<Element*>::iterator it, itBegin, itEnd;
-
- while (!overlapQueue.empty()) {
- // get first element in queue
- Element *element = overlapQueue.front();
- overlapQueue.pop();
-
- // get distance
- int distance = overlapDistance[element];
- TEST_EXIT(distance >= 0)("invalid distance\n");
-
- if (distance >= overlap)
- continue;
-
- // get all adjacent elements
- for (int i = 0; i < dim + 1; i++) {
- itBegin = vertexElements[element->getDOF(i, 0)].begin();
- itEnd = vertexElements[element->getDOF(i, 0)].end();
- for (it = itBegin; it != itEnd; ++it) {
- if (overlapDistance[*it] == -1) {
- // set distance for new member
- overlapDistance[*it] = distance + 1;
- // push neighbour to queue
- overlapQueue.push(*it);
- // mark as overlap on AMDiS mesh
- PartitionElementData *partitionData =
- dynamic_cast<PartitionElementData*>((*it)->getElementData(PARTITION_ED));
- TEST_EXIT(partitionData)("no partition data\n");
- partitionData->setPartitionStatus(OVERLAP);
- partitionData->descend(*it);
- }
- }
- }
- }
- }
-
- void ParallelProblem::fillVertexPartitions(int level, int overlap,
- bool openOverlap,
- std::map<Element*, int> &overlapDistance)
- {
- int dim = mesh->getDim();
-
- // clear partition dof vector
- vertexPartitions.clear();
-
- // first: partition elements ...
- int partition;
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- Element *element = elInfo->getElement();
- PartitionElementData *partitionData =
- dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
- if (partitionData) {
- if (partitionData->getLevel() == 0) {
- partition = partitionVec[element->getIndex()];
- }
- if (partitionData->getLevel() == level) {
- for (int i = 0; i < dim + 1; i++) {
- vertexPartitions[element->getDOF(i, 0)].insert(partition);
- }
- }
- }
- elInfo = stack.traverseNext(elInfo);
- }
-
- if (overlap > 1 || openOverlap == false) {
- // exchange mesh structure codes
- MeshStructure *structures = new MeshStructure[mpiSize];
- exchangeMeshStructureCodes(structures);
-
- // merge codes
- for (int rank = 0; rank < mpiSize; rank++) {
- if (rank != mpiRank) {
- structures[mpiRank].merge(&structures[rank]);
- }
- }
-
- MeshStructure &compositeStructure = structures[mpiRank];
- compositeStructure.reset();
-
- // get composite indices of local overlap elements
- std::map<int, Element*> indexElement;
- std::vector<int> innerOverlapElements; // not at open overlap boundary
-
- elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- Element *element = elInfo->getElement();
- PartitionElementData *partitionData =
- dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
-
- if (partitionData && partitionData->getLevel() == level) {
- int compositeIndex = compositeStructure.getCurrentElement();
- indexElement[compositeIndex] = element;
- if (partitionData->getPartitionStatus() == OVERLAP) {
- int distance = overlapDistance[element];
- if (distance < overlap || !openOverlap) {
- innerOverlapElements.push_back(compositeIndex);
- }
- }
- }
-
- if (element->isLeaf()) {
- compositeStructure.skipBranch();
- } else {
- compositeStructure.nextElement();
- }
- elInfo = stack.traverseNext(elInfo);
- }
-
- // === exchange 'inner' overlap elements ===
-
- // exchange number of overlap elements
- int *numOverlapElements = new int[mpiSize];
- int tmp = static_cast<int>(innerOverlapElements.size());
- mpiComm.Allgather(&tmp, 1, MPI_INT, numOverlapElements, 1, MPI_INT);
-
- // exchange overlap elements
- int *offset = new int[mpiSize];
- int sum = 0;
- for (int rank = 0; rank < mpiSize; rank++) {
- offset[rank] = sum;
- sum += numOverlapElements[rank];
- }
-
- int *recvBuffer = new int[sum];
- int *sendBuffer = new int[numOverlapElements[mpiRank]];
-
- int *ptr;
- std::vector<int>::iterator elemIt, elemEnd = innerOverlapElements.end();
-
- for (ptr = sendBuffer, elemIt = innerOverlapElements.begin();
- elemIt != elemEnd;
- ++elemIt, ++ptr)
- *ptr = *elemIt;
-
- mpiComm.Allgatherv(sendBuffer, numOverlapElements[mpiRank], MPI_INT,
- recvBuffer, numOverlapElements, offset, MPI_INT);
-
-
- // fill vertexPartitions for 'inner' overlap elements
- for (int rank = 0; rank < mpiSize; rank++) {
- int numElements = numOverlapElements[rank];
- int *elements = recvBuffer + offset[rank];
- for (int el = 0; el < numElements; el++) {
- Element *element = indexElement[elements[el]];
- for (int i = 0; i < dim + 1; i++) {
- vertexPartitions[element->getDOF(i, 0)].insert(rank/* + 1*/);
- }
- }
- }
-
- // free memory
- delete [] structures;
- delete [] recvBuffer;
- delete [] sendBuffer;
- delete [] numOverlapElements;
- }
- }
-
- void ParallelProblem::createOverlap(AdaptInfo *adaptInfo)
- {
- int level = localCoarseGridLevel_, overlap = 1;
- bool openOverlap = true;
-
- ParallelProblem::createOverlap(level, overlap, openOverlap, overlapDistance_);
- }
-
- // =========================================================================
- // ===== class ParallelProblemScal =========================================
- // =========================================================================
-
- ParallelProblemScal::ParallelProblemScal(std::string name,
- ProblemScal *prob,
- ProblemInstatScal *problemInstat,
- std::vector<DOFVector<double>*> vectors)
- : ParallelProblem(name,
- prob,
- problemInstat,
- vectors,
- prob->getMesh(),
- prob->getRefinementManager(),
- prob->getCoarseningManager()),
- problem(prob),
- problemInstat_(problemInstat),
- oldEstimator(NULL),
- oldMarker(NULL),
- rankSolution(mpiSize),
- usersEstimator(NULL),
- usersMarker(NULL)
- {
- for (int i = 0; i < mpiSize; i++) {
- rankSolution[i] = new DOFVector<double>(problem->getFeSpace(), "rank solution");
- }
-
- if (problemInstat_) {
- dofVectors_.push_back(problemInstat_->getOldSolution());
- }
- }
-
- ParallelProblemScal::~ParallelProblemScal()
- {
- for (int i = 0; i < mpiSize; i++) {
- delete rankSolution[i];
- }
- }
-
-
- void ParallelProblemScal::initParallelization(AdaptInfo *adaptInfo)
- {
- if (mpiSize > 1) {
- clock_t partitioningStart = clock();
-
- // create an initial partitioning of the mesh
- partitioner->createPartitionData();
-
- // set the element weights, which are 1 at the very first begin
- setElemWeights(adaptInfo);
- // and now partition the mesh
- partitionMesh(adaptInfo);
-
-
-
-#ifdef DEBUG
- //test mesh
- /*TraverseStack testStack;
- ElInfo* testElInfo = testStack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
- while(testElInfo)
- {
- Element *testElement = testElInfo->getElement();
- PartitionElementData* ped=dynamic_cast<PartitionElementData*> (testElement->getElementData(PARTITION_ED));
- if(ped!=NULL)
- std::cout<<ped->getLevel();
- else
- std::cout<<"none";
- std::cout<<"\t"<<testElement->getIndex()<<"\t"<<testElInfo->getLevel()<<std::endl;
- testElInfo=testStack.traverseNext(testElInfo);
- }*/
- //end test
-#endif
-
- globalRefineOutOfPartition(adaptInfo);
-
- refineOverlap(adaptInfo);
- createOverlap(adaptInfo);
-
- globalRefinements();
-
- oldEstimator = problem->getEstimator();
- oldMarker = problem->getMarker();
-
- ConditionalEstimator *condEstimator = NULL;
- if (usersEstimator) {
- problem->setEstimator(usersEstimator);
- } else {
- condEstimator = new ConditionalEstimator(oldEstimator);
- problem->setEstimator(condEstimator);
- }
-
- if (usersMarker) {
- problem->setMarker(usersMarker);
- } else {
- ConditionalMarker *newMarker = new ConditionalMarker("parallel marker",
- -1,
- oldMarker,
- globalCoarseGridLevel_,
- localCoarseGridLevel_);
- problem->setMarker(newMarker);
- }
-
- if (mpiRank == 0) {
- clock_t partitioningEnd = clock();
- partitioningTime = TIME_USED(partitioningStart, partitioningEnd);
- computationStart = partitioningEnd;
- }
-
- // modify file writers
- std::vector<FileWriterInterface*> fileWriters = problem->getFileWriterList();
- std::vector<FileWriterInterface*>::iterator fwIt, fwBegin, fwEnd;
- fwBegin = fileWriters.begin();
- fwEnd = fileWriters.end();
- for (fwIt = fwBegin; fwIt != fwEnd; ++fwIt) {
- (*fwIt)->setFilename((*fwIt)->getFilename() + "_proc" +
- lexical_cast<std::string>(mpiRank) + "_");
- (*fwIt)->setTraverseProperties(-1, 0, elementInPartition);
- }
- }
- }
-
- void ParallelProblemScal::exitParallelization(AdaptInfo *adaptInfo)
- {
- if (mpiSize > 1) {
- ParallelProblem::exitParallelization(adaptInfo);
-#ifdef DEBUG
- std::cout<<"write data"<<std::endl;
-#endif
- if (!timeIF_)
- problem->writeFiles(adaptInfo, true);
-
-#ifdef DEBUG
- std::cout<<"delete data"<<std::endl;
-#endif
- partitioner->deletePartitionData();
-
- if (!usersEstimator)
- delete problem->getEstimator();
- if (!usersMarker)
- delete problem->getMarker();
-
- problem->setEstimator(oldEstimator);
- problem->setMarker(oldMarker);
- }
- }
-
- int ParallelProblemScal::getNumProblems() { return 1; }
-
- ProblemStatBase *ParallelProblemScal::getProblem(int number)
- {
- return problem;
- }
-
- ProblemStatBase *ParallelProblemScal::getProblem(std::string name)
- {
- return problem;
- }
-
-
- void ParallelProblemScal::exchangeDOFVectors(AdaptInfo *adaptInfo)
- {
- std::vector<DOFVector<double>*>::iterator it, itBegin, itEnd;
- itBegin = dofVectors_.begin();
- itEnd = dofVectors_.end();
- for (it = itBegin; it != itEnd; ++it) {
- ParallelProblem::exchangeDOFVector(adaptInfo, *it);
- }
- }
-
- void ParallelProblemScal::exchangeRankSolutions(AdaptInfo *adaptInfo)
- {
- rankSolution[mpiRank]->copy(*(problem->getSolution()));
- ParallelProblem::exchangeRankSolutions(adaptInfo,
- mesh,
- rankSolution);
- }
-
- void ParallelProblemScal::buildGlobalSolution(AdaptInfo *adaptInfo)
- {
- ParallelProblem::buildGlobalSolution(adaptInfo,
- rankSolution,
- problem->getSolution());
- }
-
- void ParallelProblemScal::coarsenOutOfPartition(AdaptInfo *adaptInfo)
- {
- bool coarsenAllowed = adaptInfo->isCoarseningAllowed(0);
- adaptInfo->allowCoarsening(true, 0);
- ParallelProblem::coarsenOutOfPartition(adaptInfo);
- adaptInfo->allowCoarsening(coarsenAllowed, 0);
- }
-
-
- // =========================================================================
- // ===== class ParallelProblemVec ==========================================
- // =========================================================================
-
- ParallelProblemVec::ParallelProblemVec(std::string name,
- ProblemVec *prob,
- ProblemInstatVec *problemInstat,
- std::vector<DOFVector<double>*> vectors)
- : ParallelProblem(name,
- prob,
- problemInstat,
- vectors,
- prob->getMesh(0),
- prob->getRefinementManager(0),
- prob->getCoarseningManager(0)),
- problem(prob),
- problemInstat_(problemInstat)
- {
- nComponents = problem->getNumComponents();
-
- std::vector<FiniteElemSpace*> feSpaces(nComponents);
- for (int i = 0; i < nComponents; i++) {
- feSpaces[i] = problem->getFeSpace(i);
- }
-
- rankSolution.resize(mpiSize);
- for (int i = 0; i < mpiSize; i++) {
- rankSolution[i] = new SystemVector("rank solution", feSpaces, nComponents);
- for (int j = 0; j < nComponents; j++) {
- rankSolution[i]->setDOFVector(j, new DOFVector<double>(feSpaces[j],
- "rank solution"));
- }
- }
-
- if (problemInstat_) {
- for (int i = 0; i < nComponents; i++) {
- dofVectors_.push_back(problemInstat_->getOldSolution()->getDOFVector(i));
- }
- }
- }
-
- ParallelProblemVec::~ParallelProblemVec()
- {
- for (int i = 0; i < mpiSize; i++) {
- for (int j = 0; j < nComponents; j++) {
- delete rankSolution[i]->getDOFVector(j);
- }
- delete rankSolution[i];
- }
- }
-
- void ParallelProblemVec::initParallelization(AdaptInfo *adaptInfo)
- {
- FUNCNAME("ParallelProblem::initParallelization()");
-
- if (mpiSize > 1) {
- partitioner->createPartitionData();
- setElemWeights(adaptInfo);
- partitionMesh(adaptInfo);
-
- refineOverlap(adaptInfo);
- createOverlap(adaptInfo);
-
- oldEstimator = problem->getEstimators();
- oldMarker = problem->getMarkers();
-
- std::vector<ConditionalEstimator*> condEstimator;
- if (static_cast<int>(usersEstimator.size()) == nComponents) {
- problem->setEstimator(usersEstimator);
- } else {
- for (int i = 0; i < nComponents; i++) {
- condEstimator.push_back(new ConditionalEstimator(oldEstimator[i]));
- problem->setEstimator(condEstimator[i], i);
- }
- }
-
- if (static_cast<int>(usersMarker.size()) == nComponents) {
- for (int i = 0; i < nComponents; i++) {
- problem->setMarker(usersMarker[i], i);
- }
- } else {
- TEST_EXIT(static_cast<int>(condEstimator.size()) == nComponents)
- ("use conditional marker only together with conditional estimator\n");
- for (int i = 0; i < nComponents; i++) {
- ConditionalMarker *newMarker =
- new ConditionalMarker("parallel marker",
- i,
- oldMarker[i],
- globalCoarseGridLevel_,
- localCoarseGridLevel_);
- problem->setMarker(newMarker, i);
- }
- }
-
- // modify file writers
- std::vector<FileWriterInterface*> fileWriters = problem->getFileWriterList();
- std::vector<FileWriterInterface*>::iterator fwIt, fwBegin, fwEnd;
- fwBegin = fileWriters.begin();
- fwEnd = fileWriters.end();
- for (fwIt = fwBegin; fwIt != fwEnd; ++fwIt) {
- (*fwIt)->setFilename((*fwIt)->getFilename() + "_proc" +
- lexical_cast<std::string>(mpiComm.Get_rank()) + "_");
- (*fwIt)->setTraverseProperties(-1, 0, elementInPartition);
- }
- }
- }
-
- void ParallelProblemVec::exitParallelization(AdaptInfo *adaptInfo)
- {
- FUNCNAME("ParallelProblem::exitParallelization()");
-
- if (mpiSize > 1) {
- ParallelProblem::exitParallelization(adaptInfo);
-
- if (!timeIF_)
- problem->writeFiles(adaptInfo, true);
-
- partitioner->deletePartitionData();
-
- for (int i = 0; i < nComponents; i++) {
- if (static_cast<int>(usersEstimator.size()) == nComponents)
- delete problem->getEstimator(i);
- if (static_cast<int>(usersMarker.size()) == nComponents)
- delete problem->getMarker(i);
-
- problem->setEstimator(oldEstimator[i], i);
- problem->setMarker(oldMarker[i], i);
- }
-
- usersEstimator.resize(0);
- usersMarker.resize(0);
- }
- }
-
- void ParallelProblemVec::exchangeDOFVectors(AdaptInfo *adaptInfo)
- {
- std::vector<DOFVector<double>*>::iterator it, itBegin, itEnd;
- itBegin = dofVectors_.begin();
- itEnd = dofVectors_.end();
- for (it = itBegin; it != itEnd; ++it) {
- ParallelProblem::exchangeDOFVector(adaptInfo, *it);
- }
- }
-
- void ParallelProblemVec::exchangeRankSolutions(AdaptInfo *adaptInfo)
- {
- rankSolution[mpiRank]->copy(*(problem->getSolution()));
- std::vector<DOFVector<double>*> rankSol(mpiSize);
-
- for (int i = 0; i < nComponents; i++) {
- for (int j = 0; j < mpiSize; j++) {
- rankSol[j] = rankSolution[j]->getDOFVector(i);
- }
- ParallelProblem::exchangeRankSolutions(adaptInfo,
- mesh,
- rankSol);
- }
-
- // In debug mode, mpiRank 0 sends the map partitionVec, which is
- // equal on all ranks, to the debug server.
- if (debugMode && mpiRank == 0) {
- int sendSize = partitionVec.size() * 2;
- MPI::COMM_WORLD.Send(&sendSize, 1, MPI_INT, 0, 100);
-
- int *sendBuffer = new int[sendSize];
- int bufferPos = 0;
- for (std::map<int, int>::iterator it = partitionVec.begin();
- it != partitionVec.end();
- ++it) {
- sendBuffer[bufferPos++] = it->first;
- sendBuffer[bufferPos++] = it->second;
- }
- MPI::COMM_WORLD.Send(sendBuffer, sendSize, MPI_INT, 0, 100);
- delete [] sendBuffer;
- }
- }
-
- void ParallelProblemVec::buildGlobalSolution(AdaptInfo *adaptInfo)
- {
- std::vector<DOFVector<double>*> rankSol(mpiSize);
-
- for (int i = 0; i < nComponents; i++) {
- for (int j = 0; j < mpiSize; j++)
- rankSol[j] = rankSolution[j]->getDOFVector(i);
-
- ParallelProblem::buildGlobalSolution(adaptInfo,
- rankSol,
- problem->getSolution()->getDOFVector(i));
- }
-
- // In debug mode, every rank sends its solution to the debug server.
- if (debugMode) {
- const FiniteElemSpace *feSpace = problem->getFeSpace(0);
- const BasisFunction *basFcts = feSpace->getBasisFcts();
- DOFAdmin *admin = feSpace->getAdmin();
- int numFcts = basFcts->getNumber();
- DegreeOfFreedom *elDOFs = new DegreeOfFreedom[numFcts];
-
- int elementPartition = -1;
- std::vector<double> sendDOFs(0);
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
- while (elInfo) {
- Element *element = elInfo->getElement();
-
- PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
- (element->getElementData(PARTITION_ED));
-
- if (partitionData && (partitionData->getLevel() == 0))
- elementPartition = partitionVec[element->getIndex()];
-
- // On leaf level, each rank assemble its own solution part.
- if (element->isLeaf() && (elementPartition == mpiRank)) {
- basFcts->getLocalIndices(element, admin, elDOFs);
-
- for (int i = 0; i < numFcts; i++)
- sendDOFs.push_back((*problem->getSolution()->getDOFVector(0))[elDOFs[i]]);
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- // First, send number of values the will be send after that.
- int sendSize = sendDOFs.size();
- MPI::COMM_WORLD.Send(&sendSize, 1, MPI_INT, 0, 100);
-
- // Create send buffer of values and send it to the debug server.
- double *sendBuffer = new double[sendSize];
- for (int i = 0; i < sendSize; i++)
- sendBuffer[i] = sendDOFs[i];
- MPI::COMM_WORLD.Send(sendBuffer, sendSize, MPI_DOUBLE, 0, 100);
-
- delete [] sendBuffer;
- delete [] elDOFs;
- }
- }
-
- int ParallelProblemVec::getNumProblems()
- {
- return 1;
- }
-
- ProblemStatBase *ParallelProblemVec::getProblem(int number)
- {
- return problem;
- }
-
- ProblemStatBase *ParallelProblemVec::getProblem(std::string name)
- {
- return problem;
- }
-
- void ParallelProblemVec::coarsenOutOfPartition(AdaptInfo *adaptInfo)
- {
- std::vector<bool> coarsenAllowed(nComponents);
-
- for (int i = 0; i < nComponents; i++) {
- coarsenAllowed[i] = adaptInfo->isCoarseningAllowed(i);
- adaptInfo->allowCoarsening(true, i);
- }
-
- ParallelProblem::coarsenOutOfPartition(adaptInfo);
- for (int i = 0; i < nComponents; i++) {
- adaptInfo->allowCoarsening(coarsenAllowed[i], i);
- }
- }
-
- void ParallelProblemVec::startDebugServer(AdaptInfo *adaptInfo)
- {
- FUNCNAME("ParallelProblemVec::startDebugServer()");
-
- TEST_EXIT(debugMode)("Debug server applicable only in debug mode!\n");
-
- const FiniteElemSpace *feSpace = problem->getFeSpace(0);
- DOFAdmin *admin = feSpace->getAdmin();
- const BasisFunction *basFcts = feSpace->getBasisFcts();
- int numFcts = basFcts->getNumber();
- int mpiWorldSize = MPI::COMM_WORLD.Get_size();
- int *codeSize = new int[mpiWorldSize];
- int *nElements = new int[mpiWorldSize];
- DegreeOfFreedom *elDOFs = new DegreeOfFreedom[numFcts];
- unsigned long int *code;
- int mpiTag = 100;
-
- while (true) {
- /* ===== Communication with closeTimestep() ============== */
-
- // Get current time and timestep
- double recvTime = 0.0;
- double recvTimestep = 0.0;
- MPI::COMM_WORLD.Recv(&recvTime, 1, MPI_DOUBLE, 1, mpiTag);
- MPI::COMM_WORLD.Recv(&recvTimestep, 1, MPI_DOUBLE, 1, mpiTag);
-
- adaptInfo->setTime(recvTime);
- adaptInfo->setTimestep(recvTimestep);
-
-
-
- MeshStructure *structures = new MeshStructure[mpiWorldSize];
-
- /* ===== Communication with exchangeMeshStructureCode() ====== */
-
- // Get mesh structure code size of all meshes
- MPI::COMM_WORLD.Gather(&codeSize, 1, MPI_INT, codeSize, 1, MPI_INT, 0);
-
- // Get number of elements of all meshes
- MPI::COMM_WORLD.Gather(&nElements, 1, MPI_INT, nElements, 1, MPI_INT, 0);
-
- // Get all mesh structure codes
- for (int i = 1; i < mpiWorldSize; i++) {
- code = new unsigned long int[codeSize[i]];
-
- MPI::COMM_WORLD.Recv(code, codeSize[i], MPI_UNSIGNED_LONG, i, mpiTag);
-
- std::vector<unsigned long int> codeVec(codeSize[i]);
- for (int j = 0; j < codeSize[i]; j++)
- codeVec[j] = code[j];
- structures[i].init(codeVec, nElements[i]);
-
- delete [] code;
- }
-
-
- /* ======== Simulate synchronizeMeshes() ======= */
-
- // Build finest mesh starting from rank 1 solution
- for (int rank = 2; rank < mpiWorldSize; rank++)
- structures[1].merge(&structures[rank]);
-
- structures[1].fitMeshToStructure(mesh, refinementManager, false);
-
- delete [] structures;
-
-
- /* ======= Communication with exchangeRankSolution() ======= */
-
- // receive the vector partitionVec from COMM_WORLD.rank = 1 to know
- // the mesh partitioning on localCoarseGridLevel
- int recvSize = 0;
- MPI::COMM_WORLD.Recv(&recvSize, 1, MPI_INT, 1, mpiTag);
-
- int *recvBuffer = new int[recvSize];
- MPI::COMM_WORLD.Recv(recvBuffer, recvSize, MPI_INT, 1, mpiTag);
- partitionVec.clear();
- for (int i = 0; i < recvSize; i += 2)
- partitionVec[recvBuffer[i]] = recvBuffer[i + 1];
- delete [] recvBuffer;
-
-
- /* ======= Communication with buildGlobalSolution() ======= */
-
- // Get final solutions from all ranks
- for (int i = 1; i < mpiSize; i++) {
- // First, get number of DOFs that will be send from rank i
- MPI::COMM_WORLD.Recv(&recvSize, 1, MPI_INT, i, mpiTag);
-
- // Now, receive the solution DOFS of rank i
- double *recvBufferDouble = new double[recvSize];
- MPI::COMM_WORLD.Recv(recvBufferDouble, recvSize, MPI_DOUBLE, i, mpiTag);
-
- // Traverse through the overall mesh, get all elements that are refined from
- // an element with localCoarseGridLevel i, and set there the solutions DOFs.
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
- int elementOwner = -1;
- int dofCount = 0;
- while (elInfo) {
- Element *element = elInfo->getElement();
- if (partitionVec.count(element->getIndex()) > 0) {
- elementOwner = partitionVec[element->getIndex()];
- }
-
- if (element->isLeaf()) {
- basFcts->getLocalIndices(element, admin, elDOFs);
- if (elementOwner == (i - 1)) {
- for (int i = 0; i < numFcts; i++) {
- (*problem->getSolution()->getDOFVector(0))[elDOFs[i]] =
- recvBufferDouble[dofCount + i];
- }
- dofCount += numFcts;
- }
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
-
- delete [] recvBufferDouble;
- }
-
-
- /* ====== Finally, call debug function ====== */
- debugFunction(adaptInfo);
- }
-
- delete [] elDOFs;
- delete [] codeSize;
- delete [] nElements;
- }
-}
diff --git a/AMDiS/src/ParallelProblem.h b/AMDiS/src/ParallelProblem.h
deleted file mode 100644
index 8b271e8dbd59bf8b25c5716bce7288b6549bd7fb..0000000000000000000000000000000000000000
--- a/AMDiS/src/ParallelProblem.h
+++ /dev/null
@@ -1,486 +0,0 @@
-// ============================================================================
-// == ==
-// == AMDiS - Adaptive multidimensional simulations ==
-// == ==
-// ============================================================================
-// == ==
-// == TU Dresden ==
-// == ==
-// == Institut f�r Wissenschaftliches Rechnen ==
-// == Zellescher Weg 12-14 ==
-// == 01069 Dresden ==
-// == germany ==
-// == ==
-// ============================================================================
-// == ==
-// == https://gforge.zih.tu-dresden.de/projects/amdis/ ==
-// == ==
-// ============================================================================
-
-/** \file ParallelProblem.h */
-
-#ifndef AMDIS_PARALLELPROBLEM_H
-#define AMDIS_PARALLELPROBLEM_H
-
-#include <map>
-#include <vector>
-#include <set>
-
-#include "ProblemTimeInterface.h"
-#include "ProblemIterationInterface.h"
-#include "AdaptInfo.h"
-#include "mpi.h"
-
-namespace AMDiS {
-
- /// Interface for parallel problems
- class ParallelProblemInterface
- {
- public:
- virtual ~ParallelProblemInterface() {};
- virtual void initParallelization(AdaptInfo *adaptInfo) = 0;
- virtual void exitParallelization(AdaptInfo *adaptInfo) = 0;
- };
-
- /// Schablonen Klasse
- class ParallelProblemBase : public ParallelProblemInterface,
- public ProblemIterationInterface,
- public ProblemTimeInterface
- {
- public:
- ParallelProblemBase(std::string name,
- ProblemIterationInterface *iterationIF,
- ProblemTimeInterface *timeIF);
-
- virtual ~ParallelProblemBase() {}
-
- /** \brief
- * Must return true, if a new partitioning of the domain (due to unbalanced
- * calculation times) have to be done.
- */
- virtual bool doPartitioning(AdaptInfo *adaptInfo, double localWeightSum) = 0;
-
- virtual bool doBuildGlobalSolution(AdaptInfo *adaptInfo) = 0;
-
- /// Set for each element on the partitioning level the number of leaf elements.
- virtual double setElemWeights(AdaptInfo *adaptInfo) = 0;
-
- virtual void partitionMesh(AdaptInfo *adaptInfo) = 0;
-
- virtual void refineOverlap(AdaptInfo *adaptInfo) = 0;
-
- virtual void globalRefineOutOfPartition(AdaptInfo *adaptInfo) = 0;
-
- virtual void createOverlap(AdaptInfo *adaptInfo) = 0;
-
- virtual void exchangeDOFVectors(AdaptInfo *adaptInfo) = 0;
-
- virtual void coarsenOutOfPartition(AdaptInfo *adaptInfo) = 0;
-
- virtual void synchronizeMeshes(AdaptInfo *adaptInfo) = 0;
-
- virtual void exchangeRankSolutions(AdaptInfo *adaptInfo) = 0;
-
- virtual void buildGlobalSolution(AdaptInfo *adaptInfo) = 0;
-
- virtual void exitParallelization(AdaptInfo *adaptInfo);
-
- virtual void setTime(AdaptInfo *adaptInfo) {
- if (timeIF_)
- timeIF_->setTime(adaptInfo);
- }
-
- virtual void initTimestep(AdaptInfo *adaptInfo) {
- if (timeIF_)
- timeIF_->initTimestep(adaptInfo);
- }
-
- virtual void closeTimestep(AdaptInfo *adaptInfo);
-
- virtual void solveInitialProblem(AdaptInfo *adaptInfo) {
- if (timeIF_)
- timeIF_->solveInitialProblem(adaptInfo);
- }
-
- virtual void transferInitialSolution(AdaptInfo *adaptInfo) {
- if (timeIF_)
- timeIF_->transferInitialSolution(adaptInfo);
- }
-
- virtual void beginIteration(AdaptInfo *adaptInfo) {
- iterationIF_->beginIteration(adaptInfo);
- }
-
- virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION);
-
- virtual void endIteration(AdaptInfo *adaptInfo) {
- iterationIF_->endIteration(adaptInfo);
- }
-
- virtual void startDelayedTimestepCalculation() {}
-
- virtual bool existsDelayedCalculation() {
- return false;
- }
-
- MPI::Intracomm* getMpiComm() {
- return &mpiComm;
- }
-
- bool getDebugServerProcess() {
- return debugServerProcess;
- }
-
- protected:
- ProblemIterationInterface *iterationIF_;
-
- ProblemTimeInterface *timeIF_;
-
- clock_t computationStart;
-
- double partitioningTime;
-
- /// The rank of the current process.
- int mpiRank;
-
- /// Overall number of processes.
- int mpiSize;
-
- /** \brief
- * MPI communicator collected all processes, which should
- * be used for calculation. The Debug procces is not included
- * in this communicator.
- */
- MPI::Intracomm mpiComm;
-
- /// The MPI group, which is used for the communicator \ref mpiComm.
- MPI::Group amdisGroup;
-
- /// Defines the debug mode. If it is 1, a debug server will be started on rank 0.
- int debugMode;
-
- bool groupIsSet;
- /** \brief
- * If true, the current process is the debug server. Otherwise, the
- * processes is a working processes calculating part of the solution.
- */
- bool debugServerProcess;
- };
-
- // =========================================================================
- // ===== class ParallelProblem =============================================
- // =========================================================================
-
- class ParMetisPartitioner;
- class RefinementManager;
- class CoarseningManager;
- class ElInfo;
- class MeshStructure;
-
- class ParallelProblem : public ParallelProblemBase
- {
- public:
- ParallelProblem(std::string name,
- ProblemIterationInterface *iterationIF,
- ProblemTimeInterface *timeIF,
- std::vector<DOFVector<double>*> vectors,
- Mesh *mesh,
- RefinementManager *rm,
- CoarseningManager *cm);
-
- virtual ~ParallelProblem();
-
- virtual bool doPartitioning(AdaptInfo *adaptInfo, double localWeightSum);
-
- virtual bool doBuildGlobalSolution(AdaptInfo *adaptInfo);
-
- virtual void partitionMesh(AdaptInfo *adaptInfo);
-
- virtual void refineOverlap(AdaptInfo *adaptInfo);
-
- virtual void createOverlap(AdaptInfo *adaptInfo);
-
- virtual void globalRefineOutOfPartition(AdaptInfo *adaptInfo);
-
- virtual void coarsenOutOfPartition(AdaptInfo *adaptInfo);
-
- virtual void synchronizeMeshes(AdaptInfo *adaptInfo);
-
- virtual double setElemWeights(AdaptInfo *adaptInfo);
-
- virtual void globalRefinements();
-
- void exchangeRankSolutions(AdaptInfo *adaptInfo,
- Mesh *workMesh,
- std::vector<DOFVector<double>*> rankSolutions);
-
- void exchangeDOFVector(AdaptInfo *adaptInfo,
- DOFVector<double> *vec);
-
- void buildGlobalSolution(AdaptInfo *adaptInfo,
- std::vector<DOFVector<double>*> rankSolutions,
- DOFVector<double> *globalSolution);
-
- void createOverlap(int level, int overlap, bool openOverlap,
- std::map<Element*, int> &overlapDistance);
-
- void fillVertexPartitions(int level, int overlap, bool openOverlap,
- std::map<Element*, int> &overlapDistance);
-
- void setRepartitionSteps(int steps) {
- repartitionSteps_ = steps;
- }
-
- void puEveryTimestep(bool pu) {
- puEveryTimestep_ = pu;
- }
-
- void addDOFVector(DOFVector<double> *vec) {
- dofVectors_.push_back(vec);
- }
-
- /** \brief
- * Every process creates the mesh structure code of its mesh, and all
- * processes broadcast their mesh structure codes.
- */
- void exchangeMeshStructureCodes(MeshStructure *structures);
-
- static bool writeElement(ElInfo *elInfo);
-
- virtual void startDebugServer() {}
-
- virtual void serialize(std::ostream&) {}
-
- virtual void deserialize(std::istream&) {}
-
- protected:
- ///
- double errors2map(std::map<int, double> &errMap, int comp, bool add);
-
- ///
- std::vector<int> iList;
-
- ///
- std::string name_;
-
- /// Mesh of the problem.
- Mesh *mesh;
-
- ///
- RefinementManager *refinementManager;
-
- ///
- CoarseningManager *coarseningManager;
-
- /// Pointer to the paritioner which is used to devide a mesh into partitions.
- ParMetisPartitioner *partitioner;
-
- /** \brief
- * Stores to every coarse element index the number of the partition it
- * corresponds to.
- */
- std::map<int, int> partitionVec;
-
- /** \brief
- * Stores an old partitioning of elements. To every element index the number
- * of the parition it corresponds to is stored.
- */
- std::map<int, int> oldPartitionVec;
-
- /// Weights for the elements, i.e., the number of leaf elements within this element.
- std::map<int, double> elemWeights;
-
- /// Stores to every element the set of partitions it corresponds to.
- std::map<Element*, std::set<int> > elementPartitions_;
-
- /// Stores to every DOF the set of partitions it corresponds to.
- std::map<DegreeOfFreedom, std::set<int> > vertexPartitions;
-
- ///
- int repartitionSteps_;
-
- bool puEveryTimestep_;
-
- std::vector<DOFVector<double>*> dofVectors_;
-
- double upperPartThreshold_;
-
- double lowerPartThreshold_;
-
- int globalCoarseGridLevel_;
-
- int localCoarseGridLevel_;
-
- int globalRefinements_;
-
- std::map<Element*, int> overlapDistance_;
-
- int adaptiveThresholds_;
-
- double thresholdIncFactor_;
-
- double thresholdDecFactor_;
-
- double repartTimeFactor_;
-
- double minUpperTH_;
-
- double maxLowerTH_;
- };
-
- // =========================================================================
- // ===== class ParallelProblemScal =========================================
- // =========================================================================
-
- class FiniteElemSpace;
- class ProblemScal;
- class ProblemInstatScal;
- class Estimator;
- class Marker;
-
- class ParallelProblemScal : public ParallelProblem
- {
- public:
- ParallelProblemScal(std::string name,
- ProblemScal *problem,
- ProblemInstatScal *problemInstat,
- std::vector<DOFVector<double>*> vectors);
-
- virtual ~ParallelProblemScal();
-
- virtual void initParallelization(AdaptInfo *adaptInfo);
-
- virtual void exitParallelization(AdaptInfo *adaptInfo);
-
- virtual void exchangeDOFVectors(AdaptInfo *adaptInfo);
-
- virtual void exchangeRankSolutions(AdaptInfo *adaptInfo);
-
- virtual void buildGlobalSolution(AdaptInfo *adaptInfo);
-
- virtual void coarsenOutOfPartition(AdaptInfo *adaptInfo);
-
- virtual int getNumProblems();
-
- virtual ProblemStatBase *getProblem(int number = 0);
-
- virtual ProblemStatBase *getProblem(std::string name);
-
- void setEstimator(Estimator *est)
- {
- usersEstimator = est;
- }
-
- void setMarker(Marker *marker)
- {
- usersMarker = marker;
- }
-
- inline virtual std::string getName()
- {
- return name_;
- }
-
- protected:
- ProblemScal *problem;
-
- ProblemInstatScal *problemInstat_;
-
- Estimator *oldEstimator;
-
- Marker *oldMarker;
-
- /// Vector of all process' solution DOFVectors.
- std::vector<DOFVector<double>*> rankSolution;
-
- ///
- Estimator *usersEstimator;
-
- ///
- Marker *usersMarker;
- };
-
- // =========================================================================
- // ===== class ParallelProblemVec ==========================================
- // =========================================================================
-
- class ProblemVec;
- class ProblemInstatVec;
- class SystemVector;
-
- class ParallelProblemVec : public ParallelProblem
- {
- public:
- ParallelProblemVec(std::string name,
- ProblemVec *problem,
- ProblemInstatVec *problemInstat,
- std::vector<DOFVector<double>*> vectors);
-
- virtual ~ParallelProblemVec();
-
- virtual void initParallelization(AdaptInfo *adaptInfo);
-
- virtual void exitParallelization(AdaptInfo *adaptInfo);
-
- virtual void exchangeDOFVectors(AdaptInfo *adaptInfo);
-
- virtual void exchangeRankSolutions(AdaptInfo *adaptInfo);
-
- virtual void buildGlobalSolution(AdaptInfo *adaptInfo);
-
- virtual void coarsenOutOfPartition(AdaptInfo *adaptInfo);
-
- virtual int getNumProblems();
-
- virtual ProblemStatBase *getProblem(int number = 0);
-
- virtual ProblemStatBase *getProblem(std::string name);
-
- void setEstimator(std::vector<Estimator*> est)
- {
- usersEstimator = est;
- }
-
- void setMarker(std::vector<Marker*> marker)
- {
- usersMarker = marker;
- }
-
- inline virtual std::string getName()
- {
- return name_;
- }
-
- virtual void startDebugServer(AdaptInfo *adaptInfo);
-
- virtual void debugFunction(AdaptInfo *adaptInfo) {}
-
- protected:
- ///
- ProblemVec *problem;
-
- ///
- ProblemInstatVec *problemInstat_;
-
- ///
- std::vector<Estimator*> oldEstimator;
-
- ///
- std::vector<Marker*> oldMarker;
-
- /// Vector of all process' solution SystemVectors.
- std::vector<SystemVector*> rankSolution;
-
- ///
- std::vector<Estimator*> usersEstimator;
-
- ///
- std::vector<Marker*> usersMarker;
-
- /// Number of components of the vectorial problem.
- int nComponents;
- };
-
-}
-
-#endif
diff --git a/AMDiS/src/PollutionError.cc b/AMDiS/src/PollutionError.cc
deleted file mode 100644
index 9a89e587a0ac1997d9c5dbdb65bd89c36aa45611..0000000000000000000000000000000000000000
--- a/AMDiS/src/PollutionError.cc
+++ /dev/null
@@ -1,148 +0,0 @@
-#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;
- }
-}
-
diff --git a/AMDiS/src/PollutionError.h b/AMDiS/src/PollutionError.h
deleted file mode 100644
index 5d3229699041878a4b420353f9000cab1c82a036..0000000000000000000000000000000000000000
--- a/AMDiS/src/PollutionError.h
+++ /dev/null
@@ -1,101 +0,0 @@
-// ============================================================================
-// == ==
-// == AMDiS - Adaptive multidimensional simulations ==
-// == ==
-// ============================================================================
-// == ==
-// == TU Dresden ==
-// == ==
-// == Institut f�r 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
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index 93b33f3f066fd4c7a96a9597dd60ef98832643b2..e8fa317e0f17e6d0ddf445605b04b20b25eed488 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/AMDiS/src/ProblemScal.cc
@@ -406,13 +406,11 @@ namespace AMDiS {
// === read value file and use it for the mesh values ===
std::string valueFilename("");
GET_PARAMETER(0, mesh->getName() + "->value file name", &valueFilename);
- if (valueFilename.length()) {
+ if (valueFilename.length())
ValueReader::readValue(valueFilename,
mesh,
solution,
mesh->getMacroFileInfo());
- mesh->clearMacroFileInfo();
- }
// === do global refinements ===
if (initFlag.isSet(INIT_GLOBAL_REFINES)) {
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 36a8d2bfb84ac00581e40cc56c36f9735a532d10..03124f81f818e0fcfeda8adb34eb2825b576823e 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -208,13 +208,11 @@ namespace AMDiS {
// === read value file and use it for the mesh values ===
std::string valueFilename("");
GET_PARAMETER(0, meshes[0]->getName() + "->value file name", &valueFilename);
- if (valueFilename.length()) {
+ if (valueFilename.length())
ValueReader::readValue(valueFilename,
meshes[0],
solution->getDOFVector(0),
meshes[0]->getMacroFileInfo());
- meshes[0]->clearMacroFileInfo();
- }
// === do global refinements ===
diff --git a/AMDiS/src/RCNeighbourList.cc b/AMDiS/src/RCNeighbourList.cc
index c72991752542c326ca5577f4703d297802570418..945047d8726393899f6c6bf1ddd2bedbe9e709e8 100644
--- a/AMDiS/src/RCNeighbourList.cc
+++ b/AMDiS/src/RCNeighbourList.cc
@@ -176,7 +176,7 @@ namespace AMDiS {
if (!neighbour || neighbour > coarse_list) {
if (!el->getDOF(node + 2)) {
// face 2
- el->setDOF(node + 2, const_cast<int*>(coarse_mesh->getDOF(FACE)));
+ el->setDOF(node + 2, const_cast<int*>(coarse_mesh->getDof(FACE)));
if (neighbour)
neighbour->el->setDOF(node + coarse_list->oppVertex[0],
const_cast<int*>(el->getDOF(node + 2)));
@@ -187,7 +187,7 @@ namespace AMDiS {
if (!neighbour || neighbour > coarse_list) {
if (!el->getDOF(node + 3)) {
// face 3
- el->setDOF(node + 3, const_cast<int*>(coarse_mesh->getDOF(FACE)));
+ el->setDOF(node + 3, const_cast<int*>(coarse_mesh->getDof(FACE)));
if (neighbour)
neighbour->el->setDOF(node + coarse_list->oppVertex[1],
const_cast<int*>(el->getDOF(node + 3)));
@@ -204,7 +204,7 @@ namespace AMDiS {
if (coarse_mesh->getNumberOfDOFs(CENTER)) {
int node = coarse_mesh->getNode(CENTER);
if (!el->getDOF(node))
- el->setDOF(node, const_cast<int*>(coarse_mesh->getDOF(CENTER)));
+ el->setDOF(node, const_cast<int*>(coarse_mesh->getDof(CENTER)));
}
}
@@ -232,7 +232,7 @@ namespace AMDiS {
/****************************************************************************/
for (int i = 0; i < n_neigh; i++)
if (!rclist[i]->el->getDOF(node))
- rclist[i]->el->setDOF(node, const_cast<int*>(coarse_mesh->getDOF(CENTER)));
+ rclist[i]->el->setDOF(node, const_cast<int*>(coarse_mesh->getDof(CENTER)));
}
}
diff --git a/AMDiS/src/RefinementManager1d.cc b/AMDiS/src/RefinementManager1d.cc
index 196a5cb0e465725f02ac0f938652764d607bb57e..5dff30b39b739e1eb734816938cf32f4ba5e0a94 100644
--- a/AMDiS/src/RefinementManager1d.cc
+++ b/AMDiS/src/RefinementManager1d.cc
@@ -46,7 +46,7 @@ namespace AMDiS {
if (child[0]->getMark() > 0)
doMoreRecursiveRefine = true;
- DegreeOfFreedom *newVertexDOFs = mesh->getDOF(VERTEX);
+ DegreeOfFreedom *newVertexDOFs = mesh->getDof(VERTEX);
child[0]->setDOF(1, newVertexDOFs);
child[1]->setDOF(0, newVertexDOFs);
@@ -69,8 +69,8 @@ namespace AMDiS {
/*--------------------------------------------------------------------------*/
/* there are dofs at the barycenter of the triangles */
/*--------------------------------------------------------------------------*/
- child[0]->setDOF(mesh->getNode(CENTER), const_cast<int*>(mesh->getDOF(CENTER)));
- child[1]->setDOF(mesh->getNode(CENTER), const_cast<int*>(mesh->getDOF(CENTER)));
+ child[0]->setDOF(mesh->getNode(CENTER), const_cast<int*>(mesh->getDof(CENTER)));
+ child[1]->setDOF(mesh->getNode(CENTER), const_cast<int*>(mesh->getDof(CENTER)));
}
/*--------------------------------------------------------------------------*/
diff --git a/AMDiS/src/RefinementManager2d.cc b/AMDiS/src/RefinementManager2d.cc
index 94148d7f9fa032e41ea13df3597568080cc635b5..2be22f5f09ff000e3d309180125d7540b1c7ea37 100644
--- a/AMDiS/src/RefinementManager2d.cc
+++ b/AMDiS/src/RefinementManager2d.cc
@@ -171,15 +171,15 @@ namespace AMDiS {
// === There is one new vertex in the refinement edge. ===
- dof[0] = mesh->getDOF(VERTEX);
+ dof[0] = mesh->getDof(VERTEX);
mesh->incrementNumberOfVertices(1);
mesh->incrementNumberOfEdges(1);
if (mesh->getNumberOfDOFs(EDGE)) {
// There are two additional dofs in the refinement edge.
- dof[1] = mesh->getDOF(EDGE);
- dof[2] = mesh->getDOF(EDGE);
+ dof[1] = mesh->getDof(EDGE);
+ dof[2] = mesh->getDof(EDGE);
}
// === First refine the element. ===
@@ -275,7 +275,7 @@ namespace AMDiS {
mesh->incrementNumberOfElements(2);
if (mesh->getNumberOfDOFs(EDGE)) {
- DegreeOfFreedom* newEdgeDOFs = mesh->getDOF(EDGE);
+ DegreeOfFreedom* newEdgeDOFs = mesh->getDof(EDGE);
// There are dofs in the midpoint of the edges.
child[0]->setDOF(4, newEdgeDOFs);
@@ -294,8 +294,8 @@ namespace AMDiS {
int node = mesh->getNode(CENTER);
// There are dofs at the barycenter of the triangles.
- child[0]->setDOF(node, mesh->getDOF(CENTER));
- child[1]->setDOF(node, mesh->getDOF(CENTER));
+ child[0]->setDOF(node, mesh->getDof(CENTER));
+ child[1]->setDOF(node, mesh->getDof(CENTER));
}
}
diff --git a/AMDiS/src/RefinementManager3d.cc b/AMDiS/src/RefinementManager3d.cc
index 1088d1787d0fb17439c927f0a5956b3898d45db6..1e07ef1c4ece55682ec14d65f2ba67cf45db94ae 100644
--- a/AMDiS/src/RefinementManager3d.cc
+++ b/AMDiS/src/RefinementManager3d.cc
@@ -121,15 +121,15 @@ namespace AMDiS {
/* get new dof for the common face of child0 and child1 */
/****************************************************************************/
- DegreeOfFreedom *newDOF = mesh->getDOF(FACE);
+ DegreeOfFreedom *newDOF = mesh->getDof(FACE);
child[0]->setDOF(node, static_cast<int*>(newDOF));
child[1]->setDOF(node, static_cast<int*>(newDOF));
}
if (mesh->getNumberOfDOFs(CENTER)) {
int node = mesh->getNode(CENTER);
- child[0]->setDOF(node, const_cast<int*>(mesh->getDOF(CENTER)));
- child[1]->setDOF(node, const_cast<int*>(mesh->getDOF(CENTER)));
+ child[0]->setDOF(node, const_cast<int*>(mesh->getDof(CENTER)));
+ child[1]->setDOF(node, const_cast<int*>(mesh->getDof(CENTER)));
}
if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(FACE))
@@ -167,15 +167,15 @@ namespace AMDiS {
node0 = node1 = mesh->getNode(EDGE);
node0 += Tetrahedron::nChildEdge[el_type][0][dir];
node1 += Tetrahedron::nChildEdge[el_type][1][dir];
- DegreeOfFreedom *newDOF = mesh->getDOF(EDGE);
+ DegreeOfFreedom *newDOF = mesh->getDof(EDGE);
(const_cast<Element*>(el->getFirstChild()))->setDOF(node0, newDOF);
(const_cast<Element*>(el->getSecondChild()))->setDOF(node1, newDOF);
}
if (mesh->getNumberOfDOFs(FACE)) {
node0 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir];
- (const_cast<Element*>(el->getFirstChild()))->setDOF(node0, mesh->getDOF(FACE));
+ (const_cast<Element*>(el->getFirstChild()))->setDOF(node0, mesh->getDof(FACE));
node1 = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir];
- (const_cast<Element*>(el->getSecondChild()))->setDOF(node1, mesh->getDOF(FACE));
+ (const_cast<Element*>(el->getSecondChild()))->setDOF(node1, mesh->getDof(FACE));
}
} else { /* if (!neigh || !neigh->child[0]) */
/****************************************************************************/
@@ -316,12 +316,12 @@ namespace AMDiS {
/* get new dof's in the refinement edge */
/****************************************************************************/
- dof[0] = mesh->getDOF(VERTEX);
+ dof[0] = mesh->getDof(VERTEX);
mesh->incrementNumberOfVertices(1);
if (mesh->getNumberOfDOFs(EDGE)) {
- dof[1] = mesh->getDOF(EDGE);
- dof[2] = mesh->getDOF(EDGE);
+ dof[1] = mesh->getDof(EDGE);
+ dof[2] = mesh->getDof(EDGE);
}
for (int i = 0; i < n_neigh; i++)
diff --git a/AMDiS/src/ValueReader.cc b/AMDiS/src/ValueReader.cc
index e5c1851171928c86201e966b106118cbd308ca2d..7235ef6dd2f1195a4f27c3477568b0a02b7a629d 100644
--- a/AMDiS/src/ValueReader.cc
+++ b/AMDiS/src/ValueReader.cc
@@ -1,6 +1,7 @@
#include <cstring>
#include "ValueReader.h"
#include "MacroReader.h"
+#include "MacroInfo.h"
namespace AMDiS {
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index 3e154f5e5937ac5209f605067897d3030bef39dd..24916f89858e97edf66f93345e7cffc267aea3f9 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -186,7 +186,7 @@ namespace AMDiS {
}
ElementVector &c = tmpC[myRank];
- if (num_rows(c) < nPoints)
+ if (num_rows(c) < static_cast<unsigned int>(nPoints))
c.change_dim(nPoints);
c = 0.0;
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 50ca44af551810a601d6a69bfcfed7892c0f1fe9..a27051002547e846b19a07559405635115c39768 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -3,11 +3,12 @@
#include <fstream>
#include <boost/lexical_cast.hpp>
#include <boost/filesystem.hpp>
+#include <limits>
#include "parallel/MeshDistributor.h"
#include "parallel/ParallelDebug.h"
#include "parallel/StdMpi.h"
-#include "ParMetisPartitioner.h"
+#include "parallel/ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
@@ -64,7 +65,7 @@ namespace AMDiS {
}
- void MeshDistributor::initParallelization(AdaptInfo *adaptInfo)
+ void MeshDistributor::initParallelization()
{
FUNCNAME("MeshDistributor::initParallelization()");
@@ -93,9 +94,9 @@ namespace AMDiS {
// create an initial partitioning of the mesh
partitioner->createPartitionData();
// set the element weights, which are 1 at the very first begin
- setElemWeights(adaptInfo);
+ setInitialElementWeights();
// and now partition the mesh
- partitionMesh(adaptInfo);
+ partitionMesh();
#if (DEBUG != 0)
debug::ElementIdxToDofs elMap;
@@ -283,7 +284,7 @@ namespace AMDiS {
}
- void MeshDistributor::exitParallelization(AdaptInfo *adaptInfo)
+ void MeshDistributor::exitParallelization()
{}
@@ -504,20 +505,22 @@ namespace AMDiS {
RankToBoundMap allBound;
for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
- if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
+ if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
+ (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
allBound[it.getRank()].push_back(*it);
for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
- if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
+ if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
+ (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
allBound[it.getRank()].push_back(*it);
for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
- if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
+ if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) ||
+ (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
allBound[it.getRank()].push_back(*it);
// === Check the boundaries and adapt mesh if necessary. ===
- MSG("TEST 1: %d\n", mesh->getDOFAdmin(0).getUsedSize());
#if (DEBUG != 0)
MSG("Run checkAndAdaptBoundary ...\n");
#endif
@@ -535,8 +538,6 @@ namespace AMDiS {
#endif
} while (recvAllValues != 0);
- MSG("TEST 2: %d\n", mesh->getDOFAdmin(0).getUsedSize());
-
#if (DEBUG != 0)
debug::writeMesh(feSpace, -1, "mesh");
#endif
@@ -552,6 +553,11 @@ namespace AMDiS {
// === Update periodic mapping, if there are periodic boundaries. ===
createPeriodicMap();
+
+
+ // === The mesh has changed, so check if it is required to repartition the mesh. ===
+
+ repartitionMesh();
}
@@ -893,11 +899,10 @@ namespace AMDiS {
}
- double MeshDistributor::setElemWeights(AdaptInfo *adaptInfo)
+ void MeshDistributor::setInitialElementWeights()
{
- FUNCNAME("MeshDistributor::setElemWeights()");
+ FUNCNAME("MeshDistributor::setInitialElementWeights()");
- double localWeightSum = 0.0;
elemWeights.clear();
std::string filename = "";
@@ -915,7 +920,6 @@ namespace AMDiS {
infile >> elWeight;
elemWeights[elNum] = elWeight;
- localWeightSum += elWeight;
}
infile.close();
@@ -924,17 +928,13 @@ namespace AMDiS {
ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
elemWeights[elInfo->getElement()->getIndex()] = 1.0;
- localWeightSum++;
-
elInfo = stack.traverseNext(elInfo);
}
}
-
- return localWeightSum;
}
- void MeshDistributor::partitionMesh(AdaptInfo *adaptInfo)
+ void MeshDistributor::partitionMesh()
{
FUNCNAME("MeshDistributor::partitionMesh()");
@@ -951,6 +951,62 @@ namespace AMDiS {
}
+ void MeshDistributor::repartitionMesh()
+ {
+ FUNCNAME("MeshDistributor::repartitionMesh()");
+
+ TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
+ ("Only meshes with one DOFAdmin are supported!\n");
+
+ int nDofs = mesh->getDOFAdmin(0).getUsedDOFs();
+ int repartitioning = 0;
+
+ std::vector<int> nDofsInRank(mpiSize);
+ mpiComm.Gather(&nDofs, 1, MPI_INT, &(nDofsInRank[0]), 1, MPI_INT, 0);
+
+ if (mpiRank == 0) {
+ int nOverallDofs = 0;
+ int minDofs = std::numeric_limits<int>::max();
+ int maxDofs = std::numeric_limits<int>::min();
+ for (int i = 0; i < mpiSize; i++) {
+ nOverallDofs += nDofsInRank[i];
+ minDofs = std::min(minDofs, nDofsInRank[i]);
+ maxDofs = std::max(maxDofs, nDofsInRank[i]);
+ }
+
+ MSG("Overall DOFs: %d Min DOFs: %d Max DOFs: %d\n",
+ nOverallDofs, minDofs, maxDofs);
+
+ if (static_cast<double>(maxDofs) / static_cast<double>(minDofs))
+ repartitioning = 1;
+
+ mpiComm.Bcast(&repartitioning, 1, MPI_INT, 0);
+ } else {
+ mpiComm.Bcast(&repartitioning, 1, MPI_INT, 0);
+ }
+
+
+ if (repartitioning == 0)
+ return;
+
+ return;
+
+ elemWeights.clear();
+
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ elemWeights[elInfo->getMacroElement()->getIndex()]++;
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ partitioner->useLocalGlobalDofMap(&mapLocalGlobalDofs);
+ partitionMesh();
+
+ MSG("DONE\n");
+ }
+
+
void MeshDistributor::createInteriorBoundaryInfo()
{
FUNCNAME("MeshDistributor::createInteriorBoundaryInfo()");
diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h
index d3ce2560b44575bd7a97fbdbd3ad2ab6aa189195..5719ae74ea3821f3c35cb22415ca232be4864c2f 100644
--- a/AMDiS/src/parallel/MeshDistributor.h
+++ b/AMDiS/src/parallel/MeshDistributor.h
@@ -32,7 +32,6 @@
#include "ProblemTimeInterface.h"
#include "ProblemIterationInterface.h"
#include "FiniteElemSpace.h"
-#include "AdaptInfo.h"
#include "parallel/InteriorBoundary.h"
#include "Serializer.h"
#include "BoundaryManager.h"
@@ -80,9 +79,9 @@ namespace AMDiS {
virtual ~MeshDistributor() {}
- void initParallelization(AdaptInfo *adaptInfo);
+ void initParallelization();
- void exitParallelization(AdaptInfo *adaptInfo);
+ void exitParallelization();
void addProblemStat(ProblemVec *probVec);
@@ -106,9 +105,17 @@ namespace AMDiS {
void testForMacroMesh();
/// Set for each element on the partitioning level the number of leaf elements.
- double setElemWeights(AdaptInfo *adaptInfo);
+ void setInitialElementWeights();
- void partitionMesh(AdaptInfo *adaptInfo);
+ /// Creates a partitioning of the mesh for a given weights of the macro elements.
+ void partitionMesh();
+
+ /** \brief
+ * Checks if is required to repartition the mesh. If this is the case, a new
+ * partition will be created and the mesh will be redistributed between the
+ * ranks.
+ */
+ void repartitionMesh();
inline virtual std::string getName()
{
diff --git a/AMDiS/src/ParMetisPartitioner.cc b/AMDiS/src/parallel/ParMetisPartitioner.cc
similarity index 95%
rename from AMDiS/src/ParMetisPartitioner.cc
rename to AMDiS/src/parallel/ParMetisPartitioner.cc
index be345a3433c7d28bc9eaaf89934ec7f416f6ba33..1fc8d5b238010ad78da0fa2908892c641528bcc1 100644
--- a/AMDiS/src/ParMetisPartitioner.cc
+++ b/AMDiS/src/parallel/ParMetisPartitioner.cc
@@ -11,7 +11,8 @@
namespace AMDiS {
- ParMetisMesh::ParMetisMesh(Mesh *mesh, MPI::Intracomm *comm)
+ ParMetisMesh::ParMetisMesh(Mesh *mesh, MPI::Intracomm *comm,
+ std::map<DegreeOfFreedom, DegreeOfFreedom> *mapLocalGlobal)
: dim(mesh->getDim()),
nElements(0),
mpiComm(comm)
@@ -97,7 +98,11 @@ namespace AMDiS {
// write eind entries (element nodes)
for (int i = 0; i < dim + 1; i++) {
- *ptr_eind = element->getDOF(i, 0);
+ if (mapLocalGlobal)
+ *ptr_eind = (*mapLocalGlobal)[element->getDOF(i, 0)];
+ else
+ *ptr_eind = element->getDOF(i, 0);
+
ptr_eind++;
}
@@ -141,12 +146,19 @@ namespace AMDiS {
int ncommonnodes)
: parMetisMesh(parMesh)
{
+ FUNCNAME("ParMetisGraph::ParMetisGraph()");
+
+ TEST_EXIT(parMesh)("No ParMetisMesh defined!\n");
+ TEST_EXIT(comm)("No MPI communicator defined!\n");
+
int numflag = 0;
if (ncommonnodes == -1)
ncommonnodes = parMetisMesh->getDim();
MPI_Comm tmpComm = MPI_Comm(*comm);
+ int mpiRank = MPI::COMM_WORLD.Get_rank();
+ int mpiSize = MPI::COMM_WORLD.Get_size();
ParMETIS_V3_Mesh2Dual(parMetisMesh->getElementDist(),
parMetisMesh->getElementPtr(),
@@ -226,7 +238,7 @@ namespace AMDiS {
if (parMetisMesh)
delete parMetisMesh;
- parMetisMesh = new ParMetisMesh(mesh, mpiComm);
+ parMetisMesh = new ParMetisMesh(mesh, mpiComm, mapLocalGlobal);
int nElements = parMetisMesh->getNumElements();
@@ -280,7 +292,7 @@ namespace AMDiS {
int options[3] = {0, 0, 0}; // default options
int edgecut = -1;
int *part = new int[nElements];
-
+
if (elemWeights) {
// set tpwgts
for (int i = 0; i < mpiSize; i++)
@@ -380,7 +392,7 @@ namespace AMDiS {
// update ParMETIS mesh to new partitioning
if (!parMetisMesh)
- parMetisMesh = new ParMetisMesh(mesh, mpiComm);
+ parMetisMesh = new ParMetisMesh(mesh, mpiComm, mapLocalGlobal);
int mpiRank = mpiComm->Get_rank();
int mpiSize = mpiComm->Get_size();
diff --git a/AMDiS/src/ParMetisPartitioner.h b/AMDiS/src/parallel/ParMetisPartitioner.h
similarity index 93%
rename from AMDiS/src/ParMetisPartitioner.h
rename to AMDiS/src/parallel/ParMetisPartitioner.h
index 5ca73d47534deb693dfe562568e1211dcf1063e8..c860fb29397717a30233e742fc4b1e00ab89c820 100644
--- a/AMDiS/src/ParMetisPartitioner.h
+++ b/AMDiS/src/parallel/ParMetisPartitioner.h
@@ -45,7 +45,8 @@ namespace AMDiS {
class ParMetisMesh
{
public:
- ParMetisMesh(Mesh *mesh, MPI::Intracomm *comm);
+ ParMetisMesh(Mesh *mesh, MPI::Intracomm *comm,
+ std::map<DegreeOfFreedom, DegreeOfFreedom> *mapLocalGlobal);
~ParMetisMesh();
@@ -170,7 +171,8 @@ namespace AMDiS {
ParMetisPartitioner(Mesh *m, MPI::Intracomm *comm)
: mesh(m),
mpiComm(comm),
- parMetisMesh(NULL)
+ parMetisMesh(NULL),
+ mapLocalGlobal(NULL)
{}
void partition(std::map<int, double> *elemWeights,
@@ -194,6 +196,11 @@ namespace AMDiS {
void deletePartitionData();
+ void useLocalGlobalDofMap(std::map<DegreeOfFreedom, DegreeOfFreedom> *m)
+ {
+ mapLocalGlobal = m;
+ }
+
protected:
void distributePartitioning(int *part);
@@ -205,6 +212,8 @@ namespace AMDiS {
MPI::Intracomm *mpiComm;
ParMetisMesh *parMetisMesh;
+
+ std::map<DegreeOfFreedom, DegreeOfFreedom> *mapLocalGlobal;
};
}