Skip to content
GitLab
Commit 2c1e3427 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

NULL replaced by nullptr and a lot of small changes

parent 02ff3473
Hide whitespace changes
Inline Side-by-side
extensions/BackgroundMesh.h
View file @ 2c1e3427
......@@ -53,6 +53,20 @@ namespace experimental {
typedef std::pair<DegreeOfFreedom, PointType> DataType;
/** \brief Background-mesh structure
* Can be used to find k nearest neighbors to a given point. It uses a one
* level background-mesh to look for a box where the point lies inside and
* compares the distance to all other points in the box, respective with the
* points in the surrounding boxes.
* The point to look for are stores as DataType-objects, e.g. pairs of
* DOF-indices and coordinates. This makes it possible to evaluated DOFVectors
* at the found points.
*
* To use this class call Box::provideBackgroundMesh(feSpace) with feSpace
* a FiniteElemSpace holding all the degrees of freedom you want to look for.
* Or add an individual fillBox method that adds all the data points in the
* form of DataType to the box by calling addData for each point.
**/
struct Box {
Box(int DOW_, std::vector<int> N_);
......
extensions/BackgroundMesh2.h
View file @ 2c1e3427
......@@ -22,7 +22,9 @@
#include "AMDiS.h"
#include "Tools.h"
#ifndef USE_EXPERIMENTAL
#define USE_EXPERIMENTAL
#endif
namespace experimental {
......@@ -30,9 +32,34 @@ typedef WorldVector<double> PointType_;
typedef boost::tuple<const MacroElement*, int, unsigned long> ElInfoDataType; // data to recover elInfo
typedef std::pair<ElInfoDataType, PointType_> DataType_;
template<typename PointType, typename DataType>
// PointType ... type of coordinate vectors, e.g. WorldVector<double> in AMDiS
// DataType .... use DataType_ to look for ElInfos near coordinates
/** \brief Background-mesh structure for finding elements
* Can be used to find k nearest elements to a given point. It uses a one
* level background-mesh to look for a box where the point lies inside and
* compares the distance to all element centers in the box, respective with the
* elements in the surrounding boxes.
* The elements to look for are stores as DataType-objects, e.g. tuples of
* MacroElement*, refinementPath and refinementPathLength. This makes it
* possible to create a new ElInfo object from the mesh.
*
* To use this class call ElementBox<P,D>::provideBackgroundMesh(feSpace) with
* feSpace a FiniteElemSpace holding the mesh with all elements you want to look for.
* The template-arguments correspond to coordinate-vector type P and D a type
* in wich the elInfo is stored. Use PointType_ and DataType_ as defaults.
* Or add an individual fillBox method that adds all the elements in the
* form of DataType to the box by calling addData for each element.
**/
template<typename PointType=PointType_, typename DataType=DataType_>
struct ElementBox {
/** \brief Constructor for BackgroundMesh
* Mesh is enclosed in default box [-1,1]^2.
*
* \param DOW_ dimension of world
* \param N_ vector containing number of boxes per dimension [Nx,Ny,Nz]
**/
ElementBox(int DOW_, std::vector<int> N_)
: DOW(DOW_), N(N_), boxFilled(false)
{
......@@ -41,6 +68,14 @@ struct ElementBox {
init();
}
/** \brief Constructor for BackgroundMesh
* Mesh is enclosed in a box given by the user.
*
* \param DOW_ dimension of world
* \param min_corner_ lower-left[-front] corner of enclosing box
* \param max_corner_ upper-right[-back] corner of enclosing box
* \param N_ vector containing number of boxes per dimension [Nx,Ny,Nz]
**/
ElementBox(int DOW_, PointType min_corner_, PointType max_corner_, std::vector<int> N_)
: DOW(DOW_),
min_corner(min_corner_),
......@@ -50,15 +85,19 @@ struct ElementBox {
init();
}
/// add all elements of mesh corresponding to the feSpace to the background-mesh
void fillBox(const FiniteElemSpace* feSpace);
/// delete and clear all boxes
void clearBox()
{
boxData.clear();
}
/// calculate boxSize
int getMaxBoxSize()
{
size_t maxSize = 0;
......@@ -69,6 +108,7 @@ struct ElementBox {
}
/// test of point x lies in current box
bool inBox(PointType& x)
{
for (size_t i = 0; i < DOW; i++)
......@@ -96,9 +136,14 @@ struct ElementBox {
}
/// return ElInfo* created from DataType data nearest to x.
/// Template specialization for DataType=DataType_ and
/// PointType=PointType_ provided.
bool getNearestElInfo(PointType &x, ElInfo*& minElInfo);
/// get object of type DataType nearest to given coordinate x,
/// where data can be used to create an ElInfo from the mesh.
bool getNearestData(PointType& x, DataType& data)
{
if (!inBox(x))
......@@ -123,6 +168,9 @@ struct ElementBox {
return true;
}
/// get nData objects of type DataType nearest to given coordinate x,
/// where data can be used to create an ElInfo from the mesh.
bool getNearestData(PointType& x, std::vector<DataType>& data, int nData)
{
if (!inBox(x))
......@@ -190,6 +238,14 @@ struct ElementBox {
}
/** \brief create background-mesh structure for given feSpace
* Using fillBox() implementation to store element-information in the DataType_
* format in the background-mesh structure. When a new background mesh is created
* parameters about size are read from init-file:
* - <c>backgroundMesh->N</c> reads to \ref N
* - <c>backgroundMesh->min_corner</c> reads to \ref min_corner
* - <c>backgroundMesh->max_corner</c> reads to \ref max_corner
**/
static ElementBox<PointType, DataType>* provideBackgroundMesh(const FiniteElemSpace* feSpace)
{
if (boxMap.find(feSpace) != boxMap.end()) {
......@@ -218,6 +274,7 @@ struct ElementBox {
return boxMap[feSpace].second;
}
/// clear background mesh
static void delete_all()
{
typename std::map<const FiniteElemSpace*, std::pair<int, ElementBox<PointType, DataType>*> >::iterator boxIter;
......@@ -261,7 +318,7 @@ protected:
return nr;
}
inline void nr2idx(int nr, std::vector<int>& idx) // template specialization see below
inline void nr2idx(int nr, std::vector<int>& idx)
{
switch (DOW) {
case 1:
......@@ -286,7 +343,7 @@ protected:
}
}
/// get indices of boxes surrounding the box with given nr
void getSurroundingBoxes(int nr, std::set<int> &surrounding_nrs)
{
std::vector<int> idx;
......@@ -302,6 +359,8 @@ protected:
surrounding_nrs.erase(nr);
}
/// get indices of boxes surrounding the boxes with given set of nrs
void getSurroundingBoxes(std::set<int> &nrs, std::set<int> &surrounding_nrs)
{
std::set<int>::iterator nrsIter;
......@@ -312,7 +371,7 @@ protected:
surrounding_nrs.erase(*nrsIter);
}
/// berechne die minimale Entfernung der Ränder der Box zu x
/// get minimal distance of box-boundaries to x
double getBoxBoundaryDist(int nr, PointType& x)
{
std::vector<int> idx;
......@@ -324,10 +383,11 @@ protected:
}
double dist = 1.e15;
for (int i = 0; i < DOW; i++)
dist = std::min(dist, std::min(abs(max_c[i]-x[i]), abs(min_c[i]-x[i])));
dist = std::min(dist, std::min(std::abs(max_c[i]-x[i]), std::abs(min_c[i]-x[i])));
return dist;
}
/// add data to the background-mesh at point x
void addData(PointType x, DataType data)
{
int nr = getBox(x);
......@@ -336,7 +396,7 @@ protected:
boxData[nr].push_back(data);
}
/// used by \ref provideBackgroundMesh
static typename std::map<const FiniteElemSpace*, std::pair<int, ElementBox<PointType, DataType>*> > boxMap;
protected:
......
extensions/ExtendedProblemStat.h
View file @ 2c1e3427
......@@ -20,7 +20,7 @@
#define EXTENSIONS_EXTENDED_PROBLEM_STAT_H
#include "AMDiS.h"
#include "SingularDirichletBC.h"
#include "SingularDirichletBC2.h"
#if defined NONLIN_PROBLEM
#include "nonlin/ProblemNonLin.h"
#endif
......@@ -116,7 +116,6 @@ public:
{
ProblemStat_::buildAfterCoarsen(adaptInfo, flag, asmMatrix, asmVector);
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
// update periodic data
if (oldMeshChangeIdx != getMesh()->getChangeIndex()
|| flag.isSet(UPDATE_PERIODIC_BC)
......@@ -130,16 +129,15 @@ public:
// apply periodic boundary conditions
for (size_t k = 0; k < manualPeriodicBC.size(); k++)
applyPeriodicBC(manualPeriodicBC[k], asmMatrix, asmVector);
#endif
// update dirichlet data
if (oldMeshChangeIdx != getMesh()->getChangeIndex()
|| flag.isSet(UPDATE_DIRICHLET_BC)
|| (DirichletBcDataList.size() > 0 && singularDirichletBC.size() == 0)) {
singularDirichletBC.clear();
std::vector<DirichletBcData>::iterator it;
std::vector<DirichletBcDataBase*>::iterator it;
for (it = DirichletBcDataList.begin(); it != DirichletBcDataList.end(); it++) {
it->addToList(getFeSpace(it->row), singularDirichletBC);
(*it)->addToList(getFeSpace((*it)->row), singularDirichletBC);
}
}
......@@ -176,6 +174,29 @@ public:
}
//============================================================================
template<typename ValueContainer>
void addDirichletBC(BoundaryType type, int row, int col,
ValueContainer *values)
{
BoundaryTypeContainer *bound = new BoundaryTypeContainer(type);
DirichletBcDataList.push_back(
new DirichletBcData<BoundaryTypeContainer, ValueContainer>(
row, col, *bound, *values));
}
// general templatized version of addDirichletBC
template<typename PosType, typename ValueContainer>
void addDirichletBC(PosType *pos, int row, int col,
ValueContainer *values)
{
DirichletBcDataList.push_back(
new DirichletBcData<PosType, ValueContainer>(row, col, *pos, *values));
}
/**
* Dirichlet boundary condition at DOF-Index of vertex near to coords 'pos'.
* Value defined by AbstractFunction, that is evaluated at 'pos'
......@@ -184,7 +205,8 @@ public:
void addSingularDirichletBC(WorldVector<double> &pos, int row, int col,
ValueContainer &values)
{
DirichletBcDataList.push_back(DirichletBcData(row, col, pos, values));
WARNING("Depreciated! Use addDirichletBC instead!!!\n");
addDirichletBC(&pos, row, col, &values);
}
......@@ -192,13 +214,23 @@ public:
* Dirichlet boundary condition at DOF-Index 'idx'. Value defined by double.
**/
template<typename ValueContainer>
void addSingularDirichletBC(DegreeOfFreedom idx, int row, int col,
ValueContainer &values)
void addDirichletBC(DofIndex* idx, int row, int col,
ValueContainer *values)
{
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
if (MPI::COMM_WORLD.Get_rank() == 0)
#endif
DirichletBcDataList.push_back(DirichletBcData(row, col, idx, values));
DirichletBcDataList.push_back(
new DirichletBcData<DofIndex, ValueContainer>(row, col, *idx, *values));
}
template<typename ValueContainer>
void addSingularDirichletBC(DegreeOfFreedom idx, int row, int col,
ValueContainer &values)
{
WARNING("Depreciated! Use addDirichletBC instead!!!\n");
DofIndex* dofIndex = new DofIndex(idx);
addDirichletBC(*dofIndex, row, col, values);
}
......@@ -213,37 +245,57 @@ public:
int row, int col,
ValueContainer &values)
{
DirichletBcDataList.push_back(DirichletBcData(row, col, &signedDist, values));
WARNING("Depreciated! Use addDirichletBC instead!!!\n");
addDirichletBC(&signedDist, row, col, &values);
}
///
template<typename ValueContainer>
void addImplicitDirichletBC(AbstractFunction<double, WorldVector<double> > &signedDist,
int row, int col,
ValueContainer &values)
{
DirichletBcDataList.push_back(DirichletBcData(row, col, &signedDist, values));
WARNING("Depreciated! Use addDirichletBC instead!!!\n");
addDirichletBC(&signedDist, row, col, &values);
}
///
template<typename ValueContainer>
void addDirichletBC(BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
int row, int col,
ValueContainer *values)
{
DirichletBcDataList.push_back(
new DirichletBcData<AbstractFunction<bool, WorldVector<double> >, ValueContainer>(
row, col, nr, *meshIndicator, *values));
}
template<typename ValueContainer>
void addManualDirichletBC(BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
int row, int col,
ValueContainer &values)
{
DirichletBcDataList.push_back(DirichletBcData(row, col, nr, meshIndicator, values));
WARNING("Depreciated! Use addDirichletBC instead!!!\n");
addDirichletBC(nr, meshIndicator, row, col, &values);
}
///
template<typename ValueContainer>
void addManualDirichletBC(AbstractFunction<bool, WorldVector<double> >* meshIndicator,
int row, int col,
ValueContainer &values)
{
DirichletBcDataList.push_back(DirichletBcData(row, col, meshIndicator, values));
WARNING("Depreciated! Use addDirichletBC instead!!!\n");
addDirichletBC(meshIndicator, row, col, &values);
}
// ===========================================================================
void addManualPeriodicBC(int row,
BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap)
......@@ -292,11 +344,19 @@ protected:
if (asmMatrix) {
typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
typedef mtl::traits::range_generator<tag::nz, c_type>::type ic_type;
// if matrix-block for the identity row is NULL, create new DOFMatrix
if (getSystemMatrix(row_, col_) == NULL) {
TEST_EXIT(row_ != col_)("should have been created already\n");
(*systemMatrix)[row_][col_] = new DOFMatrix(getFeSpace(row_), getFeSpace(col_), "");
(*systemMatrix)[row_][col_]->setCoupleMatrix(true);
(*systemMatrix)[row_][col_]->getBoundaryManager()->
setBoundaryConditionMap((*systemMatrix)[row_][row_]->getBoundaryManager()->
getBoundaryConditionMap());
}
size_t numCols = static_cast<size_t>(getNumComponents());
for (size_t col = 0; col < numCols; col++) {
TEST_EXIT(getSystemMatrix(row_, col) != NULL || col != col_)
("SystemMatrix block (%d,%d) must not be NULL! Insert a Simple_ZOT(0.0) as Workaround.\n",row_,col);
if (getSystemMatrix(row_, col) == NULL)
continue;
......@@ -463,7 +523,7 @@ private:
// data for dirichlet boundary conditions
std::vector<SingularDirichletBC> singularDirichletBC;
std::vector<DirichletBcData> DirichletBcDataList;
std::vector<DirichletBcDataBase*> DirichletBcDataList;
std::map<const FiniteElemSpace*, bool> feSpaceVisited;
......
extensions/Refinement.h
View file @ 2c1e3427
......@@ -149,11 +149,11 @@ public:
markElements(markFlag);
meshChanged = refineMesh(markFlag, onlyRefine);
calcMeshSizes(minH, maxH, minLevel, maxLevel);
calcMeshSizes(minH, maxH, minLevel, maxLevel, false);
int nr = mesh->getNumberOfVertices();
meshChanged = meshChanged && oldOldNr!=nr && oldNr!=nr;
if (meshChanged) {
MSG("Mesh sizes: [%f, %f], Vs: %d, ELs: %d\n",
MSG("(local) mesh sizes: [%f, %f], Vs: %d, ELs: %d\n",
minH, maxH, nr, mesh->getNumberOfElements());
}
i++;
......@@ -161,8 +161,8 @@ public:
oldOldNr = oldNr;
oldNr = nr;
}
calcMeshSizes(minH, maxH, minLevel, maxLevel);
MSG("Final mesh: [%f, %f], Vs: %d, ELs: %d, Level: [%d, %d]\n",
calcMeshSizes(minH, maxH, minLevel, maxLevel, true);
MSG("Final (global) mesh: [%f, %f], Vs: %d, ELs: %d, Level: [%d, %d]\n",
minH, maxH, mesh->getNumberOfVertices(), mesh->getNumberOfElements(), minLevel, maxLevel);
}
......@@ -176,7 +176,7 @@ public:
return numRefinements;
}
void calcMeshSizes(double& minH, double& maxH, int& minLevel, int& maxLevel)
void calcMeshSizes(double& minH, double& maxH, int& minLevel, int& maxLevel, bool allReduce = false)
{
FixVec<WorldVector<double>, VERTEX> coords(mesh->getDim(), NO_INIT);
......@@ -204,10 +204,19 @@ public:
}
minLevel += mesh->getMacroElementLevel();
maxLevel += mesh->getMacroElementLevel();
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
if (allReduce) {
Parallel::mpi::globalMin(minH);
Parallel::mpi::globalMax(maxH);
Parallel::mpi::globalMin(minLevel);
Parallel::mpi::globalMax(maxLevel);
}
#endif
}
double calcMeshSize(ElInfo *elInfo)
double calcMeshSize(ElInfo *elInfo, bool allReduce = false)
{
FixVec<WorldVector<double>, VERTEX> coords(mesh->getDim(), NO_INIT);
coords = elInfo->getCoords();
......@@ -218,6 +227,10 @@ public:
h = std::max(h,norm(coords[i]-coords[j]));
}
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
if (allReduce)
Parallel::mpi::globalMax(h);
#endif
return h;
}
......
extensions/SingularDirichletBC.h
View file @ 2c1e3427
......@@ -47,6 +47,35 @@ struct ManualPeriodicBC {
};
namespace details {
inline bool dofOnFace(int dim, GeoIndex boundaryPos, int face, Element* el, ElementDofIterator& elDofIter)
{
// all DOFs on the face
if (elDofIter.getPosIndex() == boundaryPos &&
elDofIter.getCurrentElementPos() == face)
return true;
// Vertex DOFs on the face. Since position-number of the face
// is associated with the vertex opposite to the face all vertices
// with position != face number can be collected
if (elDofIter.getPosIndex() == VERTEX &&
elDofIter.getCurrentElementPos() != face)
return true;
// in 3d the DOFs on the edges on the face must be collected
if (dim == 3) {
if (elDofIter.getPosIndex() == EDGE) {
int localEdge = elDofIter.getCurrentElementPos();
if (el->getEdgeOfFace(face,0) == localEdge ||
el->getEdgeOfFace(face,1) == localEdge ||
el->getEdgeOfFace(face,2) == localEdge)
return true;
}
}
return false;
}
/**
* traverses the mesh and stores all DegreeOfFreedom, that lie on the boundary with index 'nr'
......@@ -62,6 +91,9 @@ namespace details {
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numBasFcts = basFcts->getNumber();
std::vector<DegreeOfFreedom> localIndices(numBasFcts);
ElementDofIterator elDofIter(feSpace, true);
GeoIndex boundaryPos = (dim == 1 ? VERTEX : (dim == 2 ? EDGE : FACE));
indices.clear();
......@@ -72,12 +104,19 @@ namespace details {
while (elInfo) {
for (int face = 0; face < dim + 1; face++) {
if (elInfo->getBoundary(face) == nr) {
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
for (int i = 0; i < numBasFcts; i++) {
elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
if ((*meshIndicator)(coord))
indices.insert(localIndices[i]);
if (!meshIndicator) {
elDofIter.reset(elInfo->getElement());
do {
if (dofOnFace(dim, boundaryPos, face, elInfo->getElement(), elDofIter))
indices.insert(elDofIter.getDof());
} while (elDofIter.next());
} else {
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
for (int i = 0; i < numBasFcts; i++) {
elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
if ((*meshIndicator)(coord))
indices.insert(localIndices[i]);
}
}
break;
}
......@@ -344,6 +383,26 @@ struct DirichletBcData {
boundary_nr(nr), meshIndicator(meshIndicator_),
val0(0.0), val1(NULL), val2(&val_) { pos.set(0.0); }
// pos = meshindicator + boundary_nr
DirichletBcData(int i_, int j_, BoundaryTypeContainer nr, double val_)
: row(i_), col(j_), posType(4), valueType(0),
idx(0), SignedDistFct(NULL), SignedDistDOF(NULL),
boundary_nr(nr.b), meshIndicator(NULL),
val0(val_), val1(NULL), val2(NULL) { pos.set(0.0); }
DirichletBcData(int i_, int j_, BoundaryTypeContainer nr, DOFVector<double> &val_)
: row(i_), col(j_), posType(4), valueType(1),
idx(0), SignedDistFct(NULL), SignedDistDOF(NULL),
boundary_nr(nr.b), meshIndicator(NULL),
val0(0.0), val1(&val_), val2(NULL) { pos.set(0.0); }
DirichletBcData(int i_, int j_, BoundaryTypeContainer nr, AbstractFunction<double, WorldVector<double> > &val_)
: row(i_), col(j_), posType(4), valueType(2),
idx(0), SignedDistFct(NULL), SignedDistDOF(NULL),
boundary_nr(nr.b), meshIndicator(NULL),
val0(0.0), val1(NULL), val2(&val_) { pos.set(0.0); }
// pos = meshindicator
DirichletBcData(int i_, int j_, AbstractFunction<bool, WorldVector<double> > *meshIndicator_, double val_)
: row(i_), col(j_), posType(4), valueType(0),
......@@ -367,41 +426,42 @@ struct DirichletBcData {
{
std::vector<DegreeOfFreedom> indices_;
std::vector<double> values_;
DegreeOfFreedom idx_ = 0;
int idx_ = -1;
switch (posType) {
case 0:
case 0: // WorldVector
if (val1 == NULL && SignedDistDOF == NULL) {
val1 = new DOFVector<double>(feSpace, "vel1");
}
if (val1 != NULL) {
if (!val1->getDofIdxAtPoint(pos, idx_))
idx_ = 0;
idx_ = -1;
} else if (SignedDistDOF != NULL) {
if (!SignedDistDOF->getDofIdxAtPoint(pos, idx_))
idx_ = 0;
idx_ = -1;
}
indices_.push_back(idx_);
if (idx >= 0)
indices_.push_back(idx_);
break;
case 1:
case 1: // DegreeOfFreedom
indices_.push_back(idx);
break;
case 2:
case 2: // AbstractFunction (signedDistFct)
details::getImplicitIndices(feSpace, SignedDistFct, indices_);
break;
case 3:
case 3: // DOFVector (signedDistDOF)
details::getImplicitIndices(feSpace, SignedDistDOF, indices_);
break;
case 4:
std::set<DegreeOfFreedom> indices__;
if (boundary_nr == 0)