From 1ef3e6b45f30a9c9fb94d9ba486b6eaae1714f38 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 27 Jun 2008 06:47:31 +0000
Subject: [PATCH] * Residual estimator around 50% faster!
---
AMDiS/src/ElInfo3d.cc | 19 +++----
AMDiS/src/ResidualEstimator.cc | 91 +++++++++++++++++-----------------
AMDiS/src/ResidualEstimator.h | 18 +++++++
3 files changed, 69 insertions(+), 59 deletions(-)
diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc
index 9114c498..2e284a2e 100644
--- a/AMDiS/src/ElInfo3d.cc
+++ b/AMDiS/src/ElInfo3d.cc
@@ -37,10 +37,6 @@ namespace AMDiS {
}
}
- // if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
- // det = calcGrdLambda(*Lambda);
- // }
-
int neighbours = mesh_->getGeo(NEIGH);
if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS) ||
@@ -81,9 +77,6 @@ namespace AMDiS {
}
}
- //int faces = mesh->getGeo(FACE);
- //int edges = mesh->getGeo(EDGE);
-
if (fillFlag_.isSet(Mesh::FILL_BOUND)) {
for (int i = 0; i < element_->getGeo(BOUNDARY); i++) {
boundary_[i] = mel->getBoundary(i);
@@ -583,14 +576,15 @@ namespace AMDiS {
boundary_[2] = elinfo_old->getBoundary(cv[2]);
boundary_[3] = elinfo_old->getBoundary(ochild);
+ int geoFace = mesh_->getGeo(FACE);
+
ce = const_cast<int*>(Tetrahedron::childEdge[el_type_local][ichild]);
for (int iedge = 0; iedge < 4; iedge++) {
- boundary_[mesh_->getGeo(FACE) + iedge] =
- elinfo_old->getBoundary(mesh_->getGeo(FACE)+ce[iedge]);
+ boundary_[geoFace + iedge] = elinfo_old->getBoundary(geoFace + ce[iedge]);
}
for (int iedge = 4; iedge < 6; iedge++) {
int i = 5 - cv[iedge - 3]; /* old vertex opposite new edge */
- boundary_[mesh_->getGeo(FACE) + iedge] = elinfo_old->getBoundary(i);
+ boundary_[geoFace + iedge] = elinfo_old->getBoundary(i);
}
if (elinfo_old->getProjection(0) &&
@@ -604,12 +598,11 @@ namespace AMDiS {
projection_[3] = elinfo_old->getProjection(ochild);
for (int iedge = 0; iedge < 4; iedge++) {
- projection_[mesh_->getGeo(FACE)+iedge] =
- elinfo_old->getProjection(mesh_->getGeo(FACE)+ce[iedge]);
+ projection_[geoFace + iedge] = elinfo_old->getProjection(mesh_->getGeo(FACE)+ce[iedge]);
}
for (int iedge = 4; iedge < 6; iedge++) {
int i = 5 - cv[iedge - 3]; /* old vertex opposite new edge */
- projection_[mesh_->getGeo(FACE) + iedge] = elinfo_old->getProjection(i);
+ projection_[geoFace + iedge] = elinfo_old->getProjection(i);
}
}
}
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index aeb70274..a0bb11f2 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -64,10 +64,12 @@ namespace AMDiS {
}
uhEl = GET_MEMORY(double*, numSystems);
+ uhNeigh = GET_MEMORY(double*, numSystems);
uhOldEl = timestep ? GET_MEMORY(double*, numSystems) : NULL;
for (int system = 0; system < numSystems; system++) {
uhEl[system] = GET_MEMORY(double, basFcts[system]->getNumber());
+ uhNeigh[system] = GET_MEMORY(double, basFcts[system]->getNumber());
if (timestep)
uhOldEl[system] = GET_MEMORY(double, basFcts[system]->getNumber());
}
@@ -104,6 +106,19 @@ namespace AMDiS {
Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_DET |
Mesh::CALL_LEAF_EL;
+
+ 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);
+ }
}
void ResidualEstimator::exit(bool output)
@@ -115,11 +130,13 @@ namespace AMDiS {
for (int system = 0; system < numSystems; system++) {
FREE_MEMORY(uhEl[system], double, basFcts[system]->getNumber());
+ FREE_MEMORY(uhNeigh[system], double, basFcts[system]->getNumber());
if (timestep)
FREE_MEMORY(uhOldEl[system], double, basFcts[system]->getNumber());
}
FREE_MEMORY(uhEl, double*, numSystems);
+ FREE_MEMORY(uhNeigh, double*, numSystems);
if (timestep) {
FREE_MEMORY(uhOldEl, double*, numSystems);
@@ -148,6 +165,8 @@ namespace AMDiS {
if (D2uhqp != NULL) {
FREE_MEMORY(D2uhqp, WorldMatrix<double>, numPoints);
}
+
+ DELETE neighInfo;
}
void ResidualEstimator::estimateElement(ElInfo *elInfo)
@@ -265,17 +284,11 @@ namespace AMDiS {
// ===== jump residuals
if (C1 && (dim > 1)) {
- int neighbours = Global::getGeo(NEIGH, dim);
- int face;
- Quadrature *surfaceQuad = Quadrature::provideQuadrature(dim - 1, degree);
- int numPointsSurface = surfaceQuad->getNumPoints();
- Vector<WorldVector<double> > jump(numPointsSurface);
-
+ int dow = Global::getGeo(WORLD);
- for (face = 0; face < neighbours; face++) {
+ for (int face = 0; face < neighbours_; face++) {
neigh = const_cast<Element*>(elInfo->getNeighbour(face));
if (neigh && neigh->getMark()) {
- ElInfo *neighInfo = mesh->createNewElInfo();
WorldVector<int> faceIndEl, faceIndNeigh;
int oppV = elInfo->getOppVertex(face);
DimVec<WorldVector<double> > LambdaNeigh(dim, NO_INIT);
@@ -287,9 +300,7 @@ namespace AMDiS {
neighInfo->setElement(const_cast<Element*>(neigh));
neighInfo->setFillFlag(Mesh::FILL_COORDS);
-
- int dow = Global::getGeo(WORLD);
-
+
int j, i1, i2;
for (int i = 0; i < dow; i++)
@@ -346,11 +357,8 @@ namespace AMDiS {
detNeigh = abs(neighInfo->calcGrdLambda(LambdaNeigh));
- Vector<WorldVector<double> > grdUhEl(numPointsSurface);
- Vector<WorldVector<double> > grdUhNeigh(numPointsSurface);
-
- for (iq = 0; iq < numPointsSurface; iq++) {
- jump[iq].set(0.0);
+ for (iq = 0; iq < nPointsSurface_; iq++) {
+ jump_[iq].set(0.0);
}
@@ -358,69 +366,62 @@ namespace AMDiS {
if (matrix[system] == NULL)
continue;
- uh[system]->getLocalVector(el, uhEl[system]);
-
- double* uhNeigh;
- uhNeigh = new double[basFcts[system]->getNumber()];
- uh[system]->getLocalVector(neigh, uhNeigh);
-
+ uh[system]->getLocalVector(el, uhEl[system]);
+ uh[system]->getLocalVector(neigh, uhNeigh[system]);
- for (iq = 0; iq < numPointsSurface; iq++) {
+ for (iq = 0; iq < nPointsSurface_; iq++) {
lambda[face] = 0.0;
for (int i = 0; i < dim; i++) {
- lambda[faceIndEl[i]] = surfaceQuad->getLambda(iq, i);
+ lambda[faceIndEl[i]] = surfaceQuad_->getLambda(iq, i);
}
basFcts[system]->evalGrdUh(lambda,
Lambda,
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[faceIndNeigh[i]] = surfaceQuad_->getLambda(iq, i);
}
basFcts[system]->evalGrdUh(lambda,
LambdaNeigh,
- uhNeigh,
- &grdUhNeigh[iq]);
+ uhNeigh[system],
+ &grdUhNeigh_[iq]);
- grdUhEl[iq] -= grdUhNeigh[iq];
+ grdUhEl_[iq] -= grdUhNeigh_[iq];
}
- delete [] uhNeigh;
-
::std::vector<double*>::iterator fac;
for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
++it, ++fac) {
- Vector<WorldVector<double> > localJump(numPointsSurface);
- for (iq = 0; iq < numPointsSurface; iq++) {
- localJump[iq].set(0.0);
+ for (iq = 0; iq < nPointsSurface_; iq++) {
+ localJump_[iq].set(0.0);
}
- (*it)->weakEvalSecondOrder(numPointsSurface,
- 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 < numPointsSurface; i++) {
- localJump[i] *= factor;
+ for (int i = 0; i < nPointsSurface_; i++) {
+ localJump_[i] *= factor;
}
}
- for (int i = 0; i < numPointsSurface; i++) {
- jump[i] += localJump[i];
+ for (int i = 0; i < nPointsSurface_; i++) {
+ jump_[i] += localJump_[i];
}
}
}
-
- for (val = iq = 0; iq < numPointsSurface; iq++) {
- val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
+ val = 0.0;
+ for (iq = 0; iq < nPointsSurface_; iq++) {
+ val += surfaceQuad_->getWeight(iq) * (jump_[iq] * jump_[iq]);
}
double d = 0.5 * (det + detNeigh);
@@ -434,8 +435,6 @@ namespace AMDiS {
neighInfo = parametric->removeParametricInfo(neighInfo);
}
- DELETE neighInfo;
-
neigh->setEstimation(neigh->getEstimation(row) + val, row);
est_el += val;
diff --git a/AMDiS/src/ResidualEstimator.h b/AMDiS/src/ResidualEstimator.h
index ccf394b4..b111703f 100644
--- a/AMDiS/src/ResidualEstimator.h
+++ b/AMDiS/src/ResidualEstimator.h
@@ -143,6 +143,8 @@ namespace AMDiS {
double **uhOldEl;
+ double **uhNeigh;
+
double *uhQP;
double *uhOldQP;
@@ -152,6 +154,22 @@ namespace AMDiS {
WorldVector<double> *grdUh_qp;
WorldMatrix<double> *D2uhqp;
+
+ ElInfo *neighInfo;
+
+ Quadrature *surfaceQuad_;
+
+ int nPointsSurface_;
+
+ Vector<WorldVector<double> > grdUhEl_;
+
+ Vector<WorldVector<double> > grdUhNeigh_;
+
+ Vector<WorldVector<double> > jump_;
+
+ Vector<WorldVector<double> > localJump_;
+
+ int neighbours_;
};
}
--
GitLab