Commit ade6bd99 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

correction of ProblemStat, i.e. an intermediate layer for standardProblemIteration added

parent 72990c0d
......@@ -360,43 +360,9 @@ namespace AMDiS {
using super::getName;
protected:
// /// Name of this problem.
// std::string name;
//
// /// Number of problem components
// int nComponents;
/// unqiue mesh-dimension for all problems
int dim;
// /** \brief
// * Number of problem meshes. If all components are defined on the same mesh,
// * this number is 1. Otherwise, this variable is the number of different meshes
// * within the problem.
// */
// int nMeshes;
// /// FE spaces of this problem.
// std::vector<FiniteElemSpace*> feSpaces;
//
// /// Meshes of this problem.
// std::vector<Mesh*> meshes;
// /** \brief
// * All actions of mesh refinement are performed by refinementManager.
// * If new refinement algorithms should be realized, one has to override
// * RefinementManager and give one instance of it to AdaptStationary.
// */
// RefinementManager *refinementManager;
//
// /** \brief
// * All actions of mesh coarsening are performed by coarseningManager.
// * If new coarsening algorithms should be realized, one has to override
// * CoarseningManager and give one instance of it to AdaptStationary.
// */
// CoarseningManager *coarseningManager;
std::vector<ProblemStatType*> problems;
};
......
......@@ -56,8 +56,7 @@ namespace AMDiS {
ProblemStatSeq::ProblemStatSeq(string nameStr,
ProblemIterationInterface *problemIteration)
: StandardProblemIteration(this),
name(nameStr),
: name(nameStr),
nComponents(-1),
nAddComponents(0),
nMeshes(0),
......@@ -826,18 +825,6 @@ namespace AMDiS {
}
Flag ProblemStatSeq::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
{
for (int i = 0; i < nComponents; i++)
if (adaptInfo->spaceToleranceReached(i))
adaptInfo->allowRefinement(false, i);
else
adaptInfo->allowRefinement(true, i);
return StandardProblemIteration::oneIteration(adaptInfo, toDo);
}
void ProblemStatSeq::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
bool asmMatrix, bool asmVector)
{
......
......@@ -55,8 +55,7 @@ namespace AMDiS {
* This class defines the stationary problem definition in sequential
* computations. For parallel computations, see \ref ParallelProblemStatBase.
**/
class ProblemStatSeq : public ProblemStatBase,
public StandardProblemIteration
class ProblemStatSeq : public ProblemStatBase
{
protected:
// Defines a mapping type from dof indices to world coordinates.
......@@ -695,7 +694,29 @@ namespace AMDiS {
};
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
typedef ProblemStatSeq ProblemStat;
struct ProblemStat : public ProblemStatSeq,
public StandardProblemIteration
{
typedef ProblemStatSeq super;
using super::getName;
ProblemStat(std::string nameStr,
ProblemIterationInterface *problemIteration = NULL)
: super(nameStr, problemIteration),
StandardProblemIteration(this)
{ }
Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo)
{
for (int i = 0; i < getNumComponents(); i++)
if (adaptInfo->spaceToleranceReached(i))
adaptInfo->allowRefinement(false, i);
else
adaptInfo->allowRefinement(true, i);
return StandardProblemIteration::oneIteration(adaptInfo, toDo);
}
};
#endif
}
......
......@@ -77,7 +77,29 @@ namespace AMDiS
} // end namespace Parallel
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
typedef Parallel::ParallelProblemStat ProblemStat;
struct ProblemStat : public Parallel::ParallelProblemStat,
public StandardProblemIteration
{
typedef Parallel::ParallelProblemStat super;
using super::getName;
ProblemStat(std::string nameStr,
ProblemIterationInterface *problemIteration = NULL)
: super(nameStr, problemIteration),
StandardProblemIteration(this)
{ }
Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo)
{
for (int i = 0; i < getNumComponents(); i++)
if (adaptInfo->spaceToleranceReached(i))
adaptInfo->allowRefinement(false, i);
else
adaptInfo->allowRefinement(true, i);
return StandardProblemIteration::oneIteration(adaptInfo, toDo);
}
};
#endif
} // end namespace AMDiS
......
#pragma once
#include <AMDiS.h>
#include <Helpers.h>
#include <functional>
#ifdef HAS_ALG_LIB
#include "interpolation.h" // from alglib
#endif
using namespace AMDiS;
// inline int maxima(const DOFVector<double> *vec, std::vector<WorldVector<double> >&x, std::vector<double> *values = NULL);
template <class Comparator>
int extrema(Comparator comp,
const DOFVector<double> *vec,
std::vector<WorldVector<double> >&x,
std::vector<double> *values = NULL)
{ FUNCNAME("extrema()");
const FiniteElemSpace* feSpace = vec->getFeSpace();
Mesh* mesh = feSpace->getMesh();
int num= 0; ///< number of extrema found
double tf= 0.7, tolF= 2.0;
// find extrema >= (<=) min_value + tf * (max_value - min_value)
Parameters::get("user parameter->threshold factor",tf);
// find extrema with distance >= tolF * h
Parameters::get("user parameter->maxima tolerance",tolF);
double minH, maxH;
int minLevel, maxLevel;
Helpers::calcMeshSizes(mesh, minH, maxH, minLevel, maxLevel) ;
double h= 2.0*minH;
double tol= tolF*h;
// selectDOF[i] = 1 ... if dof(i) is extrema
// selectDOF[i] = 0 ... otherwise
DOFVector<int> selectDOF(vec->getFeSpace(), "selection");
selectDOF.set(1);
double minNorm= min(*vec), maxNorm= max(*vec);
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Parallel::mpi::globalMin(minNorm);
Parallel::mpi::globalMax(maxNorm);
#endif
double threshold= (1.0-tf) * minNorm + tf * maxNorm;
const DOFAdmin* admin= feSpace->getAdmin();
const BasisFunction *basFcts= feSpace->getBasisFcts();
int nBasFcts= basFcts->getNumber();
std::vector<DegreeOfFreedom> localIndices(nBasFcts);
TraverseStack stack;
ElInfo *elInfo= stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
basFcts->getLocalIndices(elInfo->getElement(), admin, localIndices);
double maxValue= -1.e15; // TODO: use Comparator::init();
int maxIdx= -1;
// find vertex with maximal value that fulfills some conditions for each element
FixVec<WorldVector<double>, VERTEX> coords= elInfo->getMacroElement()->getCoord();
for (int i= 0; i < nBasFcts; ++i) {
double tempValue= (*vec)[localIndices[i]];
if (comp(tempValue, comp(maxValue, threshold) ? maxValue : threshold)) {
maxIdx= i;
maxValue= tempValue;
}
}
// mark other vertices than minimal vertex, to be definitiv no minima
for (int i= 0; i < nBasFcts; ++i) {
if (i != maxIdx)
selectDOF[localIndices[i]]= 0;
}
elInfo = stack.traverseNext(elInfo);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
// correct values at interface boundaries
Parallel::MeshDistributor::globalMeshDistributor->synchVector(functors::min_assign(), selectDOF);
#endif
x.clear();
if (values)
values->clear();
// collect coordinates to dofs with selectDOF[i] == 1
DOFVector<WorldVector<double> > coords(feSpace, "coords");
mesh->getDofIndexCoords(coords);
DOFIterator<WorldVector<double> > coordsIter(&coords, USED_DOFS);
DOFIterator<int> selectIter(&selectDOF, USED_DOFS);
DOFConstIterator<double> energyIter(vec, USED_DOFS);
for (coordsIter.reset(), selectIter.reset(), energyIter.reset();
!coordsIter.end();
++coordsIter, ++selectIter, ++energyIter)
{
if (*selectIter == 1) {
bool alreadyFound= false;
// filter minima that are too close to an other minima
for (std::vector<WorldVector<double> >::iterator it= x.begin();
it != x.end(); ++it)
{
if (norm(*it - (*coordsIter)) < tol) {
alreadyFound= true;
break;
}
}
if (!alreadyFound) {
x.push_back(*coordsIter);
if (values)
values->push_back(*energyIter);
}
}
}
num = x.size();
return num;
}
// ==========================================================================
#ifdef HAS_ALG_LIB
namespace alglib {
/// ALGLIB routines
///_______________________________________________________________________________________________
inline void dofvector2real_1d_array(DOFVector<double>& dofvector, alglib::real_1d_array& array)
{
array.setlength(dofvector.getUsedSize());
DOFIterator<double> vecIter(&dofvector,USED_DOFS);
int i = 0;
for(vecIter.reset(); !vecIter.end(); ++vecIter, ++i)
{
array(i) = (*vecIter);
}
}
inline void vector2real_1d_array(std::vector<double>& dofvector, alglib::real_1d_array& array)
{
array.setlength(dofvector.size());
for(size_t i = 0; i < dofvector.size(); ++i)
{
array(i) = dofvector[i];
}
}
inline void dofvector2real_2d_array(DOFVector<WorldVector<double> >& dofvector, alglib::real_2d_array& array)
{
array.setlength(dofvector.getUsedSize(), Global::getGeo(WORLD));
DOFIterator<WorldVector<double> > vecIter(&dofvector,USED_DOFS);
int i = 0;
for(vecIter.reset(); !vecIter.end(); ++vecIter, ++i)
{
for (int j=0; j < Global::getGeo(WORLD); ++j)
array(i,j) = (*vecIter)[j];
}
}
inline void vector2real_2d_array(std::vector<WorldVector<double> >& dofvector, alglib::real_2d_array& array)
{
array.setlength(dofvector.size(), Global::getGeo(WORLD));
for(size_t i = 0; i < dofvector.size(); ++i)
{
for (int j = 0; j < Global::getGeo(WORLD); ++j)
array(i,j) = dofvector[i][j];
}
}
inline void vector2real_2d_array(std::vector<double>& dofvector, alglib::real_2d_array& array)
{
array.setlength(dofvector.size(), 1);
for(size_t i = 0; i < dofvector.size(); ++i)
{
array(i,0) = dofvector[i];
}
}
// =====================================================================================================
inline void function_gauss(const alglib::real_1d_array &c, const alglib::real_1d_array &x, double &func, void *ptr)
{
double nrm = sqr(x[0]-c[2])+sqr(x[1]-c[3]);
func = c[0] + c[1]*exp(-nrm/c[4]);
}
inline void gradient_gauss(const alglib::real_1d_array &c, const alglib::real_1d_array &x, double &func, alglib::real_1d_array &grad, void *ptr)
{
double nrm = sqr(x[0]-c[2])+sqr(x[1]-c[3]);
func = c[0] + c[1]*exp(-nrm/c[4]);
grad[0] = 1.0;
grad[1] = exp(-nrm/c[4]);
grad[2] = -(c[1]*exp(-nrm/c[4])*(2.0*(c[2] - x[0])))/c[4];
grad[3] = -(c[1]*exp(-nrm/c[4])*(2.0*(c[3] - x[1])))/c[4];
grad[4] = (c[1]*exp(-nrm/c[4])*nrm)/sqr(c[4]);
}
inline void nonlinearGaussFit(DOFVector<double>& pfc,
WorldVector<double> &p,
double &c0, double &c1, double &c4 )
{
using namespace alglib;
extensions::KD_Tree_Dof* tree = extensions::provideKDTree(pfc.getFeSpace());
double radius = m_pi/sqrt(3.0);
std::vector<std::pair<DegreeOfFreedom, double> > indices;
tree->index->radiusSearch(p.begin(), sqr(radius), indices, nanoflann::SearchParams(10));
std::vector<WorldVector<double> > points;
std::vector<double> values;
for (size_t i = 0; i < indices.size(); i++) {
points.push_back(tree->m_data[indices[i].first]);
values.push_back(pfc[indices[i].first]);
}
real_2d_array x;
vector2real_2d_array(points, x);
real_1d_array y;
vector2real_1d_array(values, y);
int numParams = 5;
real_1d_array c; c.setlength(numParams);
c(0) = c0;
c(1) = c1;
c(2) = p[0];
c(3) = p[1];
c(4) = c4;
double tol = abs(radius*0.2);
real_1d_array lb; lb.setlength(numParams);
lb(0) = 0.0;
lb(1) = 1.e-4;
lb(2) = c(2) - tol;
lb(3) = c(3) - tol;
lb(4) = 1.e-4;
real_1d_array ub; ub.setlength(numParams);
ub(0) = numeric_limits<double>::infinity();
ub(1) = numeric_limits<double>::infinity();
ub(2) = c(2) + tol;
ub(3) = c(3) + tol;
ub(4) = numeric_limits<double>::infinity();
double epsf = 0.0;
double epsx = 0.0;
ae_int_t maxits = 0;
ae_int_t info;
lsfitstate state;
lsfitreport rep;
// double diffstep = 1.e-6;
lsfitcreatefg(x, y, c, points.size(), Global::getGeo(WORLD), numParams, true, state);
lsfitsetcond(state, epsf, epsx, maxits);
lsfitsetbc(state, lb, ub);
alglib::lsfitfit(state, alglib::function_gauss, alglib::gradient_gauss);
lsfitresults(state, info, c, rep);
WorldVector<double> p0 = p;
p[0] = c(2); p[1] = c(3);
c0 = c(0); c1 = c(1); c4 = c(4);
// std::cout<< "x: [" << p0[0] << ", " << p0[1] << "] => [" << p[0] << ", " << p[1] << "]\n";
}
// =====================================================================================================
inline void function_gauss1D(const alglib::real_1d_array &c, const alglib::real_1d_array &x, double &func, void *ptr)
{
double nrm = sqr(x[0]-c[2]);
func = c[0] + c[1]*exp(-nrm/c[3]);
}
inline void gradient_gauss1D(const alglib::real_1d_array &c, const alglib::real_1d_array &x, double &func, alglib::real_1d_array &grad, void *ptr)
{
double nrm = sqr(x[0]-c[2]);
func = c[0] + c[1]*exp(-nrm/c[3]);
grad[0] = 1.0;
grad[1] = exp(-nrm/c[3]);
grad[2] = -(c[1]*exp(-nrm/c[3])*(2.0*(c[2] - x[0])))/c[3];
grad[3] = (c[1]*exp(-nrm/c[3])*nrm)/sqr(c[3]);
}
inline void nonlinearGaussFit1D(DOFVector<double>& pfc,
WorldVector<double> &p,
double &c0, double &c1, double &c3 )
{
using namespace alglib;
extensions::KD_Tree_Dof* tree = extensions::provideKDTree(pfc.getFeSpace());
double radius = m_pi/sqrt(3.0);
std::vector<std::pair<DegreeOfFreedom, double> > indices;
tree->index->radiusSearch(p.begin(), sqr(radius), indices, nanoflann::SearchParams(10));
std::vector<WorldVector<double> > points;
std::vector<double> values;
for (size_t i = 0; i < indices.size(); i++) {
points.push_back(tree->m_data[indices[i].first]);
values.push_back(pfc[indices[i].first]);
}
real_2d_array x;
vector2real_2d_array(points, x);
real_1d_array y;
vector2real_1d_array(values, y);
int numParams = 4;
real_1d_array c; c.setlength(numParams);
c(0) = c0;
c(1) = c1;
c(2) = p[0];
c(3) = c3;
double tol = abs(radius*0.2);
real_1d_array lb; lb.setlength(numParams);
lb(0) = 0.0;
lb(1) = 1.e-4;
lb(2) = c(2) - tol;
lb(3) = 1.e-4;
real_1d_array ub; ub.setlength(numParams);
ub(0) = numeric_limits<double>::infinity();
ub(1) = numeric_limits<double>::infinity();
ub(2) = c(2) + tol;
ub(3) = numeric_limits<double>::infinity();
double epsf = 0.0;
double epsx = 0.0;
ae_int_t maxits = 0;
ae_int_t info;
lsfitstate state;
lsfitreport rep;
// double diffstep = 1.e-6;
lsfitcreatefg(x, y, c, points.size(), Global::getGeo(WORLD), numParams, true, state);
lsfitsetcond(state, epsf, epsx, maxits);
lsfitsetbc(state, lb, ub);
alglib::lsfitfit(state, alglib::function_gauss, alglib::gradient_gauss);
lsfitresults(state, info, c, rep);
WorldVector<double> p0 = p;
p[0] = c(2);
c0 = c(0); c1 = c(1); c3 = c(3);
}
} // end namespace alglib
#endif
inline std::string getFilename(std::string name, double time, std::string ext = ".dat")
{
std::string filename = "";
int indexLength, indexDecimals;
Parameters::get(name + "->filename", filename);
Parameters::get(name + "->index length", indexLength);
Parameters::get(name + "->index decimals", indexDecimals);
TEST_EXIT(indexLength <= 99)("index lenght > 99\n");
TEST_EXIT(indexDecimals <= 97)("index decimals > 97\n");
TEST_EXIT(indexDecimals < indexLength)("index length <= index decimals\n");
char formatStr[9];
char timeStr[20];
sprintf(formatStr, "%%0%d.%df", indexLength, indexDecimals);
sprintf(timeStr, formatStr, time);
filename += timeStr + ext;
return filename;
}
inline void unique_value(std::vector<WorldVector<double> > &vec,
std::vector<double> &val,
double tol,
std::vector<unsigned> *ind = NULL)
{
compareTol<WorldVector<double> > comp(tol);
size_t newVec = 0;
for (size_t i = 0; i < vec.size(); ++i) {
bool inNew = false;
for (size_t j = 0; j < newVec; ++j) {
if (comp(vec[i], vec[j])) {
if (val[i] > val[j]+FLT_TOL) {
vec[j] = vec[i];
val[j] = val[i];
}
inNew = true;
break;
}
}
if (!inNew) {
vec[newVec] = vec[i];
val[newVec] = val[i];
newVec++;
if (ind)
ind->push_back(i);
}
}
vec.erase(vec.begin()+newVec,vec.end());
}
inline void writeMaxima(std::string name, AdaptInfo* adaptInfo, DOFVector<double>* rho, std::vector<WorldVector<double> >& x)
{
bool useFitting = false;
bool calcMaxima = false;
Parameters::get(name + "->calculate maxima", calcMaxima);
Parameters::get(name + "->use maxima fitting", useFitting);
x.clear();
if (!calcMaxima)
return;
std::vector<double> val;
extrema(std::greater<double>(), rho, x, &val);
if (useFitting) {
double c0=0.0, c1=1.0, c4=1.0; // parameters for nonlinear fitting
for (size_t i = 0; i < x.size(); i++) {
#ifdef HAS_ALG_LIB
WorldVector<double> new_x = x[i];
try {
if (Global::getGeo(WORLD) == 2)
alglib::nonlinearGaussFit(*vec, new_x, c0, c1, c4);
else if (Global::getGeo(WORLD) == 1)
alglib::nonlinearGaussFit1D(*vec, new_x, c0, c1, c4);
else {
WARNING("No fitting implemented for 3d!\n");
break;
}
x[i] = new_x;
} catch(...) {}
#endif
}
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Parallel::StdMpi<std::vector<WorldVector<double> > > stdMpi(MPI::COMM_WORLD, true);