From c2ebc43542c8616213f901c1471cd131795b6933 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 19 Jun 2009 13:21:09 +0000
Subject: [PATCH] Solved bug, when coupled problems have different basis
functions degree.
---
AMDiS/libtool | 16 +--
AMDiS/src/Element.h | 8 +-
AMDiS/src/FiniteElemSpace.cc | 34 +++----
AMDiS/src/Lagrange.cc | 9 +-
AMDiS/src/MacroReader.cc | 75 ++++++--------
AMDiS/src/Mesh.cc | 20 ++--
AMDiS/src/ProblemVec.cc | 32 +++---
AMDiS/src/ResidualEstimator.cc | 172 ++++++++++++++++-----------------
AMDiS/src/ResidualEstimator.h | 26 ++---
AMDiS/src/Triangle.cc | 22 ++---
AMDiS/src/Triangle.h | 10 +-
11 files changed, 191 insertions(+), 233 deletions(-)
diff --git a/AMDiS/libtool b/AMDiS/libtool
index f7da7548..41ace2e6 100755
--- a/AMDiS/libtool
+++ b/AMDiS/libtool
@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
-# Libtool was configured on host NWRW03:
+# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -328,7 +328,7 @@ link_all_deplibs=unknown
sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/"
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-3.0.1 /usr/lib/qt-3.3/lib "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
@@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
-# Libtool was configured on host NWRW03:
+# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -7041,7 +7041,7 @@ link_all_deplibs=unknown
sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/"
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-3.0.1 /usr/lib/qt-3.3/lib "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
@@ -7065,7 +7065,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
-# Libtool was configured on host NWRW03:
+# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
@@ -7109,7 +7109,7 @@ LTCC="gcc"
LTCFLAGS="-g -O2"
# A language-specific compiler.
-CC="f95"
+CC="g77"
# Is the compiler the GNU C compiler?
with_gcc=yes
@@ -7346,10 +7346,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
-sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/"
+sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../ /lib/i386-redhat-linux/3.4.6/ /lib/ /usr/lib/i386-redhat-linux/3.4.6/ /usr/lib/"
# Run-time system search path for libraries
-sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-3.0.1 /usr/lib/qt-3.3/lib "
+sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h
index 305d9fa2..15a06b95 100644
--- a/AMDiS/src/Element.h
+++ b/AMDiS/src/Element.h
@@ -328,7 +328,6 @@ namespace AMDiS {
/** \} */
- // ===== pure virtual methods =================================================
/** \name pure virtual methods
* \{
@@ -339,12 +338,9 @@ namespace AMDiS {
* Used by Estimator for the jumps => same quadrature nodes from both sides!
*/
virtual const FixVec<int,WORLD>&
- sortFaceIndices(int face, FixVec<int,WORLD> *vec) const = 0;
+ sortFaceIndices(int face, FixVec<int, WORLD> *vec) const = 0;
- /** \brief
- * Returns a copy of itself. Needed by Mesh to create Elements by a
- * prototype.
- */
+ /// Returns a copy of itself. Needed by Mesh to create Elements by a prototype.
virtual Element *clone() = 0;
/** \brief
diff --git a/AMDiS/src/FiniteElemSpace.cc b/AMDiS/src/FiniteElemSpace.cc
index 8fd00e8f..b2582bf1 100644
--- a/AMDiS/src/FiniteElemSpace.cc
+++ b/AMDiS/src/FiniteElemSpace.cc
@@ -10,12 +10,12 @@ namespace AMDiS {
std::vector<FiniteElemSpace*> FiniteElemSpace::feSpaces;
FiniteElemSpace::FiniteElemSpace(DOFAdmin* admin_,
- const BasisFunction* bas_fcts_,
+ const BasisFunction* bas_fcts,
Mesh* aMesh,
const std::string& aString)
: name(aString),
admin(admin_),
- basFcts(bas_fcts_),
+ basFcts(bas_fcts),
mesh(aMesh)
{
FUNCNAME("FiniteElemSpace::FiniteElemSpace()");
@@ -30,17 +30,16 @@ namespace AMDiS {
for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
admin_local = &(mesh->getDOFAdmin(i));
int j = 0;
- for (; j <= mesh->getDim(); j++) {
+ for (; j <= mesh->getDim(); j++)
if (admin_local->getNumberOfDOFs(j) != (*ndof)[j])
break;
- }
if (j > mesh->getDim())
break;
admin_local = NULL;
}
- if (!admin_local) {
+ if (!admin_local)
admin_local = mesh->createDOFAdmin(name, *ndof);
- }
+
admin = const_cast<DOFAdmin*>(admin_local);
}
@@ -76,24 +75,13 @@ namespace AMDiS {
{
int numSpaces = static_cast<int>(feSpaces.size());
- for (int i = 0; i < numSpaces; i++) {
- if (admin) {
- if (feSpaces[i]->admin == admin) {
- return feSpaces[i];
- }
- } else {
- if (feSpaces[i]->basFcts == basFcts && feSpaces[i]->mesh == mesh) {
- return feSpaces[i];
- }
- }
- }
+ for (int i = 0; i < numSpaces; i++)
+ if (feSpaces[i]->basFcts == basFcts &&
+ feSpaces[i]->mesh == mesh &&
+ (!admin || (admin && feSpaces[i]->admin == admin)))
+ return feSpaces[i];
- if (admin) {
- ERROR_EXIT("no fespace found for this admin\n");
- return NULL;
- } else {
- return new FiniteElemSpace(NULL, basFcts, mesh, name_);
- }
+ return new FiniteElemSpace(admin, basFcts, mesh, name_);
}
int FiniteElemSpace::calcMemoryUsage()
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 5de5cba6..5a32a601 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -1046,11 +1046,14 @@ namespace AMDiS {
int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(CENTER);
Element *el = list->getElement(0);
int node = drv->getFESpace()->getMesh()->getNode(CENTER);
- DegreeOfFreedom dof0 = el->getDOF(node, n0); /* Parent center */
- DegreeOfFreedom dof_new = el->getChild(0)->getDOF(node, n0); /* newest vertex is center */
+ // Parent center
+ DegreeOfFreedom dof0 = el->getDOF(node, n0);
+ // Newest vertex is center
+ DegreeOfFreedom dof_new = el->getChild(0)->getDOF(node, n0);
(*drv)[dof_new] = (*drv)[dof0];
- dof_new = el->getChild(1)->getDOF(node, n0); /* newest vertex is center */
+ // Newest vertex is center
+ dof_new = el->getChild(1)->getDOF(node, n0);
(*drv)[dof_new] = (*drv)[dof0];
}
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index 5d826df2..41815863 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -158,9 +158,8 @@ namespace AMDiS {
associated = new VertexVector(mesh->getVertexAdmin(), "vertex vector");
mesh->periodicAssociations[boundaryType] = associated;
VertexVector::Iterator it(associated, ALL_DOFS);
- for (it.reset2(); !it.end(); ++it) {
+ for (it.reset2(); !it.end(); ++it)
*it = it.getDOFIndex();
- }
}
for (int j = 0; j < dim; j++) {
@@ -213,11 +212,10 @@ namespace AMDiS {
}
}
- for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
- if (periodicMap.getEntry(i) != -1) {
- MSG("identification : vertex %d is now vertex %d\n", i, periodicMap.getEntry(i));
- }
- }
+ for (int i = 0; i < mesh->getNumberOfVertices(); i++)
+ if (periodicMap.getEntry(i) != -1)
+ MSG("identification : vertex %d is now vertex %d\n",
+ i, periodicMap.getEntry(i));
}
// =========================================================
@@ -268,7 +266,6 @@ namespace AMDiS {
boundaryDOFs(mesh);
// initial boundary projections
- //if(dim > 1) {
int numFaces = mesh->getGeo(FACE);
int dim = mesh->getDim();
mel = mesh->firstMacroElement();
@@ -291,7 +288,6 @@ namespace AMDiS {
}
}
}
- //}
macroInfo->fillBoundaryInfo(mesh);
@@ -1409,46 +1405,35 @@ namespace AMDiS {
&(nei->projection[mesh->getGeo(FACE)+edge_no]);
if (mesh->getNumberOfDOFs(EDGE))
- {
- // if(periodic) {
- // nei->element->setDOF(node+edge_no, mesh->getDOF(EDGE));
- // } else {
- nei->element->setDOF(node+edge_no,edge_dof);
- // }
- }
+ nei->element->setDOF(node+edge_no,edge_dof);
- if (next_el[edge_no][0] != opp_v)
- {
- if(nei->getBoundary(next_el[edge_no][0])) {
- lbound = newBound(nei->getBoundary(next_el[edge_no][0]), lbound);
- Projection *neiProject = nei->getProjection(next_el[edge_no][0]);
- if(!lproject)
+ if (next_el[edge_no][0] != opp_v) {
+ if(nei->getBoundary(next_el[edge_no][0])) {
+ lbound = newBound(nei->getBoundary(next_el[edge_no][0]), lbound);
+ Projection *neiProject = nei->getProjection(next_el[edge_no][0]);
+ if (!lproject) {
+ lproject = neiProject;
+ } else {
+ if (neiProject && (lproject->getID() < neiProject->getID()))
lproject = neiProject;
- else {
- if(neiProject && (lproject->getID() < neiProject->getID())) {
- lproject = neiProject;
- }
- }
}
- opp_v = nei->getOppVertex(next_el[edge_no][0]);
- nei = nei->getNeighbour(next_el[edge_no][0]);
}
- else
- {
- if(nei->getBoundary(next_el[edge_no][1])) {
- lbound = newBound(nei->getBoundary(next_el[edge_no][1]), lbound);
- Projection *neiProject = nei->getProjection(next_el[edge_no][1]);
- if(!lproject)
+ opp_v = nei->getOppVertex(next_el[edge_no][0]);
+ nei = nei->getNeighbour(next_el[edge_no][0]);
+ } else {
+ if (nei->getBoundary(next_el[edge_no][1])) {
+ lbound = newBound(nei->getBoundary(next_el[edge_no][1]), lbound);
+ Projection *neiProject = nei->getProjection(next_el[edge_no][1]);
+ if (!lproject) {
+ lproject = neiProject;
+ } else {
+ if (neiProject && (lproject->getID() < neiProject->getID()))
lproject = neiProject;
- else {
- if(neiProject && (lproject->getID() < neiProject->getID())) {
- lproject = neiProject;
- }
- }
}
- opp_v = nei->getOppVertex(next_el[edge_no][1]);
- nei = nei->getNeighbour(next_el[edge_no][1]);
}
+ opp_v = nei->getOppVertex(next_el[edge_no][1]);
+ nei = nei->getNeighbour(next_el[edge_no][1]);
+ }
}
if (!nei)
{
@@ -1861,10 +1846,10 @@ namespace AMDiS {
{
FUNCNAME("MacroReader::checkMesh()");
- int i, nused, nfree;
+ int i, nused, nfree;
DOFAdmin *localAdmin=mesh->admin[0];
- Flag fill_flag;
- int error_detected = 0;
+ Flag fill_flag;
+ int error_detected = 0;
MSG("checking mesh\n");
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 57c15af7..98b3535e 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -387,14 +387,15 @@ namespace AMDiS {
localAdmin->setMesh(this);
- std::vector<DOFAdmin*>::iterator dai = std::find(admin.begin(),admin.end(),localAdmin);
+ std::vector<DOFAdmin*>::iterator dai =
+ std::find(admin.begin(),admin.end(),localAdmin);
TEST_EXIT(dai == admin.end())
("admin %s is already associated to mesh %s\n",
localAdmin->getName().c_str(), this->getName().c_str());
// if this will be required, see the untested code in revision < 224
- TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work yet!\n");
+ // TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work yet!\n");
admin.push_back(localAdmin);
@@ -435,9 +436,8 @@ namespace AMDiS {
}
node[CENTER] = nNodeEl;
- if (nDOF[CENTER] > 0) {
+ if (nDOF[CENTER] > 0)
nNodeEl += 1;
- }
}
void Mesh::dofCompress()
@@ -1082,25 +1082,25 @@ namespace AMDiS {
// If there is no value file which should be written, we can delete
// the information of the macro file.
- if (!valueFilename.length()) {
+ if (!valueFilename.length())
clearMacroFileInfo();
- }
}
initialized = true;
}
- bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
+ bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2)
+ {
std::map<BoundaryType, VertexVector*>::iterator it;
std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
- for (it = periodicAssociations.begin(); it != end; ++it) {
+ for (it = periodicAssociations.begin(); it != end; ++it)
if ((*(it->second))[dof1] == dof2)
return true;
- }
return false;
}
- bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
+ bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2)
+ {
std::vector<DegreeOfFreedom> associatedToDOF1;
std::map<BoundaryType, VertexVector*>::iterator it;
std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 4fe4075f..e832ec49 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -72,17 +72,9 @@ namespace AMDiS {
WARNING("feSpaces already created\n");
} else {
if (initFlag.isSet(INIT_FE_SPACE) ||
- (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
+ (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE)))
+ createFESpace(NULL);
- if (!adoptFlag.isSet(INIT_MESH)) {
- createFESpace(NULL);
- } else {
- TEST_EXIT(meshes.size() == 1)("I have to think about it!\n");
- TEST_EXIT(meshes[0]->getNumberOfDOFAdmin() == 1)("I have to think about it!\n");
-
- createFESpace(&(const_cast<DOFAdmin&>(meshes[0]->getDOFAdmin(0))));
- }
- }
if (adoptProblem &&
(adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
feSpaces = adoptProblem->getFESpaces();
@@ -106,9 +98,9 @@ namespace AMDiS {
WARNING("no feSpace created\n");
// === create system ===
- if (initFlag.isSet(INIT_SYSTEM)) {
+ if (initFlag.isSet(INIT_SYSTEM))
createMatricesAndVectors();
- }
+
if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
solution = adoptProblem->getSolution();
rhs = adoptProblem->getRHS();
@@ -189,9 +181,13 @@ namespace AMDiS {
&globalRefinements);
// Initialize the meshes if there is no serialization file.
- for (int i = 0; i < static_cast<int>(meshes.size()); i++)
- if (initFlag.isSet(INIT_MESH) && meshes[i] && !(meshes[i]->isInitialized()))
+ for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
+ bool initMesh = initFlag.isSet(INIT_MESH) ||
+ (adoptProblem && adoptFlag.isSet(INIT_MESH));
+
+ if (initMesh && meshes[i] && !(meshes[i]->isInitialized()))
meshes[i]->initialize();
+ }
// === read value file and use it for the mesh values ===
std::string valueFilename("");
@@ -233,9 +229,9 @@ namespace AMDiS {
sprintf(number, "%d", i);
int refSet = -1;
GET_PARAMETER(0, name + "->refinement set[" + number + "]", "%d", &refSet);
- if (refSet < 0) {
+ if (refSet < 0)
refSet = 0;
- }
+
if (meshForRefinementSet[refSet] == NULL) {
Mesh *newMesh = new Mesh(meshName, dim);
meshForRefinementSet[refSet] = newMesh;
@@ -296,9 +292,9 @@ namespace AMDiS {
}
for (int i = 0; i < nComponents; i++) {
- for (int j = 0; j < nComponents; j++) {
+ for (int j = 0; j < nComponents; j++)
traverseInfo.getMatrix(i, j).setFESpace(componentSpaces[i], componentSpaces[j]);
- }
+
traverseInfo.getVector(i).setFESpace(componentSpaces[i]);
}
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index f01ac356..fab865ff 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -12,9 +12,11 @@ namespace AMDiS {
: Estimator(name, r),
C0(1.0),
C1(1.0),
- C2(1.0),
+ C2(0.0),
C3(1.0)
{
+ FUNCNAME("ResidualEstimator::ResidualEstimator()");
+
GET_PARAMETER(0, name + "->C0", "%f", &C0);
GET_PARAMETER(0, name + "->C1", "%f", &C1);
GET_PARAMETER(0, name + "->C2", "%f", &C2);
@@ -24,6 +26,8 @@ namespace AMDiS {
C1 = C1 > 1.e-25 ? sqr(C1) : 0.0;
C2 = C2 > 1.e-25 ? sqr(C2) : 0.0;
C3 = C3 > 1.e-25 ? sqr(C3) : 0.0;
+
+ TEST_EXIT(C2 == 0.0)("C2 is not used! Please remove it or set it to 0.0!\n");
}
void ResidualEstimator::init(double ts)
@@ -46,22 +50,19 @@ namespace AMDiS {
basFcts[system] = uh[system]->getFESpace()->getBasisFcts();
degree = std::max(degree, basFcts[system]->getDegree());
}
-
degree *= 2;
quad = Quadrature::provideQuadrature(dim, degree);
nPoints = quad->getNumPoints();
Flag flag = INIT_PHI | INIT_GRD_PHI;
- if (degree > 2) {
- flag |= INIT_D2_PHI;
- }
+ if (degree > 2)
+ flag |= INIT_D2_PHI;
- for (int system = 0; system < nSystems; system++) {
+ for (int system = 0; system < nSystems; system++)
quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system],
*quad,
- flag);
- }
+ flag);
uhEl = new double*[nSystems];
uhNeigh = new double*[nSystems];
@@ -76,9 +77,7 @@ namespace AMDiS {
uhQP = timestep ? new double[nPoints] : NULL;
uhOldQP = timestep ? new double[nPoints] : NULL;
-
riq = new double[nPoints];
-
grdUh_qp = NULL;
D2uhqp = NULL;
@@ -110,16 +109,16 @@ namespace AMDiS {
neighInfo = mesh->createNewElInfo();
// prepare date for computing jump residual
- if (C1 && (dim > 1)) {
- surfaceQuad_ = Quadrature::provideQuadrature(dim - 1, degree);
- nPointsSurface_ = surfaceQuad_->getNumPoints();
- grdUhEl_.resize(nPointsSurface_);
- grdUhNeigh_.resize(nPointsSurface_);
- jump_.resize(nPointsSurface_);
- localJump_.resize(nPointsSurface_);
- neighbours_ = Global::getGeo(NEIGH, dim);
- lambdaNeigh_ = new DimVec<WorldVector<double> >(dim, NO_INIT);
- lambda_ = new DimVec<double>(dim, NO_INIT);
+ if (C1 > 0.0 && dim > 1) {
+ surfaceQuad = Quadrature::provideQuadrature(dim - 1, degree);
+ nPointsSurface = surfaceQuad->getNumPoints();
+ grdUhEl.resize(nPointsSurface);
+ grdUhNeigh.resize(nPointsSurface);
+ jump.resize(nPointsSurface);
+ localJump.resize(nPointsSurface);
+ nNeighbours = Global::getGeo(NEIGH, dim);
+ lambdaNeigh = new DimVec<WorldVector<double> >(dim, NO_INIT);
+ lambda = new DimVec<double>(dim, NO_INIT);
}
}
@@ -165,8 +164,8 @@ namespace AMDiS {
delete [] D2uhqp;
if (C1 && (dim > 1)) {
- delete lambdaNeigh_;
- delete lambda_;
+ delete lambdaNeigh;
+ delete lambda;
}
delete neighInfo;
@@ -187,43 +186,40 @@ namespace AMDiS {
double est_el = el->getEstimation(row);
double h2 = h2_from_det(det, dim);
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
riq[iq] = 0.0;
- }
for (int system = 0; system < nSystems; system++) {
if (matrix[system] == NULL)
continue;
- // init assemblers
+ // === init assemblers ===
+
for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
- ++it, ++itfac) {
- if (*itfac == NULL || **itfac != 0.0) {
+ ++it, ++itfac)
+ if (*itfac == NULL || **itfac != 0.0)
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
- }
- }
- if (C0) {
+ if (C0 > 0.0)
for (it = const_cast<DOFVector<double>*>(fh[system])->getOperatorsBegin();
it != const_cast<DOFVector<double>*>(fh[system])->getOperatorsEnd();
- ++it) {
+ ++it)
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
- }
- }
if (timestep && uhOld[system]) {
TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
uhOld[system]->getLocalVector(el, uhOldEl[system]);
- // ===== time and element residuals
- if (C0 || C3) {
+ // === Compute time error. ===
+
+ if (C0 > 0.0 || C3 > 0.0) {
uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
- if (C3 && uhOldQP && system == std::max(row, 0)) {
+ if (C3 > 0.0 && uhOldQP && system == std::max(row, 0)) {
val = 0.0;
for (int iq = 0; iq < nPoints; iq++) {
double tiq = (uhQP[iq] - uhOldQP[iq]);
@@ -236,27 +232,32 @@ namespace AMDiS {
}
}
- if (C0) {
+ // === Compute element residual. ===
+
+ if (C0 > 0.0) {
for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
++it, ++itfac) {
if (*itfac == NULL || **itfac != 0.0) {
- if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
+ if (uhQP == NULL && (*it)->zeroOrderTerms()) {
uhQP = new double[nPoints];
uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
}
- if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
+ if (grdUh_qp == NULL &&
+ ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
grdUh_qp = new WorldVector<double>[nPoints];
uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
}
- if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) {
+ if (D2uhqp == NULL && degree > 2 && (*it)->secondOrderTerms()) {
D2uhqp = new WorldMatrix<double>[nPoints];
uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);
}
}
}
+ // === Compute the element residual and store it in irq. ===
+
r(elInfo,
nPoints,
uhQP,
@@ -284,17 +285,19 @@ namespace AMDiS {
est_el += val;
- // ===== jump residuals
+
+ // === Compute jump residuals. ===
+
if (C1 && (dim > 1)) {
int dow = Global::getGeo(WORLD);
- for (int face = 0; face < neighbours_; face++) {
+ for (int face = 0; face < nNeighbours; face++) {
Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
if (neigh && neigh->getMark()) {
int oppV = elInfo->getOppVertex(face);
- el->sortFaceIndices(face, &faceIndEl_);
- neigh->sortFaceIndices(oppV, &faceIndNeigh_);
+ el->sortFaceIndices(face, &faceIndEl);
+ neigh->sortFaceIndices(oppV, &faceIndNeigh);
neighInfo->setElement(const_cast<Element*>(neigh));
neighInfo->setFillFlag(Mesh::FILL_COORDS);
@@ -315,8 +318,8 @@ namespace AMDiS {
for (it = infoList.begin(); it != infoList.end(); ++it) {
if (it->elementSide == face) {
for (int i = 0; i < dim; i++) {
- int i1 = faceIndEl_[i];
- int i2 = faceIndNeigh_[i];
+ int i1 = faceIndEl[i];
+ int i2 = faceIndNeigh[i];
int j = 0;
for (; j < dim; j++) {
@@ -340,24 +343,21 @@ namespace AMDiS {
if (!periodicCoords) {
for (int i = 0; i < dim; i++) {
- int i1 = faceIndEl_[i];
- int i2 = faceIndNeigh_[i];
+ int i1 = faceIndEl[i];
+ int i2 = faceIndNeigh[i];
for (int j = 0; j < dow; j++)
neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
}
}
Parametric *parametric = mesh->getParametric();
- if (parametric) {
- neighInfo = parametric->addParametricInfo(neighInfo);
- }
+ if (parametric)
+ neighInfo = parametric->addParametricInfo(neighInfo);
- double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh_));
+ double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh));
- for (int iq = 0; iq < nPointsSurface_; iq++) {
- jump_[iq].set(0.0);
- }
-
+ for (int iq = 0; iq < nPointsSurface; iq++)
+ jump[iq].set(0.0);
for (int system = 0; system < nSystems; system++) {
if (matrix[system] == NULL)
@@ -366,28 +366,26 @@ namespace AMDiS {
uh[system]->getLocalVector(el, uhEl[system]);
uh[system]->getLocalVector(neigh, uhNeigh[system]);
- for (int iq = 0; iq < nPointsSurface_; iq++) {
- (*lambda_)[face] = 0.0;
- for (int i = 0; i < dim; i++) {
- (*lambda_)[faceIndEl_[i]] = surfaceQuad_->getLambda(iq, i);
- }
+ for (int iq = 0; iq < nPointsSurface; iq++) {
+ (*lambda)[face] = 0.0;
+ for (int i = 0; i < dim; i++)
+ (*lambda)[faceIndEl[i]] = surfaceQuad->getLambda(iq, i);
- basFcts[system]->evalGrdUh(*lambda_,
+ basFcts[system]->evalGrdUh(*lambda,
grdLambda,
uhEl[system],
- &grdUhEl_[iq]);
+ &grdUhEl[iq]);
- (*lambda_)[oppV] = 0.0;
- for (int i = 0; i < dim; i++) {
- (*lambda_)[faceIndNeigh_[i]] = surfaceQuad_->getLambda(iq, i);
- }
+ (*lambda)[oppV] = 0.0;
+ for (int i = 0; i < dim; i++)
+ (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
- basFcts[system]->evalGrdUh(*lambda_,
- *lambdaNeigh_,
+ basFcts[system]->evalGrdUh(*lambda,
+ *lambdaNeigh,
uhNeigh[system],
- &grdUhNeigh_[iq]);
+ &grdUhNeigh[iq]);
- grdUhEl_[iq] -= grdUhNeigh_[iq];
+ grdUhEl[iq] -= grdUhNeigh[iq];
}
std::vector<double*>::iterator fac;
@@ -398,31 +396,26 @@ namespace AMDiS {
++it, ++fac) {
if (*fac == NULL || **fac != 0.0) {
- for (int iq = 0; iq < nPointsSurface_; iq++) {
- localJump_[iq].set(0.0);
- }
+ for (int iq = 0; iq < nPointsSurface; iq++)
+ localJump[iq].set(0.0);
- (*it)->weakEvalSecondOrder(nPointsSurface_,
- grdUhEl_.getValArray(),
- localJump_.getValArray());
+ (*it)->weakEvalSecondOrder(nPointsSurface,
+ grdUhEl.getValArray(),
+ localJump.getValArray());
double factor = *fac ? **fac : 1.0;
- if (factor != 1.0) {
- for (int i = 0; i < nPointsSurface_; i++) {
- localJump_[i] *= factor;
- }
- }
+ if (factor != 1.0)
+ for (int i = 0; i < nPointsSurface; i++)
+ localJump[i] *= factor;
- for (int i = 0; i < nPointsSurface_; i++) {
- jump_[i] += localJump_[i];
- }
+ for (int i = 0; i < nPointsSurface; i++)
+ jump[i] += localJump[i];
}
}
}
val = 0.0;
- for (int iq = 0; iq < nPointsSurface_; iq++) {
- val += surfaceQuad_->getWeight(iq) * (jump_[iq] * jump_[iq]);
- }
+ for (int iq = 0; iq < nPointsSurface; iq++)
+ val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
double d = 0.5 * (det + detNeigh);
@@ -431,9 +424,8 @@ namespace AMDiS {
else
val *= C1 * d;
- if (parametric) {
+ if (parametric)
neighInfo = parametric->removeParametricInfo(neighInfo);
- }
neigh->setEstimation(neigh->getEstimation(row) + val, row);
est_el += val;
diff --git a/AMDiS/src/ResidualEstimator.h b/AMDiS/src/ResidualEstimator.h
index 97ddbc1c..a2ac7a85 100644
--- a/AMDiS/src/ResidualEstimator.h
+++ b/AMDiS/src/ResidualEstimator.h
@@ -94,7 +94,7 @@ namespace AMDiS {
/// Constant in front of edge/face residual
double C1;
- /// Constant in front of coarsening term
+ /// Not used! Was thought to be the constant in front of coarsening term.
double C2;
/// Constant in front of the time
@@ -126,6 +126,7 @@ namespace AMDiS {
double *uhOldQP;
+ /// Stores the element residual computed at the quadrature points of the element.
double *riq;
WorldVector<double> *grdUh_qp;
@@ -134,27 +135,28 @@ namespace AMDiS {
ElInfo *neighInfo;
- Quadrature *surfaceQuad_;
+ Quadrature *surfaceQuad;
- int nPointsSurface_;
+ int nPointsSurface;
- Vector<WorldVector<double> > grdUhEl_;
+ Vector<WorldVector<double> > grdUhEl;
- Vector<WorldVector<double> > grdUhNeigh_;
+ Vector<WorldVector<double> > grdUhNeigh;
- Vector<WorldVector<double> > jump_;
+ Vector<WorldVector<double> > jump;
- Vector<WorldVector<double> > localJump_;
+ Vector<WorldVector<double> > localJump;
- WorldVector<int> faceIndEl_;
+ WorldVector<int> faceIndEl;
- WorldVector<int> faceIndNeigh_;
+ WorldVector<int> faceIndNeigh;
- DimVec<WorldVector<double> > *lambdaNeigh_;
+ DimVec<WorldVector<double> > *lambdaNeigh;
- DimVec<double> *lambda_;
+ DimVec<double> *lambda;
- int neighbours_;
+ /// Maximal number of neighbours an element may have in the used dimension.
+ int nNeighbours;
};
}
diff --git a/AMDiS/src/Triangle.cc b/AMDiS/src/Triangle.cc
index 3dadb368..1ef94d5d 100644
--- a/AMDiS/src/Triangle.cc
+++ b/AMDiS/src/Triangle.cc
@@ -49,15 +49,15 @@ namespace AMDiS {
FixVec<int,WORLD> *val = NULL;
vof = vertexOfEdge[face];
- if (NULL==sorted_2d) {
- sorted_2d = new MatrixOfFixVecs<FixVec<int,WORLD> >(2,3,2,NO_INIT);
+ if (NULL == sorted_2d) {
+ sorted_2d = new MatrixOfFixVecs<FixVec<int,WORLD> >(2, 3, 2, NO_INIT);
- (*sorted_2d)[1][0][1]=(*sorted_2d)[1][1][0]=
- (*sorted_2d)[2][0][0]=(*sorted_2d)[2][1][1]=0;
- (*sorted_2d)[0][0][0]=(*sorted_2d)[0][1][1]=
- (*sorted_2d)[2][0][1]=(*sorted_2d)[2][1][0]=1;
- (*sorted_2d)[0][0][1]=(*sorted_2d)[0][1][0]=
- (*sorted_2d)[1][0][0]=(*sorted_2d)[1][1][1]=2;
+ (*sorted_2d)[1][0][1] = (*sorted_2d)[1][1][0] =
+ (*sorted_2d)[2][0][0] = (*sorted_2d)[2][1][1] = 0;
+ (*sorted_2d)[0][0][0] = (*sorted_2d)[0][1][1] =
+ (*sorted_2d)[2][0][1] = (*sorted_2d)[2][1][0] = 1;
+ (*sorted_2d)[0][0][1] = (*sorted_2d)[0][1][0] =
+ (*sorted_2d)[1][0][0] = (*sorted_2d)[1][1][1] = 2;
}
if (dof[vof[0]][0] < dof[vof[1]][0])
@@ -68,11 +68,11 @@ namespace AMDiS {
if (vec) {
val = vec;
*val = (*sorted_2d)[face][no];
- }
- else
+ } else {
val = &((*sorted_2d)[face][no]);
+ }
- return(*(const_cast<const FixVec<int,WORLD>* >(val)));
+ return *(const_cast<const FixVec<int,WORLD>* >(val));
}
}
diff --git a/AMDiS/src/Triangle.h b/AMDiS/src/Triangle.h
index 4a279cc4..fa8b2f64 100644
--- a/AMDiS/src/Triangle.h
+++ b/AMDiS/src/Triangle.h
@@ -94,8 +94,7 @@ namespace AMDiS {
bool hasSide(Element* sideElem) const;
/// implements Element::sortFaceIndices
- const FixVec<int,WORLD>& sortFaceIndices(int face,
- FixVec<int,WORLD> *vec) const;
+ const FixVec<int,WORLD>& sortFaceIndices(int face, FixVec<int,WORLD> *vec) const;
/// implements Element::isLine. Returns false because this element is a Triangle
inline bool isLine() const
@@ -109,10 +108,7 @@ namespace AMDiS {
return true;
}
- /** \brief
- * implements Element::isTetrahedron. Returns false because this element is a
- * Triangle
- */
+ /// implements Element::isTetrahedron. Returns false because this element is a Triangle
inline bool isTetrahedron() const
{
return false;
@@ -148,7 +144,7 @@ namespace AMDiS {
TEST_EXIT_DBG(edge >= 0 && edge < 3)("invalid edge\n");
return edge;
}
-
+
std::string getTypeName() const
{
return "Triangle";
--
GitLab