Commit e2b785f0 authored by Praetorius, Simon's avatar Praetorius, Simon

PeriodicBC corrected, GeometryTools moved to new namespace, BaseProblem extended

parent 544ddc37
......@@ -37,17 +37,35 @@ public:
for (int i = 0; i < nComponents; ++i)
exactSolutions[i] = NULL;
}
template<typename SubProblemType>
ExtendedProblemStat(std::string nameStr, SubProblemType *subProblem)
:
ProblemStat_(nameStr, subProblem)
, oldMeshChangeIdx(0)
{
exactSolutions.resize(nComponents);
for (int i = 0; i < nComponents; ++i)
exactSolutions[i] = NULL;
}
~ExtendedProblemStat()
{
for (int i = 0; i < nComponents; ++i)
if (exactSolutions[i]) {
delete exactSolutions[i];
exactSolutions[i] = NULL;
}
}
void initialize(Flag initFlag,
ProblemStatSeq *adoptProblem = NULL,
Flag adoptFlag = INIT_NOTHING)
{
ProblemStat_::initialize(initFlag, adoptProblem, adoptFlag);
// if (initFlag.isSet(INIT_EXACT_SOLUTION)) {
for (int i = 0; i < nComponents; ++i)
exactSolutions[i] = new DOFVector< double >(getFeSpace(i), "exact_solution");
// }
for (int i = 0; i < nComponents; ++i)
exactSolutions[i] = new DOFVector< double >(getFeSpace(i), "exact_solution");
}
inline DOFVector<double>* getRhsVector(int i = 0)
......@@ -85,7 +103,7 @@ public:
manualPeriodicBC.clear();
std::vector<PeriodicBcData>::iterator it;
for (it = PeriodicBcDataList.begin(); it != PeriodicBcDataList.end(); it++)
it->addToList(getFeSpace(), manualPeriodicBC);
it->addToList(getFeSpace(it->row), manualPeriodicBC);
}
// update dirichlet data
......
......@@ -2,7 +2,7 @@
#define EPSILON FLT_EPSILON
namespace meshconv {
namespace meshconv2 {
// ###############
// ## General ##
......
......@@ -5,7 +5,7 @@
#include "FixVec.h"
namespace meshconv {
namespace meshconv2 {
/**
* \defgroup General Dimension independet methods
......@@ -227,7 +227,7 @@ void coordinate_transform(double x[], double shift[], double scale[], double alp
double volume_tetrahedron(double tet0[], double tet1[], double tet2[], double tet3[]);
/**@}*/
} // end namespace meshconv
} // end namespace meshconv2
#include "GeometryTools.hh"
......
namespace meshconv {
namespace meshconv2 {
template<int dow>
void coordinate_transform(double x[], double shift[], double scale[], double alpha, double beta)
......@@ -38,4 +38,4 @@ namespace meshconv {
x[i] = coords[i]+shift[i];
}
} // end namespace meshconv
} // end namespace meshconv2
......@@ -87,7 +87,7 @@ void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> sh
T_ = Rz*Ry*Rx*S;
WorldMatrix<double> T;
vector_operations::fillMatrix(T,T_);
vector_operations::fillMatrix(T_,T);
FixVec<WorldVector<double>, VERTEX> coords;
deque<MacroElement*>::iterator macro;
......
......@@ -179,7 +179,7 @@ struct PhaseFieldToSignedDist : AbstractFunction<double, double>
double operator()(const double &phi) const
{
double z = std::max(-1.0 + 1.e-10, std::min(1.0 - 1.e-10, -phi));
double z = std::max(-1.0 + 1.e-12, std::min(1.0 - 1.e-12, -phi));
return (epsilon/scalingFactor) * atanh(z);
}
......
......@@ -24,7 +24,7 @@ using namespace AMDiS;
* - PlaneRotation(shift, angle)
* - Propeller(radius, angle)
* - Torus(radius1, radius2)
* - Polygon(v0, v1, v2, v2), or Polygon(vec(v_i))
* - Polygon(v0, v1, v2, v3), or Polygon(vec(v_i))
**/
// signed distance functions for surfaces in polar coords
......@@ -32,7 +32,7 @@ static double signedDist2D(const WorldVector<double> x, const WorldVector<double
AbstractFunction<double, double> *radius, double eps)
{
FUNCNAME("signedDist2D");
using vector_operations::epsNorm;
using vector_operations::norm;
WorldVector<double> x_trans;
double norm_xy;
......@@ -45,7 +45,7 @@ static double signedDist2D(const WorldVector<double> x, const WorldVector<double
x_trans[k] = x[k] - midPoint[k];
}
norm_xy = epsNorm(x_trans, eps);
norm_xy = norm(x_trans, eps);
if (x_trans[0] >= 0.0) {
alpha = acos(x_trans[0] / norm_xy);
......@@ -60,7 +60,7 @@ static double signedDist3D(const WorldVector<double> x, const WorldVector<double
BinaryAbstractFunction<double, double, double> *radius, double eps)
{
FUNCNAME("signedDist3D");
using vector_operations::epsNorm;
using vector_operations::norm;
WorldVector<double> x_trans;
double norm_xyz, norm_xy;
......@@ -73,7 +73,7 @@ static double signedDist3D(const WorldVector<double> x, const WorldVector<double
x_trans[k] = x[k] - midPoint[k];
}
norm_xyz = epsNorm(x_trans, eps);
norm_xyz = norm(x_trans, eps);
norm_xy = sqrt(sqr(x_trans[0]) + sqr(x_trans[1]) + sqr(eps));
if(x_trans[1]>=0) {
......@@ -568,8 +568,8 @@ public:
{
double result = 1.e15;
for (size_t i = 0; i < vertices.size()-1; i++)
result = std::min(result, meshconv::distance_point_line_2d(x.begin(), vertices[i].begin(), vertices[i+1].begin()));
return factor * result * (meshconv::point_in_polygon(x.begin(), vertices) ? -1.0 : 1.0);
result = std::min(result, meshconv2::distance_point_line_2d(x.begin(), vertices[i].begin(), vertices[i+1].begin()));
return factor * result * (meshconv2::point_in_polygon(x.begin(), vertices) ? -1.0 : 1.0);
}
void refine(int np)
......
......@@ -39,7 +39,7 @@ namespace details {
inline void getBoundaryIndices(const FiniteElemSpace* feSpace,
BoundaryType nr,
AbstractFunction<bool, WorldVector<double> >* meshIndicator,
std::vector<DegreeOfFreedom> &indices)
std::set<DegreeOfFreedom> &indices)
{
Mesh* mesh = feSpace->getMesh();
int dim = mesh->getDim();
......@@ -61,7 +61,7 @@ namespace details {
for (int i = 0; i < numBasFcts; i++) {
elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
if ((*meshIndicator)(coord))
indices.push_back(localIndices[i]);
indices.insert(localIndices[i]);
}
break;
}
......@@ -77,7 +77,7 @@ namespace details {
**/
inline void getBoundaryIndices(const FiniteElemSpace* feSpace,
AbstractFunction<bool, WorldVector<double> >* meshIndicator,
std::vector<DegreeOfFreedom> &indices)
std::set<DegreeOfFreedom> &indices)
{
Mesh* mesh = feSpace->getMesh();
int dim = mesh->getDim();
......@@ -99,7 +99,7 @@ namespace details {
for (int i = 0; i < numBasFcts; i++) {
elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
if ((*meshIndicator)(coord))
indices.push_back(localIndices[i]);
indices.insert(localIndices[i]);
}
break;
}
......@@ -191,23 +191,82 @@ namespace details {
AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
{
std::vector<DegreeOfFreedom> indices;
std::set<DegreeOfFreedom> indices;
details::getBoundaryIndices(feSpace, nr, meshIndicator, indices);
association.clear();
DOFVector<WorldVector<double> > coords(feSpace, "coords");
feSpace->getMesh()->getDofIndexCoords(coords);
for (size_t i = 0; i < indices.size(); i++)
std::set<DegreeOfFreedom>::iterator it;
for (it = indices.begin(); it != indices.end(); it++)
{
DegreeOfFreedom idx2;
WorldVector<double> p;
p = (*periodicMap)(coords[indices[i]]);
p = (*periodicMap)(coords[*it]);
TEST_EXIT(coords.getDofIdxAtPoint(p,idx2))("periodic association not found!\n");
association.push_back(std::make_pair(indices[i], idx2));
association.push_back(std::make_pair(*it, idx2));
}
}
inline bool getDofIdxAtPoint(const FiniteElemSpace* feSpace,
WorldVector<double> &p,
DegreeOfFreedom &idx,
ElInfo *oldElInfo = NULL,
bool useOldElInfo = false)
{
FUNCNAME("DOFVector::getDofIdxAtPoint()");
Mesh *mesh = feSpace->getMesh();
const BasisFunction *basFcts = feSpace->getBasisFcts();
int dim = mesh->getDim();
int numBasFcts = basFcts->getNumber();
std::vector<DegreeOfFreedom> localIndices(numBasFcts);
DimVec<double> lambda(dim, NO_INIT);
ElInfo *elInfo = mesh->createNewElInfo();
idx = 0;
bool 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)
return false;
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
WorldVector<double> coord;
int minIdx = -1;
double minDist = 1.e15;
for (int i = 0; i < numBasFcts; i++) {
elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
WorldVector<double> dist = coord - p;
double newDist = norm(dist);
if (newDist < minDist) {
minIdx = i;
minDist = newDist;
}
}
if (minIdx >= 0)
idx = localIndices[minIdx];
if(!oldElInfo) delete elInfo;
return inside;
}
/**
* get al list of pairs that describes wich DegreeOfFreedom are assiciated in the given mesh
......@@ -216,32 +275,34 @@ namespace details {
* the 'periodicMap' associates a second DegreeOfFreedom.
**/
inline void getPeriodicAssociation(const FiniteElemSpace* feSpace,
AbstractFunction<bool, WorldVector<double> >* meshIndicator,
AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
AbstractFunction<bool, WorldVector<double> >* meshIndicator,
AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
{
std::vector<DegreeOfFreedom> indices;
std::set<DegreeOfFreedom> indices;
details::getBoundaryIndices(feSpace, meshIndicator, indices);
MSG_DBG("Nr of indices: %d\n", indices.size());
association.clear();
DOFVector<WorldVector<double> > coords(feSpace, "coords");
feSpace->getMesh()->getDofIndexCoords(coords);
for (size_t i = 0; i < indices.size(); i++)
std::set<DegreeOfFreedom>::iterator it;
for (it = indices.begin(); it != indices.end(); it++)
{
DegreeOfFreedom idx2;
WorldVector<double> p;
p = (*periodicMap)(coords[indices[i]]);
p = (*periodicMap)(coords[*it]);
TEST_EXIT(coords.getDofIdxAtPoint(p,idx2))("periodic association not found!\n");
association.push_back(std::make_pair(indices[i], idx2));
association.push_back(std::make_pair(*it, idx2));
}
MSG_DBG("Nr of associations: %d\n", association.size());
}
inline void getDOFValues(const FiniteElemSpace* feSpace,
AbstractFunction<double, WorldVector<double> >* values,
std::vector<DegreeOfFreedom> &indices,
std::vector<double> &dofValues)
inline void getDOFValues(const FiniteElemSpace* feSpace, AbstractFunction<double, WorldVector<double> >* values, std::vector<DegreeOfFreedom> &indices, std::vector<double> &dofValues)
{
WorldVector<double> x;
for (size_t i = 0; i < indices.size(); i++) {
......@@ -250,13 +311,11 @@ namespace details {
}
}
inline void getDOFValues(const FiniteElemSpace* feSpace,
DOFVector<double>* values,
std::vector<DegreeOfFreedom> &indices,
std::vector<double> &dofValues)
inline void getDOFValues(const FiniteElemSpace* feSpace, DOFVector<double>* values, std::vector<DegreeOfFreedom> &indices, std::vector<double> &dofValues)
{
for (size_t i = 0; i < indices.size(); i++)
for (size_t i = 0; i < indices.size(); i++) {
dofValues.push_back((*values)[indices[i]]);
}
}
} // end namespace
......@@ -416,10 +475,14 @@ struct DirichletBcData {
break;
case 4:
std::set<DegreeOfFreedom> indices__;
if (boundary_nr == 0)
details:: getBoundaryIndices(feSpace, meshIndicator, indices_);
details:: getBoundaryIndices(feSpace, meshIndicator, indices__);
else
details:: getBoundaryIndices(feSpace, boundary_nr, meshIndicator, indices_);
details:: getBoundaryIndices(feSpace, boundary_nr, meshIndicator, indices__);
std::set<DegreeOfFreedom>::iterator it;
for (it = indices__.begin(); it != indices__.end(); it++)
indices_.push_back(*it);
break;
}
......@@ -441,8 +504,9 @@ struct DirichletBcData {
feSpace->getMesh()->getDofIndexCoords(coords);
for (size_t i = 0; i < indices_.size(); i++)
for (size_t i = 0; i < indices_.size(); i++) {
list.push_back(SingularDirichletBC(row, col, indices_[i], values_[i]));
}
MSG_DBG("Dirichle BC at %d DOFs added.\n",indices_.size());
}
......@@ -490,10 +554,10 @@ struct PeriodicBcData {
details::getPeriodicAssociation(feSpace, nr, meshIndicator, periodicMap, associations);
list.push_back(ManualPeriodicBC(row, associations));
}
int row;
protected:
int row;
BoundaryType nr;
AbstractFunction<bool, WorldVector<double> >* meshIndicator;
......
......@@ -117,7 +117,7 @@ namespace tools {
struct SequenceQuad
{
// integration zwischen zwei Punkten mittels Trapezregel
template<typename Vector1, typename Vector2, typename T> /// requires Callable1<Functor, T>, Addable<T>, Multiplicable<T>, Nullifyable<T>
template<typename Vector1, typename Vector2, typename T> /// requires Addable<T>, Multiplicable<T>, Nullifyable<T>
T quad(Vector1& y, Vector2& x)
{
T I; nullify(I);
......
......@@ -18,6 +18,7 @@
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp>
#include <boost/utility/enable_if.hpp>
#include <boost/type_traits.hpp>
using namespace AMDiS;
......@@ -108,6 +109,7 @@ namespace vector_operations {
type operator()(const Value[Rows][Cols]) { return Rows; }
};
/// size implementation for AMDiS WorldVectors
template <typename Value>
struct num_rows< WorldVector<Value> >
{
......@@ -115,6 +117,7 @@ namespace vector_operations {
type operator()(const WorldVector<Value>& v) { return static_cast<size_t>(v.getSize()); }
};
/// size implementation for AMDiS WorldMatrices
template <typename Value>
struct num_rows< WorldMatrix<Value> >
{
......@@ -122,6 +125,7 @@ namespace vector_operations {
type operator()(const WorldMatrix<Value>& v) { return static_cast<size_t>(v.getNumRows()); }
};
/// size implementation for arithmetic types (double, float, int, ...)
template <typename Value>
struct num_rows< Value, typename boost::enable_if< boost::is_arithmetic<Value> >::type >
{
......@@ -135,6 +139,7 @@ namespace vector_operations {
template <typename Collection, class Enable = void>
struct num_cols {};
/// size implementation for AMDiS WorldMatrices
template <typename Value>
struct num_cols< WorldMatrix<Value> >
{
......@@ -148,6 +153,7 @@ namespace vector_operations {
template <typename Collection, class Enable = void>
struct resize {};
/// change_dim implementation for AMDiS WorldVectors
template <typename Value>
struct resize< WorldVector<Value> >
{
......@@ -158,6 +164,7 @@ namespace vector_operations {
}
};
/// change_dim implementation for STL vectors
template <typename Value>
struct resize< std::vector<Value> >
{
......@@ -167,6 +174,7 @@ namespace vector_operations {
}
};
/// change_dim implementation for MTL4 vectors
template <typename Value>
struct resize< mtl::dense_vector<Value> >
{
......@@ -178,6 +186,7 @@ namespace vector_operations {
// _________________________________________________________________________________
/// change_dim implementation for AMDiS WorldMatrices
template <typename Value>
struct resize< WorldMatrix<Value> >
{
......@@ -188,6 +197,7 @@ namespace vector_operations {
}
};
/// change_dim implementation for MTL4 matrix
template <typename Value, typename Param>
struct resize< mtl::matrix::base_matrix<Value, Param> >
{
......@@ -197,6 +207,7 @@ namespace vector_operations {
}
};
/// change_dim implementation for MTL4 matrix
template <typename Value>
struct resize< mtl::matrix::dense2D<Value> >
{
......@@ -206,6 +217,7 @@ namespace vector_operations {
}
};
/// change_dim implementation for MTL4 matrix
template <typename Value>
struct resize< mtl::matrix::compressed2D<Value> >
{
......@@ -214,12 +226,67 @@ namespace vector_operations {
v.change_dim(r,c);
}
};
// _________________________________________________________________________________
/// General declaration, used to disable unsupported types
template <typename Collection, class Enable = void>
struct at {};
template <typename Value>
struct at< WorldMatrix<Value> >
{
typedef Value type;
type& operator()(WorldMatrix<Value>& v, size_t r, size_t c) {
return v[r][c];
}
type const& operator()(const WorldMatrix<Value>& v, size_t r, size_t c) {
return v[r][c];
}
};
template <typename Value, typename Param>
struct at< mtl::matrix::base_matrix<Value, Param> >
{
typedef Value type;
type& operator()(mtl::matrix::base_matrix<Value, Param>& v, size_t r, size_t c) {
return v(r,c);
}
type const& operator()(const mtl::matrix::base_matrix<Value, Param>& v, size_t r, size_t c) {
return v(r,c);
}
};
template <typename Value>
struct at< mtl::matrix::dense2D<Value> >
{
typedef Value type;
type& operator()(mtl::matrix::dense2D<Value>& v, size_t r, size_t c) {
return v(r,c);
}
type const& operator()(const mtl::matrix::dense2D<Value>& v, size_t r, size_t c) {
return v(r,c);
}
};
template <typename Value>
struct at< mtl::matrix::compressed2D<Value> >
{
typedef Value type;
type& operator()(mtl::matrix::compressed2D<Value>& v, size_t r, size_t c) {
return v(r,c);
}
type const& operator()(const mtl::matrix::compressed2D<Value>& v, size_t r, size_t c) {
return v(r,c);
}
};
}
/// num_rows function for non-MTL types (uses implicit enable_if)
template <typename Collection>
typename traits::num_rows<Collection>::type
inline num_rows(const Collection& c)
num_rows(const Collection& c)
{
return traits::num_rows<Collection>()(c);
}
......@@ -227,41 +294,63 @@ namespace vector_operations {
/// num_cols function for non-MTL types (uses implicit enable_if)
template <typename Collection>
typename traits::num_cols<Collection>::type
inline num_cols(const Collection& c)
num_cols(const Collection& c)
{
return traits::num_cols<Collection>()(c);
}
/// resize function for vectors
template <typename Collection>
typename traits::num_cols<Collection>::vector_void_type
inline resize(const Collection& c, size_t rows)
typename traits::resize<Collection>::vector_void_type
resize(Collection& c, size_t rows)
{
traits::resize<Collection>()(c, rows);
}
/// resize function for matrices
template <typename Collection>
typename traits::num_cols<Collection>::matrix_void_type
inline resize(const Collection& c, size_t rows, size_t cols)
typename traits::resize<Collection>::matrix_void_type
resize(Collection& c, size_t rows, size_t cols)
{
traits::resize<Collection>()(c, rows, cols);
}
template<typename T>
T norm(mtl::dense_vector<T> &b) { return two_norm(b); };
template<typename T>
T norm(std::vector<T> &b)
/// at function for matrices
template <typename Collection>
typename traits::at<Collection>::type const&
at(const Collection& c, size_t rows, size_t cols)
{
return traits::at<Collection>()(c, rows, cols);
}
/// at function for matrices
template <typename Collection>
typename traits::at<Collection>::type&
at(Collection& c, size_t rows, size_t cols)
{
return traits::at<Collection>()(c, rows, cols);
}
/// 2-norm for vectors
template<typename Vector>
typename boost::enable_if< is_vector<Vector>, typename ValueType<Vector>::type >::type
norm(const Vector &b)
{
T value = 0.0;
for (size_t i = 0; i < b.size(); ++i)
typename ValueType<Vector>::type value;
nullify(value);
for (size_t i = 0; i < num_rows(b); ++i)
value += sqr(b[i]);
return sqrt(value);
};
/// 2-norm for MTL4 vectors
template<typename T>
T norm(const mtl::dense_vector<T> &b) { return two_norm(b); };
/// 1-norm for matrices
template<typename Matrix>
double matrix_norm(Matrix &M)
typename boost::enable_if< is_matrix<Matrix>, double >::type
matrix_norm(const Matrix &M)
{