Skip to content
GitLab
Commit 1ef3e6b4 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Residual estimator around 50% faster!

parent ffc67701
Hide whitespace changes
Inline Side-by-side
AMDiS/src/ElInfo3d.cc
View file @ 1ef3e6b4
......@@ -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);
}
}
}
......
AMDiS/src/ResidualEstimator.cc
View file @ 1ef3e6b4
......@@ -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;
......
AMDiS/src/ResidualEstimator.h
View file @ 1ef3e6b4
......@@ -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_;
};
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment