From 6374f084a08571f5bbaf8a624a16a1b9d8c68c6d Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Wed, 15 Jun 2011 07:42:04 +0000
Subject: [PATCH] addBoundary*Operator, static updateAnimationFile and small
changes in Initfile
---
AMDiS/src/Global.h | 2 +
AMDiS/src/Initfile.cc | 96 ++++----
AMDiS/src/Initfile.h | 364 ++++++++++++++++------------
AMDiS/src/ProblemStat.cc | 26 ++
AMDiS/src/ProblemStat.h | 20 ++
AMDiS/src/io/VtkWriter.h | 2 +-
AMDiS/src/reinit/ElementLevelSet.cc | 2 +-
demo/init/bunny.init | 1 +
demo/init/ellipt.dat.2d | 27 ++-
demo/src/bunny.cc | 6 +-
demo/src/ellipt.cc | 302 ++++++++++++++++++++++-
11 files changed, 626 insertions(+), 222 deletions(-)
diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h
index 0be17a4b..8bdb675c 100644
--- a/AMDiS/src/Global.h
+++ b/AMDiS/src/Global.h
@@ -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
diff --git a/AMDiS/src/Initfile.cc b/AMDiS/src/Initfile.cc
index 595795aa..d0aef6f9 100644
--- a/AMDiS/src/Initfile.cc
+++ b/AMDiS/src/Initfile.cc
@@ -8,25 +8,29 @@ 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));
// remove everything after the %
pos = value.find('%');
if (pos != string::npos)
- value = value.substr(0, pos);
+ value = value.substr(0, pos);
}
string name;
string value;
};
+
Initfile* Initfile::singlett= NULL;
std::set<std::string> Initfile::fn_include_list;
@@ -39,7 +43,7 @@ namespace AMDiS {
fn_include_list.clear();
singlett->read(in);
singlett->getInternalParameters();
-
+
// initialize global strcutures using parameters
Global::init();
}
@@ -53,8 +57,8 @@ namespace AMDiS {
std::ifstream inputFile;
inputFile.open(fn.c_str(), std::ios::in);
if (!inputFile.is_open())
- throw runtime_error("init-file cannot be opened for reading");
-
+ throw runtime_error("init-file cannot be opened for reading");
+
fn_include_list.insert(fn);
read(inputFile);
}
@@ -71,34 +75,37 @@ namespace AMDiS {
std::string whitespaces = " \t\r\f";
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) {
- // 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) {
- // include file by '#include "filename"' or '#include <filename>'
- size_t pos = sw.find_first_not_of(whitespaces, std::string("#include").size() + 1);
- size_t epos = 0;
- std::string fn = "";
- std::stringstream errorMsg;
- switch (char c= swap[pos++]) {
- case '<':
- c= '>';
- case '\"':
- whitespaces += c;
- epos = sw.find_first_of(whitespaces, pos);
- fn = sw.substr(pos, epos - pos);
-
- if (sw[epos]!=c) {
- errorMsg << "filename in #include not terminated by " << c;
- throw std::runtime_error(errorMsg.str());
- }
- break;
- default:
- throw std::runtime_error("no filename given for #include");
- }
- read(fn);
+
+ 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) {
+ // include file by '#include "filename"' or '#include <filename>'
+ size_t pos = sw.find_first_not_of(whitespaces,
+ std::string("#include").size() + 1);
+ size_t epos = 0;
+ std::string fn = "";
+ std::stringstream errorMsg;
+ switch (char c= swap[pos++]) {
+ case '<':
+ c= '>';
+ case '\"':
+ whitespaces += c;
+ epos = sw.find_first_of(whitespaces, pos);
+ fn = sw.substr(pos, epos - pos);
+
+ if (sw[epos]!=c) {
+ errorMsg << "filename in #include not terminated by " << c;
+ throw std::runtime_error(errorMsg.str());
+ }
+ break;
+ default:
+ throw std::runtime_error("no filename given for #include");
+ }
+ read(fn);
}
in.getline(swap, line_length);
}
@@ -109,8 +116,8 @@ namespace AMDiS {
{
for (int i = 0; i < argc; i++) {
if (strcmp("-rs", argv[i]) == 0) {
- std::string input(argv[i + 1]);
- add("argv->rs", input, 0);
+ std::string input(argv[i + 1]);
+ add("argv->rs", input, 0);
}
}
}
@@ -122,15 +129,15 @@ namespace AMDiS {
int val = 0;
get("level of information", val, 0);
msgInfo = val;
-
+
val = 1;
get("WAIT", val, 0);
msgWait = val;
-
+
val = 1;
get("parameter information", val, 0);
paramInfo = val;
-
+
val = 0;
get("break on missing tag", val, 0);
breakOnMissingTag = val;
@@ -145,7 +152,7 @@ namespace AMDiS {
{
initIntern();
for (Initfile::iterator it = singlett->begin(); it != singlett->end(); it++)
- std::cout << (*it).first << " => " << (*it).second << std::endl;
+ std::cout << (*it).first << " => " << (*it).second << std::endl;
}
@@ -153,7 +160,7 @@ namespace AMDiS {
void Initfile::write(ostream& out)
{
for (Initfile::iterator it = begin() ; it!=end(); it++)
- out << (*it).first << ": " << (*it).second << std::endl;
+ out << (*it).first << ": " << (*it).second << std::endl;
}
@@ -164,8 +171,7 @@ namespace AMDiS {
outFile.open(fn.c_str(), std::ios::out);
if (!outFile.is_open())
throw runtime_error("init-file cannot be opened for writing");
-
+
write(outFile);
}
-
-}
+} // end namespace AMDiS
diff --git a/AMDiS/src/Initfile.h b/AMDiS/src/Initfile.h
index 90538216..e23c0210 100644
--- a/AMDiS/src/Initfile.h
+++ b/AMDiS/src/Initfile.h
@@ -1,3 +1,14 @@
+//
+// 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
@@ -25,23 +36,25 @@
namespace AMDiS {
namespace InitfileInternal {
-
+
/// Exceptions
struct WrongVectorSize : std::runtime_error {
WrongVectorSize(std::string m)
- : std::runtime_error(m)
+ : std::runtime_error(m)
{}
};
+
struct NoDelim : std::runtime_error {
NoDelim(std::string m)
- : std::runtime_error(m)
+ : std::runtime_error(m)
{}
};
+
struct WrongVectorFormat : std::runtime_error {
WrongVectorFormat(std::string m)
- : std::runtime_error(m)
+ : std::runtime_error(m)
{}
};
@@ -49,33 +62,33 @@ namespace AMDiS {
struct WrongValueFormat : std::runtime_error {
static std::string name(int)
{
- return "int";
+ return "int";
}
-
+
static std::string name(bool)
{
- return "bool";
+ return "bool";
}
-
+
static std::string name(double)
{
- return "double";
+ return "double";
}
-
+
static std::string name(unsigned int)
{
- return "unsigned int";
+ return "unsigned int";
}
template<typename G>
static std::string name(G)
{
- return std::string(typeid(G).name());
+ return std::string(typeid(G).name());
}
-
+
WrongValueFormat(std::string value)
- : std::runtime_error(std::string("cannot convert '") +
- value + std::string("' into <") + name(T()) + ">")
+ : std::runtime_error(std::string("cannot convert '") +
+ value + std::string("' into <") + name(T()) + ">")
{}
};
@@ -83,38 +96,38 @@ namespace AMDiS {
struct BadArithmeticExpression : std::runtime_error {
static std::string name(int)
{
- return "int";
+ return "int";
}
-
+
static std::string name(bool)
{
- return "bool";
+ return "bool";
}
-
+
static std::string name(double)
{
- return "double";
+ return "double";
}
-
+
static std::string name(unsigned int)
{
- return "unsigned int";
+ return "unsigned int";
}
-
+
template<typename G>
static std::string name(G)
{
- return std::string(typeid(G).name());
+ return std::string(typeid(G).name());
}
-
+
BadArithmeticExpression(std::string m, std::string value)
- : std::runtime_error(std::string("cannot evaluate expression '") +
- value + std::string("' into <") + name(T()) +
- ">\nParser message: '" + m + "'")
+ : std::runtime_error(std::string("cannot evaluate expression '") +
+ value + std::string("' into <") + name(T()) +
+ ">\nParser message: '" + m + "'")
{}
};
-
+
/// trim std::string
inline std::string trim(const std::string& oldStr)
{
@@ -124,14 +137,15 @@ 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);
for (unsigned i = 0; i < delims.length(); i++) {
- pos = value.find(delims[i]);
- if (pos != std::string::npos)
- return i;
+ pos = value.find(delims[i]);
+ if (pos != std::string::npos)
+ return i;
}
throw NoDelim("cannot detect the delimiter in " + value);
return 0;
@@ -139,12 +153,12 @@ namespace AMDiS {
/** parse an container from tag tag. The Container must have the properties:
- * - type value_type
- * - member function push_back
- */
+ * - type value_type
+ * - member function push_back
+ */
template< typename Container >
inline void getContainer(const std::string val_, Container& c)
- {
+ {
// accepted brackets and delimiters for vector input
std::string begBrackets= "{[(";
std::string endBrackets= "}])";
@@ -154,35 +168,37 @@ 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;
ValueType swap;
try {
- curDelim = checkDelim(val, delims);
- pos = val.find(delims[curDelim], oldPos);
- while (pos != std::string::npos) {
- std::string curWord = val.substr(oldPos, pos - oldPos);
- oldPos = pos + 1;
- convert(curWord, swap);
- c.push_back(swap);
- pos= val.find(delims[curDelim], oldPos);
- }
- //last entry
- std::string curWord = val.substr(oldPos, val.length() - 1 - oldPos);
- convert(curWord, swap);
- c.push_back(swap);
+ curDelim = checkDelim(val, delims);
+ pos = val.find(delims[curDelim], oldPos);
+ while (pos != std::string::npos) {
+ std::string curWord = val.substr(oldPos, pos - oldPos);
+ oldPos = pos + 1;
+ convert(curWord, swap);
+ c.push_back(swap);
+ pos = val.find(delims[curDelim], oldPos);
+ }
+ //last entry
+ std::string curWord = val.substr(oldPos, val.length() - 1 - oldPos);
+ convert(curWord, swap);
+ c.push_back(swap);
} catch (NoDelim nd) {
- std::string curWord = val.substr(1, val.length() - 2);
- curWord = trim(curWord);
- if (curWord.length() > 0) {
- // container with one entry
- convert(curWord, swap);
- c.push_back(swap);
- }
+ std::string curWord = val.substr(1, val.length() - 2);
+ curWord = trim(curWord);
+ if (curWord.length() > 0) {
+ // container with one entry
+ convert(curWord, swap);
+ c.push_back(swap);
+ }
}
}
@@ -197,43 +213,62 @@ namespace AMDiS {
/// convert string to intrinsic type
template<typename T>
inline void convert(const std::string valStr, T& value ,
- typename boost::enable_if<boost::is_pod<T> >::type* p = NULL ,
- typename boost::disable_if<boost::is_enum<T> >::type* p2 = NULL)
+ typename boost::enable_if<boost::is_pod<T> >::type* p = NULL ,
+ typename boost::disable_if<boost::is_enum<T> >::type* p2 = NULL)
{
using boost::lexical_cast;
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) {
- throw WrongValueFormat< T >(valStr);
- } catch(boost::bad_numeric_cast e) {
- throw WrongValueFormat< T >(valStr);
+ parser.SetExpr(valStr);
+ value = numeric_cast< T >(parser.Eval());
+ } catch (boost::bad_lexical_cast e) {
+ throw WrongValueFormat< T >(valStr);
+ } catch (boost::bad_numeric_cast e) {
+ throw WrongValueFormat< T >(valStr);
} catch (mu::Parser::exception_type &e) {
- throw BadArithmeticExpression<T>(e.GetMsg(), valStr);
+ throw BadArithmeticExpression<T>(e.GetMsg(), valStr);
}
}
template<typename T>
inline void convert(const std::string valStr, T& value,
- typename boost::enable_if< boost::is_enum< T > >::type* p = NULL)
+ 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) {
- throw WrongValueFormat< T >(valStr);
+ 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)
@@ -241,10 +276,10 @@ namespace AMDiS {
std::vector<T> temp_vec;
getContainer(valStr, temp_vec);
if (static_cast<int>(temp_vec.size()) != c.getSize())
- throw WrongVectorSize("wrong number of entries for WorldVector");
-
+ throw WrongVectorSize("wrong number of entries for WorldVector");
+
for (unsigned i = 0; i < temp_vec.size(); i++)
- c[i] = temp_vec[i];
+ c[i] = temp_vec[i];
}
@@ -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();
}
@@ -280,24 +317,23 @@ namespace AMDiS {
{
std::vector<T> temp_vec(c.getSize());
for (unsigned i = 0; i < temp_vec.size(); i++)
- temp_vec[i] = c[i];
+ temp_vec[i] = c[i];
convert(temp_vec, valStr);
}
} // 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
std::string valStr;
/// initialize with value as string
InitEntry(std::string v = "")
- : valStr(v)
+ : valStr(v)
{}
/// cast string to type T
@@ -341,59 +377,71 @@ 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)
- : std::invalid_argument(m)
+ : std::invalid_argument(m)
{}
};
- 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)
+ : std::invalid_argument(m)
{}
};
-
- /** initialize init-file from file with filename in, read data and save it to singleton-map
- * @param in: filename string
- */
+
+
+ /** 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");
- init(filename);
+ 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.
- * 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]
- */
+
+
+ /** 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]
+ */
template<typename T>
static void get(const std::string tag, T& value, int debugInfo = -1)
{
using namespace InitfileInternal;
initIntern();
if (debugInfo == -1)
- debugInfo = singlett->getMsgInfo();
-
+ debugInfo = singlett->getMsgInfo();
+
try {
- std::string valStr(singlett->checkedGet(tag));
- valStr = trim(valStr);
- convert(valStr, value);
- if (debugInfo == 2)
- std::cout << "Parameter '" << tag << "' initialized with: " << value << std::endl;
+ std::string valStr(singlett->checkedGet(tag));
+ valStr = trim(valStr);
+ convert(valStr, value);
+ if (debugInfo == 2) {
+ std::cout << "Parameter '" << tag << "'"
+ << " initialized with: " << value << std::endl;
+ }
} catch (TagNotFound ia) {
- if (debugInfo >= 1)
- std::cout << ia.what() << std::endl;
+ if (debugInfo >= 1)
+ std::cout << ia.what() << std::endl;
}
}
@@ -405,16 +453,16 @@ namespace AMDiS {
int debugInfo = singlett->getMsgInfo();
InitEntry result;
try {
- std::string valStr(singlett->checkedGet(tag));
- valStr = trim(valStr);
- result = InitEntry(valStr);
+ std::string valStr(singlett->checkedGet(tag));
+ valStr = trim(valStr);
+ result = InitEntry(valStr);
} catch (TagNotFound ia) {
- if (debugInfo >= 1)
- std::cout << ia.what() << std::endl;
+ if (debugInfo >= 1)
+ std::cout << ia.what() << std::endl;
}
return result;
}
-
+
/// update map tag->value_old to tag->value in singleton
template<typename T>
@@ -423,17 +471,19 @@ namespace AMDiS {
using namespace InitfileInternal;
initIntern();
if (debugInfo == -1)
- debugInfo = singlett->getMsgInfo();
-
+ debugInfo = singlett->getMsgInfo();
+
std::string swap = "";
convert(value, swap);
(*singlett)[trim(tag)] = swap;
// 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,43 +491,51 @@ 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()
{
initIntern();
singlett->clear();
}
-
+
+
/// save singlett-data to file with filename fn
static void save(std::string fn)
{
@@ -485,64 +543,70 @@ namespace AMDiS {
singlett->write(fn);
}
- protected:
+protected:
Initfile()
- : msgInfo(0),
- msgWait(1),
- paramInfo(1),
- breakOnMissingTag(0)
+ : msgInfo(0),
+ msgWait(1),
+ paramInfo(1),
+ breakOnMissingTag(0)
{}
-
+
+
static void initIntern()
{
if (singlett == NULL)
- singlett = new Initfile;
+ 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 + "'");
- else
- throw TagNotFoundBreak("required tag '"+ tag + "' not found");
+ if (breakOnMissingTag == 0 || msgInfo <= 2)
+ throw TagNotFound("there is no tag '" + tag + "'");
+ else
+ 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
- * Comment char: percent '%'
- * Include files: #include "filename" or #include <filename>
- */
+ * @param in: the stream to fill the data from.
+ * Current dataformat: tag:value
+ * Comment char: percent '%'
+ * Include files: #include "filename" or #include <filename>
+ */
void read(std::istream& in);
/// Fill the initfile from a file with filename fn
void read(std::string fn);
-
+
/** Write data-map to initfile
- * @param out: the stream to fill the data in.
- */
+ * @param out: the stream to fill the data in.
+ */
void write(std::ostream& out);
/// Write data-map to initfile with filename fn
void write(std::string fn);
-
+
/// read parameters for msgInfo, msgWait, paramInfo
void getInternalParameters();
-
+
int msgInfo, msgWait, paramInfo, breakOnMissingTag;
};
-
+
typedef Initfile Parameters;
} // end namespace AMDiS
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index 047c9f77..97e56a52 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -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)
diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h
index 6e870973..a04a3aa5 100644
--- a/AMDiS/src/ProblemStat.h
+++ b/AMDiS/src/ProblemStat.h
@@ -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.
diff --git a/AMDiS/src/io/VtkWriter.h b/AMDiS/src/io/VtkWriter.h
index a10939b0..430a3efc 100644
--- a/AMDiS/src/io/VtkWriter.h
+++ b/AMDiS/src/io/VtkWriter.h
@@ -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:
diff --git a/AMDiS/src/reinit/ElementLevelSet.cc b/AMDiS/src/reinit/ElementLevelSet.cc
index 9c3e6cc7..64ae1956 100644
--- a/AMDiS/src/reinit/ElementLevelSet.cc
+++ b/AMDiS/src/reinit/ElementLevelSet.cc
@@ -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) {
diff --git a/demo/init/bunny.init b/demo/init/bunny.init
index 0e496538..7c67faf2 100644
--- a/demo/init/bunny.init
+++ b/demo/init/bunny.init
@@ -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
diff --git a/demo/init/ellipt.dat.2d b/demo/init/ellipt.dat.2d
index e01f7688..3a5ce968 100644
--- a/demo/init/ellipt.dat.2d
+++ b/demo/init/ellipt.dat.2d
@@ -1,32 +1,35 @@
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
diff --git a/demo/src/bunny.cc b/demo/src/bunny.cc
index 883e7d18..43fae19c 100644
--- a/demo/src/bunny.cc
+++ b/demo/src/bunny.cc
@@ -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);
diff --git a/demo/src/ellipt.cc b/demo/src/ellipt.cc
index b52698c1..a5ad31d1 100644
--- a/demo/src/ellipt.cc
+++ b/demo/src/ellipt.cc
@@ -1,4 +1,6 @@
#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");
+// applySingularPertubation(prob->getSystemMatrix(0,0)->getBaseMatrix());
+// applySingularDBC(prob->getSystemMatrix(0,0)->getBaseMatrix());
+
+ SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;
+ solverMatrix.setMatrix(*prob->getSystemMatrix());
+ prob->setSolverMatrix(solverMatrix);
+
+ debug::writeMatlabMatrix(*(prob->getSystemMatrix()), "ellipt_mat.m");
+ prob->solve(adaptInfo, false);
+ }
+
+ if (toDo.isSet(ESTIMATE)) {
+ prob->estimate(adaptInfo);
+ }
+ return flag;
+ };
+
+ template <typename Matrix>
+ void applySingularPertubation(Matrix& m)
+ {
+ using namespace mtl;
+ // Type of m's elements
+ typedef typename mtl::Collection<Matrix>::value_type value_type;
+ int nnz= m.nnz_local(m.num_rows()-1);
+ // Create inserter for matrix m
+ // Existing values are not overwritten but inserted
+ matrix::inserter<Matrix, update_times<value_type> > ins(m, nnz);
+
+ ins[m.num_rows()-1][m.num_rows()-1] << 1.0+1.e-5;
+ }
+
+ template <typename Matrix>
+ void applySingularDBC(Matrix& m)
+ {
+ using namespace mtl;
+
+ typename traits::row<Matrix>::type r(m);
+ typename traits::col<Matrix>::type c(m);
+ typename traits::value<Matrix>::type v(m);
+ typedef typename traits::range_generator<tag::row, Matrix>::type c_type;
+ typedef typename traits::range_generator<tag::nz, c_type>::type ic_type;
+// typedef typename Collection<Matrix>::value_type value_type;
+
+ for (c_type cursor(begin<tag::row>(m)), cend(end<tag::row>(m)); cursor != cend; ++cursor) {
+ for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
+ v(*icursor, (r(*icursor)==c(*icursor)?1.0:0.0));
+ }
+ break;
+ }
+ }
+
+ template <typename Matrix, typename Vector>
+ void applyDbc(Matrix& m, Vector& vec, DegreeOfFreedom idx, double value)
+ {
+ using namespace mtl;
+ typename Matrix::size_type idx_= idx;
+
+ typename traits::row<Matrix>::type r(m);
+ typename traits::col<Matrix>::type c(m);
+ typename traits::value<Matrix>::type v(m);
+ typedef typename traits::range_generator<tag::row, Matrix>::type c_type;
+ typedef typename traits::range_generator<tag::nz, c_type>::type ic_type;
+// typedef typename Collection<Matrix>::value_type value_type;
+
+ for (c_type cursor(begin<tag::row>(m)+idx_), cend(begin<tag::row>(m)+idx_+1); cursor != cend; ++cursor) {
+ for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
+ v(*icursor, (r(*icursor)==c(*icursor)?1.0:0.0));
+ if(r(*icursor)==c(*icursor)) {
+ MSG("DBC set to (%d, %d)\n", r(*icursor), c(*icursor));
+ }
+ }
+ break;
+ }
+
+ vec[idx]= value;
+ }
+
+
+DegreeOfFreedom getDOFidxAtCoords(DOFVector<double> *vec, WorldVector<double> &p, bool &inside, ElInfo *oldElInfo=NULL, bool useOldElInfo=false)
+{
+ const FiniteElemSpace *feSpace = vec->getFeSpace();
+ Mesh *mesh = vec->getFeSpace()->getMesh();
+ const BasisFunction *basFcts = feSpace->getBasisFcts();
+
+ int dim = mesh->getDim();
+ int numBasFcts = basFcts->getNumber();
+
+ DegreeOfFreedom *localIndices = new DegreeOfFreedom[numBasFcts];
+ DimVec<double> lambda(dim, NO_INIT);
+
+ ElInfo *elInfo = mesh->createNewElInfo();
+ int dofIdx=0;
+ inside=false;
+
+ if(oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
+ inside=mesh->findElInfoAtPoint(p,elInfo,lambda,oldElInfo->getMacroElement(),NULL,NULL);
+ delete oldElInfo;
+ } else {
+ inside=mesh->findElInfoAtPoint(p,elInfo,lambda,NULL,NULL,NULL);
+ }
+ if(oldElInfo) oldElInfo=elInfo;
+ if(inside) {
+ basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
+ FixVec<WorldVector<double>, VERTEX> coords = elInfo->getMacroElement()->getCoord();
+ int minIdx=-1;
+ double minDist=1.e15;
+ for(int i=0;i<coords.size();++i) {
+ WorldVector<double> dist = coords[i]-p;
+ double newDist = AMDiS::norm(dist);
+ if(newDist<minDist) {
+ minIdx=i;
+ minDist=newDist;
+ }
+ }
+ if(minIdx>=0) {
+ dofIdx = localIndices[minIdx];
+ }
+ }
+ delete [] localIndices;
+ if(!oldElInfo) delete elInfo;
+ return dofIdx;
+};
+protected:
+MyProblemVec *prob;
+};
+
+/// < factor*phase*grad(u) , grad(psi) >
+class Phase_SOT : public SecondOrderTerm
+{
+public:
+ /// Constructor.
+ Phase_SOT(DOFVectorBase<double>* phaseDV_, double fac_=1.0) :
+ SecondOrderTerm(6), phaseDV(phaseDV_), fac(fac_)
+ {
+ TEST_EXIT(phaseDV_)("no vector phase!\n");
+ auxFeSpaces.insert(phaseDV_->getFeSpace());
+ };
+
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL) {
+ getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
+ };
+
+ void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad = NULL) {
+ getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
+ };
+
+ inline void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const
+ {
+ const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
+ for (int iq = 0; iq < nPoints; iq++)
+ l1lt(Lambda, LALt[iq], f(iq) * phase[iq] * fac);
+ };
+
+ inline void eval(int nPoints,
+ const AMDiS::ElementVector&,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ ElementVector& result,
+ double opFactor)
+ {
+ if (D2UhAtQP) {
+ for (int iq = 0; iq < nPoints; iq++) {
+ double feval = f(iq) * phase[iq] * opFactor * fac;
+ double resultQP = 0.0;
+ for (int i = 0; i < dimOfWorld; i++)
+ resultQP += D2UhAtQP[iq][i][i];
+ result[iq] += feval * resultQP;
+ }
+ }
+ };
+
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result)
+ {
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = f(iq) * phase[iq] * fac;
+ axpy(factor, grdUhAtQP[iq], result[iq]);
+ }
+ };
+
+ virtual double f(const int iq) const
+ {
+ return 1.0;
+ };
+ void setFactor(const double fac_) { fac=fac_; }
+
+protected:
+ DOFVectorBase<double> *phaseDV;
+ mtl::dense_vector<double> phase;
+ double fac;
+};
+
+/// simple circle: phi_0(x,y) := sqrt(x^2+y^2)-r
+class Circle : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+ Circle(double radius_, WorldVector<double> midPoint_) :
+ radius(radius_), midPoint(midPoint_) {}
+ Circle(double radius_) :
+ radius(radius_) { midPoint.set(0.0); }
+ double operator()(const WorldVector<double>& x) const
+ {
+ double result = 0.0;
+
+ for(int k=0; k<x.getSize(); k++) {
+ result += sqr(x[k]-midPoint[k]);
+ }
+ result = sqrt(result)-radius;
+ return result;
+ };
+private:
+ double radius;
+ WorldVector<double> midPoint;
+};
+
/// Dirichlet boundary function
class G : public AbstractFunction<double, WorldVector<double> >
{
@@ -15,7 +295,7 @@ public:
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
- return exp(-10.0 * (x * x));
+ return 0.0;
}
};
@@ -29,10 +309,7 @@ public:
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
- int dow = Global::getGeo(WORLD);
- double r2 = (x * x);
- double ux = exp(-10.0 * r2);
- return -(400.0 * r2 - 20.0 * dow) * ux;
+ return sin(M_PI*x[0]);
}
};
@@ -53,22 +330,25 @@ int main(int argc, char* argv[])
// ===== create and init the scalar problem =====
- ProblemStat ellipt("ellipt");
+ MyProblemVec ellipt("ellipt");
ellipt.initialize(INIT_ALL);
-
+ MyIteration elIter(&ellipt);
// === create adapt info ===
- AdaptInfo adaptInfo("ellipt->adapt");
+ AdaptInfo adaptInfo("ellipt->adapt", 2);
// === create adapt ===
- AdaptStationary adapt("ellipt->adapt", ellipt, adaptInfo);
+ AdaptStationary adapt("ellipt->adapt", elIter, adaptInfo);
+// DOFVector<double> phase(ellipt.getFeSpace(), "phase");
+// phase.interpol(new Circle(0.25));
// ===== create matrix operator =====
Operator matrixOperator(ellipt.getFeSpace());
matrixOperator.addSecondOrderTerm(new Simple_SOT);
ellipt.addMatrixOperator(matrixOperator, 0, 0);
+ ellipt.addMatrixOperator(matrixOperator, 1, 1);
// ===== create rhs operator =====
@@ -76,10 +356,12 @@ int main(int argc, char* argv[])
Operator rhsOperator(ellipt.getFeSpace());
rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
ellipt.addVectorOperator(rhsOperator, 0);
+ ellipt.addVectorOperator(rhsOperator, 1);
// ===== add boundary conditions =====
- ellipt.addDirichletBC(1, 0, 0, new G);
+ ellipt.addDirichletBC(3, 0, 0, new G);
+ ellipt.addDirichletBC(5, 1, 1, new G);
// ===== start adaption loop =====
--
GitLab