Commit 6374f084 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

addBoundary*Operator, static updateAnimationFile and small changes in Initfile

parent 52419cef
......@@ -56,6 +56,7 @@
#include <mpi.h>
#endif
#include <boost/algorithm/string.hpp>
#include "boost/tuple/tuple.hpp"
#include "AMDiS_fwd.h"
......@@ -461,6 +462,7 @@ namespace AMDiS {
GRD_PSI,
GRD_PHI
};
}
#endif // AMDIS_GLOBAL_H
......
......@@ -8,13 +8,16 @@ using namespace std;
namespace AMDiS {
/// the small parser for the initfile. see description of read(Initfile&, istream&)
/// the small parser for the initfile. see description of
/// read(Initfile&,istream&)
struct Parser {
Parser(const string& line)
{
size_t pos = line.find(':');
if (pos == string::npos)
throw runtime_error("cannot find the delimiter ':' in line '" + line + "'");
if (pos == string::npos) {
throw runtime_error("cannot find the delimiter ':' in line "
"'" + line + "'");
}
name = line.substr(0, pos);
value = line.substr(pos + 1, line.length() - (pos + 1));
......@@ -27,6 +30,7 @@ namespace AMDiS {
string value;
};
Initfile* Initfile::singlett= NULL;
std::set<std::string> Initfile::fn_include_list;
......@@ -72,13 +76,16 @@ namespace AMDiS {
std::string sw(swap);
size_t pos0 = sw.find_first_not_of(whitespaces);
if (pos0 != std::string::npos && sw[pos0] != '%' && sw[pos0] != '#' && sw[pos0] != 0) {
if (pos0 != std::string::npos && sw[pos0] != '%'
&& sw[pos0] != '#' && sw[pos0] != 0) {
// parse line and extract map: tag->value
Parser parser(sw);
operator[](parser.name) = parser.value;
} else if (sw[pos0] == '#' && static_cast<size_t>(sw.find("#include")) == pos0) {
} else if (sw[pos0] == '#'
&& static_cast<size_t>(sw.find("#include")) == pos0) {
// include file by '#include "filename"' or '#include <filename>'
size_t pos = sw.find_first_not_of(whitespaces, std::string("#include").size() + 1);
size_t pos = sw.find_first_not_of(whitespaces,
std::string("#include").size() + 1);
size_t epos = 0;
std::string fn = "";
std::stringstream errorMsg;
......@@ -167,5 +174,4 @@ namespace AMDiS {
write(outFile);
}
}
} // end namespace AMDiS
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.
#ifndef INITFILE_H
#define INITFILE_H
......@@ -33,12 +44,14 @@ namespace AMDiS {
{}
};
struct NoDelim : std::runtime_error {
NoDelim(std::string m)
: std::runtime_error(m)
{}
};
struct WrongVectorFormat : std::runtime_error {
WrongVectorFormat(std::string m)
: std::runtime_error(m)
......@@ -124,7 +137,8 @@ namespace AMDiS {
}
/// return the delimiter or throw an exception if there is no known delimiter in value
/// return the delimiter or throw an exception if there is no known
/// delimiter in value
inline size_t checkDelim(const std::string& value, const std::string& delims)
{
size_t pos(std::string::npos);
......@@ -154,9 +168,11 @@ namespace AMDiS {
std::string val = trim(val_);
size_t pos = begBrackets.find(val[0]);
if (pos == std::string::npos)
throw WrongVectorFormat("cannot convert '" + val + "' into a list. No leading bracket found!");
throw WrongVectorFormat("cannot convert "
"'" + val + "' into a list. No leading bracket found!");
if (val[val.length() - 1] != endBrackets[pos])
throw WrongVectorFormat("begin and end bracket are different in value '" + val + "'");
throw WrongVectorFormat("begin and end bracket are different in"
" value '" + val + "'");
size_t oldPos = 1;
size_t curDelim = 0;
typedef typename Container::value_type ValueType;
......@@ -169,7 +185,7 @@ namespace AMDiS {
oldPos = pos + 1;
convert(curWord, swap);
c.push_back(swap);
pos= val.find(delims[curDelim], oldPos);
pos = val.find(delims[curDelim], oldPos);
}
//last entry
std::string curWord = val.substr(oldPos, val.length() - 1 - oldPos);
......@@ -204,15 +220,15 @@ namespace AMDiS {
using boost::numeric_cast;
mu::Parser parser;
parser.DefineConst(_T("M_PI"), m_pi);
parser.DefineConst(_T("M_E"), m_e);
try {
parser.SetExpr(valStr);
value = numeric_cast< T >(parser.Eval());
} catch(boost::bad_lexical_cast e) {
} catch (boost::bad_lexical_cast e) {
throw WrongValueFormat< T >(valStr);
} catch(boost::bad_numeric_cast e) {
} catch (boost::bad_numeric_cast e) {
throw WrongValueFormat< T >(valStr);
} catch (mu::Parser::exception_type &e) {
throw BadArithmeticExpression<T>(e.GetMsg(), valStr);
......@@ -224,16 +240,35 @@ namespace AMDiS {
inline void convert(const std::string valStr, T& value,
typename boost::enable_if< boost::is_enum< T > >::type* p = NULL)
{
unsigned int swap = 0;
int swap = 0;
try {
swap = boost::lexical_cast<unsigned int>(trim(valStr));
} catch(boost::bad_lexical_cast e) {
swap = boost::lexical_cast<int>(trim(valStr));
} catch (boost::bad_lexical_cast e) {
throw WrongValueFormat< T >(valStr);
}
value = static_cast< T >(swap);
}
/// convert special enums
inline void convert(const std::string valStr, Norm& value)
{
std::string swapStr = boost::to_upper_copy(valStr);
if (swapStr == "NO_NORM")
value = static_cast< Norm >(NO_NORM);
else if (swapStr == "H1_NORM")
value = static_cast< Norm >(H1_NORM);
else if (swapStr == "L2_NORM")
value = static_cast< Norm >(L2_NORM);
else {
int swap= 0;
convert(valStr, swap);
value = static_cast< Norm >(swap);
}
}
/// convert string to WorldVector
template< typename T >
inline void convert(const std::string valStr, WorldVector<T>& c)
......@@ -264,13 +299,15 @@ namespace AMDiS {
}
/// convert value of arbitrary type to string using stringstream and operator<< for type
/// convert value of arbitrary type to string using stringstream and
/// operator<< for type
template<typename T>
inline void convert(const T value, std::string& valStr)
{
std::stringstream ss;
ss.precision(6);
ss << value;
valStr= ss.str();
valStr = ss.str();
}
......@@ -286,10 +323,9 @@ namespace AMDiS {
} // end namespace InitfileInternal
/** The entry in an initfile. This helper class was constructed to allow calls like
* val = data.get(tag)
* for arbitrary types of val. At current stage, only double and bool is supported
/** The entry in an initfile. This helper class was constructed to allow calls
* like val = data.get(tag) for arbitrary types of val. At current stage, only
* double and bool is supported
*/
struct InitEntry {
///the value as string
......@@ -341,12 +377,14 @@ namespace AMDiS {
}
/** Basis data container as a map of tag on a value as strings. The container throws an exception, if the tag was not found.
/** Basis data container as a map of tag on a value as strings. The container
* throws an exception, if the tag was not found.
*/
struct Initfile : public std::map<std::string, std::string>
{
typedef std::map< std::string, std::string > super;
/// Exceptions
struct TagNotFound : std::invalid_argument {
TagNotFound(std::string m)
......@@ -354,28 +392,36 @@ namespace AMDiS {
{}
};
struct TagNotFoundBreak : std::invalid_argument { // print 'tag not found' and exit
struct TagNotFoundBreak : std::invalid_argument {
// print 'tag not found' and exit
TagNotFoundBreak(std::string m)
: std::invalid_argument(m)
{}
};
/** initialize init-file from file with filename in, read data and save it to singleton-map
/** initialize init-file from file with filename in, read data and save it
* to singleton-map
* @param in: filename string
*/
static void init(std::string in);
static void init(int print, string filename, const char *flags = NULL)
{
WARNING("Parameters::init(int,std::string,const char*) is depreciated. Use Parameters::init(std::string) instead!\n");
WARNING("Parameters::init(int,std::string,const char*) is depreciated. "
"Use Parameters::init(std::string) instead!\n");
init(filename);
}
/** Static get routine for getting parameter-values from init-file initialized in init()-method.
/** Static get routine for getting parameter-values from init-file
* initialized in init()-method.
* Cast the value to the desired type using std::stringstream.
* @param tag: The tag to look for
* @param value: The result.
* @param debugInfo: msgInfo for current parameter. (0..no printing, 1..print missing parameter info, 2..print parameter value) [optional]
* @param debugInfo: msgInfo for current parameter. (0..no printing,
* 1..print missing parameter info, 2..print parameter value) [optional]
*/
template<typename T>
static void get(const std::string tag, T& value, int debugInfo = -1)
......@@ -389,8 +435,10 @@ namespace AMDiS {
std::string valStr(singlett->checkedGet(tag));
valStr = trim(valStr);
convert(valStr, value);
if (debugInfo == 2)
std::cout << "Parameter '" << tag << "' initialized with: " << value << std::endl;
if (debugInfo == 2) {
std::cout << "Parameter '" << tag << "'"
<< " initialized with: " << value << std::endl;
}
} catch (TagNotFound ia) {
if (debugInfo >= 1)
std::cout << ia.what() << std::endl;
......@@ -431,9 +479,11 @@ namespace AMDiS {
// update msg parameters msgInfo, msgWait, paramInfo
singlett->getInternalParameters();
if (debugInfo == 2)
std::cout << "Parameter '" << tag << "' set to: " << value << std::endl;
std::cout << "Parameter '" << tag << "'"
<< " set to: " << value << std::endl;
}
/// add map tag->value to data in singleton
template< typename T >
static void add(const std::string tag, T& value, int debugInfo = -1)
......@@ -441,36 +491,43 @@ namespace AMDiS {
set(tag, value, debugInfo);
}
/// rescheduling parameter
static void readArgv(int argc, char **argv);
/// Returns specified info level
static int getMsgInfo()
{
return (singlett != NULL) ? singlett->msgInfo : 0;
}
/// Returns specified wait value
static int getMsgWait()
{
return (singlett != NULL) ? singlett->msgWait : 0;
}
/// Checks whether parameters are initialized. if not, call init()
static bool initialized()
{
return (singlett != NULL);
}
/// return pointer to singleton
static Initfile *getSingleton()
{
return singlett;
}
/// print all data in singleton to std::cout
static void printParameters();
/// clear data in singleton
static void clearData()
{
......@@ -478,6 +535,7 @@ namespace AMDiS {
singlett->clear();
}
/// save singlett-data to file with filename fn
static void save(std::string fn)
{
......@@ -485,7 +543,7 @@ namespace AMDiS {
singlett->write(fn);
}
protected:
protected:
Initfile()
: msgInfo(0),
msgWait(1),
......@@ -493,31 +551,37 @@ namespace AMDiS {
breakOnMissingTag(0)
{}
static void initIntern()
{
if (singlett == NULL)
singlett = new Initfile;
}
/// list of processed files
static std::set< std::string > fn_include_list;
/// pointer to the singleton that contains the data
static Initfile* singlett;
/// return the value of the given tag or throws an exception if the tag does not exist
/// return the value of the given tag or throws an exception if the tag
/// does not exist
std::string checkedGet(const std::string& tag) const
{
super::const_iterator it = find(tag);
if (it == end()) {
if (breakOnMissingTag == 0 || msgInfo <= 2)
throw TagNotFound("there is no tag '"+ tag + "'");
throw TagNotFound("there is no tag '" + tag + "'");
else
throw TagNotFoundBreak("required tag '"+ tag + "' not found");
throw TagNotFoundBreak("required tag '" + tag + "' not found");
}
return it->second;
}
/** Fill the initfile from an input stream.
* @param in: the stream to fill the data from.
* Current dataformat: tag:value
......
......@@ -1535,6 +1535,32 @@ namespace AMDiS {
}
void ProblemStatSeq::addBoundaryMatrixOperator(BoundaryType type,
Operator *op, int row, int col)
{
FUNCNAME("ProblemStat::addBoundaryOperator()");
RobinBC *robin =
new RobinBC(type, NULL, op, componentSpaces[row], componentSpaces[col]);
if (systemMatrix && (*systemMatrix)[row][col])
(*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(robin);
}
void ProblemStatSeq::addBoundaryVectorOperator(BoundaryType type,
Operator *op, int row)
{
FUNCNAME("ProblemStat::addBoundaryOperator()");
RobinBC *robin =
new RobinBC(type, op, NULL, componentSpaces[row]);
if (rhs)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
}
void ProblemStatSeq::assembleOnOneMesh(FiniteElemSpace *feSpace,
Flag assembleFlag,
DOFMatrix *matrix, DOFVector<double> *vector)
......
......@@ -289,6 +289,26 @@ namespace AMDiS {
/// Adds a periodic boundary condition.
virtual void addPeriodicBC(BoundaryType type, int row, int col);
/// add boundary operator to matrix side
virtual void addBoundaryMatrixOperator(BoundaryType type,
Operator *op, int row, int col);
virtual void addBoundaryMatrixOperator(BoundaryType type,
Operator &op, int row, int col)
{
addBoundaryMatrixOperator(type, &op, row, col);
}
/// add boundary operator to rhs (vector) side
virtual void addBoundaryVectorOperator(BoundaryType type,
Operator *op, int row);
virtual void addBoundaryVectorOperator(BoundaryType type,
Operator &op, int row)
{
addBoundaryVectorOperator(type, &op, row);
}
/** \brief
* This function assembles a DOFMatrix and a DOFVector for the case,
* the meshes from row and col FE-space are equal.
......
......@@ -102,7 +102,7 @@ namespace AMDiS {
}
/// Adds a new entry to a ParaView animation file.
int updateAnimationFile(std::string valueFilename,
static int updateAnimationFile(std::string valueFilename,
std::vector< std::string > *paraViewAnimationFrames,
std::string animationFilename);
protected:
......
......@@ -404,7 +404,7 @@ ElementLevelSet::calcIntersecNormal_3d(WorldVector<double> &normalVec)
int
ElementLevelSet::getVertexPos(const DimVec<double> barCoords)
{
double vertex_val;
double vertex_val= 0.0;
for (int i=0; i<=dim; ++i) {
if (barCoords[i] > 1-1.e-15) {
......
......@@ -6,6 +6,7 @@ bunnyMesh->global refinements: 0
bunny->mesh: bunnyMesh
bunny->dim: 2
bunny->components: 1
bunny->polynomial degree[0]: 1
bunny->solver: cg
......
dimension of world: 2
elliptMesh->macro file name: ./macro/macro.stand.2d
elliptMesh->macro file name: ./macro/macro.square.2d
elliptMesh->global refinements: 0
ellipt->mesh: elliptMesh
ellipt->dim: 2
ellipt->polynomial degree[0]: 1
ellipt->components: 1
ellipt->polynomial degree[1]: 1
ellipt->components: 2
ellipt->solver: cg
ellipt->solver: umfpack %cg
ellipt->solver->symmetric strategy: 1
ellipt->solver->store symbolic: 0
ellipt->solver->max iteration: 1000
ellipt->solver->tolerance: 1.e-8
ellipt->solver->info: 2
ellipt->solver->left precon: diag
ellipt->solver->right precon: no
ellipt->estimator[0]: residual
ellipt->estimator[0]->error norm: 1 % 1: H1_NORM, 2: L2_NORM
ellipt->estimator[0]->C0: 0.1 % constant of element residual
ellipt->estimator[0]->C1: 0.1 % constant of jump residual
%ellipt->estimator[0]: residual
%ellipt->estimator[0]->error norm: H1_NORM % 1: H1_NORM, 2: L2_NORM
%ellipt->estimator[0]->C0: 0.1 % constant of element residual
%ellipt->estimator[0]->C1: 0.1 % constant of jump residual
ellipt->marker[0]->strategy: 2 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
ellipt->marker[0]->MSGamma: 0.5
%ellipt->marker[0]->strategy: 0 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
%ellipt->marker[0]->MSGamma: 0.5
ellipt->adapt[0]->tolerance: 1e-4
ellipt->adapt[0]->refine bisections: 2
%ellipt->adapt[0]->tolerance: 1e-4
%ellipt->adapt[0]->refine bisections: 2
ellipt->adapt->max iteration: 10
ellipt->adapt->max iteration: 0
ellipt->output->filename: output/ellipt
ellipt->output->ParaView format: 1
......@@ -56,10 +56,10 @@ int main(int argc, char* argv[])
bunny.initialize(INIT_ALL);
// === create adapt info ===
AdaptInfo *adaptInfo = new AdaptInfo("bunny->adapt");
AdaptInfo adaptInfo("bunny->adapt");
// === create adapt ===
AdaptStationary *adapt = new AdaptStationary("bunny->adapt", &bunny, adaptInfo);
AdaptStationary adapt("bunny->adapt", bunny, adaptInfo);
// ===== create matrix operator =====
Operator matrixOperator(bunny.getFeSpace());
......@@ -75,7 +75,7 @@ int main(int argc, char* argv[])
bunny.addVectorOperator(&rhsOperator, 0);
// ===== start adaption loop =====
adapt->adapt();
adapt.adapt();
bunny.writeFiles(adaptInfo, true);
......
#include "AMDiS.h"
#include "ProblemStatDbg.h"
#include "Debug.h"
using namespace AMDiS;
using namespace std;
......@@ -7,6 +9,284 @@ using namespace std;
// ===== function definitions ================================================
// ===========================================================================
class MyProblemVec : public ProblemStatDbg
{
public:
MyProblemVec(std::string nameStr, ProblemIterationInterface *problemIteration = NULL)
: ProblemStatDbg(nameStr,problemIteration) {}
SolverMatrix<Matrix<DOFMatrix*> >* getSolverMatrix() { return &solverMatrix; }
void setSolverMatrix(SolverMatrix<Matrix<DOFMatrix*> >& solverMatrix_) { solverMatrix= solverMatrix_; }
inline DOFVector<double>* getRhsVector(int i = 0)
{
return rhs->getDOFVector(i);
}
};
class MyIteration : public StandardProblemIteration
{
public:
MyIteration(MyProblemVec *prob_) : StandardProblemIteration(prob_), prob(prob_) {}
~MyIteration() {}
Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION)
{ FUNCNAME("PIteration::oneIteration()");
Flag flag, markFlag;
if (toDo.isSet(MARK))
markFlag = prob->markElements(adaptInfo);
else
markFlag = 3;
// refine
if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
flag= prob->refineMesh(adaptInfo);
// coarsen
if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED))
flag|= prob->coarsenMesh(adaptInfo);
if (toDo.isSet(SOLVE)) {
prob->buildAfterCoarsen(adaptInfo, markFlag, true, true);
MSG("apply Dirichlet BC...\n");
bool inside= false;
WorldVector<double> p1;
p1[0]= 1.0;
p1[1]= -1.0;
for (int i= 0; i<prob->getSolution()->getSize(); ++i) {
DegreeOfFreedom idx1= getDOFidxAtCoords(prob->getSolution()->getDOFVector(i), p1, inside);
if(inside) {
MSG("set DBC to row %d...\n",idx1);
applyDbc(prob->getSystemMatrix(i,i)->getBaseMatrix(), *prob->getRhsVector(i), idx1, 0.0);
} else {
MSG("idx not found!!!\n");
}
}
// prob->writeDbgMatrix("ellipt_mat.dat");