/*
$Id: adapt.cc,v 1.4 1996/10/11 15:15:13 bzferdma Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "adapt.h"
#include "utils.h"
#include "numerics.h"
#include "problem.h"
#include "problemtr.h"
#include "elements.h"
#include "triang.h"
#include "triangA.h"
#include "int.h"
#include "cmdpars.h"
extern CmdPars Cmd;
//-------------------------------------------------------------------------
ErrorEstimator:: ErrorEstimator()
{
elementError = 0;
minRefinementRatio = 0.05; Cmd.get("minRefRatio", &minRefinementRatio);
infoErrorEstimator = 0; Cmd.get("infoErrorEstimator", &infoErrorEstimator);
timeErrorEstimator = 0; Cmd.get("timeErrorEstimator", &timeErrorEstimator);
accTime = 0; Cmd.get("accTimeErrorEstimator", &accTime);
refStrategy = uniform;
if (Cmd.isSet("refStrategy","maxValue")) refStrategy=maxValue;
else if (Cmd.isSet("refStrategy","extraPol")) refStrategy=extrapolation;
else if (Cmd.isSet("refStrategy","uniform")) refStrategy=uniform;
else if (Cmd.isSet("refStrategy","random")) refStrategy=random;
else MissingParameter("refStrategy");
}
//-------------------------------------------------------------------------
ErrorEstimator:: ~ErrorEstimator() { }
//-------------------------------------------------------------------------
Bool ErrorEstimator:: estimateError(Problem& /*problem*/)
{
if (infoErrorEstimator) cout << "\n No Error Estimator";
EnergyError = machMax(Real(0.0));
return False;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
DLY:: DLY(Element* elementDLY0, Interface* interfaceDLY0)
: ErrorEstimator(), elementDLY(elementDLY0), interfaceDLY(interfaceDLY0)
{ }
//-------------------------------------------------------------------------
DLY:: ~DLY()
{
delete elementDLY;
delete interfaceDLY;
}
//-------------------------------------------------------------------------
EFDLY:: EFDLY(Element* elementDLY0, Interface* interface0)
: DLY(elementDLY0, interface0)
{ }
//-------------------------------------------------------------------------
ResidualError:: ResidualError() : ErrorEstimator() { }
//-------------------------------------------------------------------------
QuadStdTriangleError:: QuadStdTriangleError() : ErrorEstimator() { }
//-------------------------------------------------------------------------
void ErrorEstimator:: adapt(Problem& problem, Real* energyError,
Real reqAbsGlobalPrec, Real* linENorm)
{
int i, markedElements = 0;
static Real accTimeError = 0.0;
Timer timer, accTimer;
elementError = new Vector<Real>(problem.Mesh()->noOfElements());
FORALL(*elementError,i) (*elementError)[i] = 0.0;
if (estimateError(problem) == True)
setRefinementFlags(*problem.Mesh(), reqAbsGlobalPrec, &markedElements);
if (markedElements == 0)
UniformRefinement(*problem.Mesh(), reqAbsGlobalPrec, &markedElements);
delete elementError;
elementError = 0;
*energyError = EnergyError;
if (linENorm != 0) *linENorm = LinENorm;
//if (Cmd.isTrue("L2Norm")) cout << "\tL2-Error: " << L2Error << "\n";
if (timeErrorEstimator)
{ cout << "\n\tError Estimation: "; timer.cpu(); cout.flush(); }
if (accTime)
{
accTimeError += accTimer.cpu(False);
cout << Form("\tAccumulated time: %1.2f sec.\n", accTimeError);
}
}
//-------------------------------------------------------------------------
void ErrorEstimator:: setRefinementFlags(MESH& mesh, Real reqGlobalPrec,
int* markedElems)
{
switch (refStrategy)
{
case maxValue: MaxValue(mesh, reqGlobalPrec, markedElems);
break;
case random: RandomRefinement(mesh, reqGlobalPrec, markedElems);
break;
case extrapolation: ExtrapolRefinement(mesh, reqGlobalPrec, markedElems);
break;
case uniform: UniformRefinement(mesh, reqGlobalPrec, markedElems);
break;
}
}
//-------------------------------------------------------------------------
void ErrorEstimator:: MaxValue(MESH& mesh, Real reqAbsGlobalPrec, int* marked)
{
int i, loops=0;
int markedElements;
PATCH* patch;
Real maxError = 0.0, limit;
FORALL (*elementError, i)
if ((*elementError)[i] > maxError) maxError = (*elementError)[i];
Real limitFactor = 0.25; Cmd.get("limitFactor", &limitFactor);
limit = limitFactor*maxError;
if (!(EnergyError == 0))
limit = Max(limit, reqAbsGlobalPrec/(EnergyError)*sqrt(limit*maxError));
do
{
markedElements = 0;
mesh.resetElemIter();
i = 0;
while ((patch = mesh.elemIterAll()))
{
++i;
if ( (*elementError)[i] > limit || DoRefine(i) )
{
patch->setMark();
++markedElements;
}
else patch->unMark();
}
limit *= 0.9;
if (loops < 100) ++loops;
else break;
}
while (markedElements < minRefinementRatio * mesh.noOfElements());
*marked = markedElements;
if (infoErrorEstimator)
{
cout << "\t -- Refinement: MaxValue\n";
cout << "\tMarked Elements: " << markedElements
<< " (" << int(100*markedElements/mesh.noOfElements()) << " %), "
<< " Limit reduction loops: " << loops <<"\n";
}
}
//-------------------------------------------------------------------------
void ErrorEstimator:: ExtrapolRefinement(MESH& mesh, Real reqAbsGlobalPrec,
int* marked)
{
int i, loops=0;
int markedElements;
PATCH* elem, *elemFather;
Real err, extrapolErr, fatherErr, limit = 0.0, maxError = 0.0;
mesh.resetElemIter();
i = 0;
while ((elem = mesh.elemIterAll()))
{
err = (*elementError)[++i];
elem->setError(err);
elemFather = elem->Father();
if (elemFather != 0)
{
fatherErr = elemFather->getError();
if (equal(fatherErr,0.0)) extrapolErr = 0.0;
else extrapolErr = err*err/fatherErr;
}
else extrapolErr = 0.0; // err
if (extrapolErr > err) extrapolErr = err; // rather unlikely
if (extrapolErr > limit) limit = extrapolErr;
if (err > maxError) maxError = err;
}
if (mesh.maxDepth==0) limit=0.5*maxError;
Real limitFactor = 0.5; Cmd.get("limitFactor", &limitFactor);
limit *= limitFactor; // 0.5 is empirical
if (!(EnergyError==0.0))
limit = Max(limit, reqAbsGlobalPrec/(EnergyError)*sqrt(limit*maxError));
do
{
markedElements = 0;
mesh.resetElemIter();
i = 0;
while ((elem = mesh.elemIterAll()))
{
++i;
if ((*elementError)[i] > limit || DoRefine(i))
{
elem->setMark();
++markedElements;
}
else elem->unMark();
}
limit *= 0.9;
if (loops < 100) ++loops;
else break;
}
while (markedElements < minRefinementRatio * mesh.noOfElements());
*marked = markedElements;
if (infoErrorEstimator)
{
cout << "\t -- Refinement: Extrapolation\n";
cout << "\tMarked Elements: " << markedElements
<< " (" << int(100*markedElements/mesh.noOfElements()) << " %), "
<< " Limit reduction loops: " << loops <<"\n";
}
}
//-------------------------------------------------------------------------
void ErrorEstimator:: UniformRefinement(MESH& mesh, Real /*reqAbsGlobalPrec*/,
int* marked)
{
PATCH* patch;
mesh.resetElemIter();
while ((patch = mesh.elemIterAll())) patch->setMark();
*marked = mesh.noOfElements();
if (infoErrorEstimator)
{
cout << "\t -- Refinement: uniform\n";
cout << "\tMarked Elements: " << mesh.noOfElements() << "\n";
}
}
//-------------------------------------------------------------------------
void ErrorEstimator:: RandomRefinement(MESH& mesh, Real /*reqAbsGlobalPrec*/,
int* marked)
{
int i;
PATCH* patch;
mesh.resetElemIter();
i = 0;
while ((patch = mesh.elemIterAll()))
{
++i;
if (i%10 == 1) patch->setMark();
else patch->unMark();
}
*marked = i;
if (infoErrorEstimator)
{
cout << "\t -- Refinement: random\n";
cout << "\tMarked Elements: "<< float(i)/10<<"\n";
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Bool ResidualError:: estimateError(Problem& /*problem*/)
{
cout << "\n\tError Est.: ResidualError not implemented\n" ;
cout.flush(); abort();
/*
int i, depth, noTri, noTri1, noTri2, ip;
TR* tr;
if (infoErrorEstimator) cout << "\n\tError Est.: Residual Error";
MESH& mesh = *problem.Mesh();
Element& element = *problem.element;
Interface& interface = *problem.interface;
Vector<Num>& u = problem.u;
noTri = 0;
mesh.setAllInnerNodes(&noTri); // Number triangles:
int spaceDim = element.spaceDim;
EDG* ed;
Vector<Num> gradElem(spaceDim); // local gradient
Vector<Num> uElem(element.noOfNodes); // local solution
Vector(int) node (element.noOfNodes);
Matrix<Num> grad (spaceDim, mesh.noOfTriangles);
Vector<Real> mat(mesh.noOfTriangles);
// calculate gradients:
Jacobian Jac;
ip = 1;
FORALL (mesh.trList, depth)
for (TrListProc tp1(mesh.trList[depth],all,&tr); tr; tr=tp1.All())
{
interface.getGlobalNodes(*tr, node);
FORALL(uElem,i) uElem[i] = u[node[i]];
Jac.compJinv(*tr);
element.gradient(*tr, Jac, uElem, gradElem, ip);
mat[tr->node] = problem.material->E(tr->classA);
for (i=1; i<=spaceDim; ++i)
grad(i,tr->node) = mat[tr->node] * gradElem[i];
// cout << "\n--- gradient: " << gradElem[1] << " " << gradElem[2];
}
Num flux1, flux2, js;
Real normal[4], length;
FORALL (mesh.edgList, depth)
for (EdgListProc ep(mesh.edgList[depth],all,&ed); ed; ed=ep.All())
{
if (ed->boundP == DIRICHLET) continue;
ed->compNormalVector(normal, &length);
noTri2 = 0;
if (ed->t1 == 0) noTri1 = ed->t2->node;
else
{
noTri1 = ed->t1->node;
if (ed->t2) noTri2 = ed->t2->node;
}
flux1 = flux2 = 0.0;
for (i=1; i<=spaceDim; ++i) flux1 += normal[i] * grad(i,noTri1);
if (noTri2)
for (i=1; i<=spaceDim; ++i) flux2 += normal[i] * grad(i,noTri2);
js = flux1-flux2;
js = js*js * length*length / 24.;
if (noTri1 && noTri2) js *= 0.5;
(*elementError)[noTri1] += abs(js / mat[noTri1]);
if (noTri2)
(*elementError)[noTri2] += abs(js / mat[noTri2]);
}
EnergyError = 0.0;
FORALL (*elementError, i) EnergyError += (*elementError)[i];
*/
return False;
}
//-------------------------------------------------------------------------
// returns RELATIVE error in per cent !
Bool QuadStdTriangleError:: estimateError(Problem& /*problem*/)
{
cout << "\n\tError Est.: QuadStdTriangleError not implemented\n" ;
cout.flush(); abort();
/*
int i, j, k, depth;
TR* t;
Num TotalLinEnergy=0.0, LinEnergy, TotalQuadEnergy=0.0, QuadEnergy;
MESH& mesh = *problem.Mesh();
Element& element = *problem.element;
Interface& interface = *problem.interface;
Vector<Num>& u = problem.u;
if (infoErrorEstimator) cout << "\n\tError Est.: QuadStdTriangleError";
int noTri = 0;
mesh.setAllInnerNodes(&noTri); // Number triangles
// assemble elements and calculate energies:
int dimElem = element.noOfNodes;
Matrix<Num> AElem(dimElem, dimElem);
Vector<Num> Au(dimElem), uElem(dimElem);
Vector(int) globalNodes(element.noOfNodes);
Jacobian Jac;
FORALL (mesh.trList, depth)
for (TrListProc trp2(mesh.trList[depth],all,&t); t; t=trp2.All())
{
Jac.compJinv(*t);
element.initAb(AElem);
interface.getGlobalNodes(*t, globalNodes);
problem.assembleErrorEstimator(*t, Jac, AElem);
if (problem.symmetry == sym) // set upper triangle
{
FORALL_ROWS(AElem,i)
for (j=1; j<i; ++j) AElem(j,i) = AElem(i,j);
}
FORALL(Au,i) Au[i] = 0.0;
FORALL(uElem,i) {
if (globalNodes[i]) uElem[i] = u[globalNodes[i]];
else uElem[i] = 0.0;
}
for (i=1; i<=3; ++i) // compute linear Au = A*u
{
if (globalNodes[i]==0) continue;
for (j=1; j<=3; ++j) Au[i] += AElem(i,j) * uElem[j];
}
LinEnergy = 0.0;
for (i=1; i<=3; ++i) LinEnergy += conj(uElem[i]) * Au[i];
LinEnergy *= 0.5;
for (i=1; i<=3; ++i) // compute quadratic Au = A*u
{
if (globalNodes[i]==0) continue;
for (j=4; j<=6; ++j) Au[i] += AElem(i,j) * uElem[j];
}
for (i=4; i<=6; ++i)
{
if (globalNodes[i]==0) continue;
for (j=1; j<=6; ++j) Au[i] += AElem(i,j) * uElem[j];
}
QuadEnergy = 0.0;
for (i=1; i<=6; ++i) QuadEnergy += Au[i] * uElem[i];
QuadEnergy *= 0.5;
TotalLinEnergy += LinEnergy;
TotalQuadEnergy += QuadEnergy;
(*elementError)[t->node] = 0.25 * abs(real(LinEnergy-QuadEnergy));
}
EnergyError = 0.25 * (TotalLinEnergy-TotalQuadEnergy);
EnergyError = abs(100 * EnergyError / TotalQuadEnergy);
LinENorm = abs(TotalLinEnergy);
if (infoErrorEstimator)
{
cout << "\n\nTotalQuadEnergy = " << TotalQuadEnergy
<< " TotalLinEnergy = " << TotalLinEnergy
<< " Error = " << EnergyError << " %\n";
}
// cout << "\n elementError: " << *elementError <<"\n";
*/
return False;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Bool DLY:: estimateError(Problem& problem)
{
if (infoErrorEstimator)
{
if (Cmd.isSet("ErrorEstim","ModDLY")) cout << "\n\tError Est.: ModDLY";
else cout << "\n\tError Est.: DLY";
}
Vector<Num> error(1);
Bool status = solveQuadDefectProblem(problem, error);
if (status == True) distributeEdgeErrors(problem, error);
return status;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Bool DLY:: solveQuadDefectProblem(Problem& problem, Vector<Num>& b)
{
int i, j;
Num sum;
PATCH* patch;
const Vector<Num>& u = problem.u;
const int edNodeM1 = u.high();
int edNode = edNodeM1;
interfaceDLY->setHighOrderNodes(&edNode);
if (edNode == edNodeM1) // only dirichlet-edges !
{
EnergyError = 0.0;
return False;
}
b.resize(edNodeM1+1,edNode);
Vector<Num> ADiag(edNodeM1+1,edNode);
FORALL(b,i) b[i] = ADiag[i] = 0.0;
// -- local assembling for DLY:
Vector<Num> uPrevOnNewMesh(1);
TransientProblem* trP;
if ((trP = problem.castToTransientProblem()))
{
if (trP->transientMode)
{
uPrevOnNewMesh.resize(b.low(),b.high());
trP->solutionToNewMesh(uPrevOnNewMesh, interfaceDLY);
}
}
const int low = problem.element->NoOfNodes()+1;
const int high = elementDLY->NoOfNodes();
Matrix<Num> AElem(high, high);
Vector<Num> bElem(high);
Vector<int> globalNodes(high);
Jacobian Jac;
problem.Mesh()->resetElemIter();
while ((patch = problem.Mesh()->elemIterAll()))
{
patch->compJinv(Jac);
elementDLY->initAb(AElem,bElem);
problem.assemble(*elementDLY, *patch, Jac, AElem, bElem,
elementDLY->DLYPattern, True);
interfaceDLY->getGlobalNodes(patch, globalNodes);
for (i=low; i<=high; ++i) // sum up diagonal and right hand side
{
if (globalNodes[i] == 0) continue;
ADiag[globalNodes[i]] += AElem(i,i);
b [globalNodes[i]] += bElem[i];
sum = 0.0;
for (j=1; j<low; ++j) // subtract 'coupling terms' from rhs
{
if (globalNodes[j]) sum += AElem(i,j) * u[globalNodes[j]];
}
b[globalNodes[i]] -= sum;
}
if (trP) // transient problem (implicit Euler assumed)
{
if (trP->transientMode)
{
for (i=low; i<=high; ++i) AElem(i,i) = 0.0;
elementDLY->assembleLumpedP(*patch, AElem, Jac,
elementDLY->DLYPattern);
for (i=low; i<=high; ++i)
if (globalNodes[i])
b[globalNodes[i]] += AElem(i,i) * uPrevOnNewMesh[globalNodes[i]];
}
}
}
FORALL(b,i) b[i] = conj(b[i])* b[i] / ADiag[i]; // error on edge
EnergyError = 0.0;
FORALL(b,i) EnergyError += Abs(b[i]); // sum up global Error
EnergyError *= 0.5;
return True;
}
//-------------------------------------------------------------------------
void DLY:: distributeEdgeErrors(Problem& problem, Vector<Num>& error)
{
int i, n, node, no=0;
PATCH* patch;
MESH* mesh = problem.Mesh();
mesh->resetElemIter();
const int nComp = problem.NComp();
if (problem.SpaceDim() == 1)
{
while ((patch = mesh->elemIterAll()))
{
++no;
if ((node=patch->getNode()))
for (n=1; n<=nComp; ++n) (*elementError)[no] += Abs(error[node++]);
}
}
else
{
while ((patch = mesh->elemIterAll()))
{
++no;
const Vector<EDG*>& edges = patch->getEdges();
FORALL(edges,i)
{
if (edges[i])
{
if ((node = edges[i]->getNode()))
{
for (n=1; n<=nComp; ++n)
(*elementError)[no] += Abs(error[node++]);
}
}
}
}
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Bool EFDLY:: estimateError(Problem& problem)
{
int i, node;
PATCH* patch;
if (infoErrorEstimator) cout << "\n\tError Est.: EF-DLY";
if (problem.spaceDim != 3)
{
cout << "\n*** Estimator EFDLY: dimension not 3\n\n"; cout.flush(); abort();
}
Vector<Num> error(1);
Bool status = solveQuadDefectProblem(problem, error);
if (status == True)
{
// -- distribute errors to elements
int noOfFaces = problem.element->NoOfFaces();
Vector<int> faceNodes(noOfFaces);
problem.Mesh()->resetElemIter();
int no = 0;
while ((patch = problem.Mesh()->elemIterAll()))
{
++no;
const Vector<EDG*>& edges = patch->getEdges();
FORALL(edges,i)
{
if (edges[i])
{
node = edges[i]->getNode();
if (node) (*elementError)[no] += Abs(error[node]);
// support: ~ 6 tetrahedra)
}
}
patch->getFaceNodes(faceNodes);
FORALL(faceNodes,i)
{
node = faceNodes[i];
if (node) (*elementError)[no] += 3.0 * Abs(error[node]);
// different support: 2 tetrahedra)
}
}
}
return status;
}
//-------------------------------------------------------------------------
syntax highlighted by Code2HTML, v. 0.9.1