Commit 34287643 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Introduction of multi-mesh method for gradients of basis function. Some code...

Introduction of multi-mesh method for gradients of basis function. Some code refactoring, removed NEW and DELETE macros.
parent 42675765
......@@ -3,11 +3,10 @@
const Flag HL_SignedDist::VEL_EXT = 0X01L;
const Flag HL_SignedDist::VEL_EXT_FROM_VEL_FIELD = 0X02L;
void
HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_,
const DOFVector<double> *lS_DOF_,
DOFVector<double> *sD_DOF_,
ElementFunction<double> *elFct)
void HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_,
const DOFVector<double> *lS_DOF_,
DOFVector<double> *sD_DOF_,
ElementFunction<double> *elFct)
{
adaptInfo = adaptInfo_;
......@@ -31,8 +30,8 @@ HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_,
// ===== Create DOF vector to mark boundary vertices. =====
if (bound_DOF)
DELETE bound_DOF;
bound_DOF = NEW DOFVector<double>(feSpace, "bound_DOF");
delete bound_DOF;
bound_DOF = new DOFVector<double>(feSpace, "bound_DOF");
bound_DOF->set(0.0);
// ===== Boundary vertex and boundary value initialization. =====
......@@ -52,30 +51,28 @@ HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_,
printSignedDistFct();
}
void
HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_,
DOFVector<double> *lS_DOF_)
void HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_, DOFVector<double> *lS_DOF_)
{
TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");
sD_DOF = NEW DOFVector<double>(lS_DOF_->getFESpace(), "sD_DOF");
sD_DOF = new DOFVector<double>(lS_DOF_->getFESpace(), "sD_DOF");
calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, NULL);
// lS_DOF_->copy(*sD_DOF);
*lS_DOF_ = *sD_DOF;
DELETE sD_DOF;
delete sD_DOF;
}
void
HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
DOFVector<double> *origVel_DOF_,
DOFVector<double> *vel_DOF_,
const DOFVector<double> *lS_DOF_,
DOFVector<double> *sD_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
void HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
DOFVector<double> *origVel_DOF_,
DOFVector<double> *vel_DOF_,
const DOFVector<double> *lS_DOF_,
DOFVector<double> *sD_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
{
TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
TEST_EXIT(velExt)("velExt not defined !\n");
......@@ -86,27 +83,26 @@ HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
if (calcSDFct || sD_DOF_ != NULL) {
calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF_, elFct);
}
else {
} else {
TEST_EXIT(vel_DOF_)("illegal level set function lS_DOF_ !\n");
sD_DOF = NEW DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
sD_DOF = new DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
DELETE sD_DOF;
delete sD_DOF;
}
velExt->printVelDOF(adaptInfo_);
}
void
HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
DOFVector<double> *vel_DOF_,
DOFVector<double> *lS_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
void HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
DOFVector<double> *vel_DOF_,
DOFVector<double> *lS_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
{
TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
TEST_EXIT(velExt)("velExt not defined !\n");
......@@ -116,42 +112,40 @@ HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");
DOFVector<double> *newVel_DOF_ =
NEW DOFVector<double>(vel_DOF_->getFESpace(), "vel_DOF_");
new DOFVector<double>(vel_DOF_->getFESpace(), "vel_DOF_");
velExt->setVelocity(vel_DOF_, newVel_DOF_);
velExt->printOrigVelDOF(adaptInfo_);
if (calcSDFct && elFct == NULL)
if (calcSDFct && elFct == NULL) {
calcSignedDistFct(adaptInfo_, lS_DOF_);
else {
sD_DOF = NEW DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
} else {
sD_DOF = new DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
if (calcSDFct)
lS_DOF_->copy(*sD_DOF);
DELETE sD_DOF;
delete sD_DOF;
}
vel_DOF_->copy(*newVel_DOF_);
velExt->printVelDOF(adaptInfo_);
DELETE newVel_DOF_;
delete newVel_DOF_;
}
void
HL_SignedDist::calcVelocityExt(
AdaptInfo *adaptInfo_,
std::vector<DOFVector<double> *> origVel_DOF_,
std::vector<DOFVector<double> *> vel_DOF_,
const DOFVector<double> *lS_DOF_,
DOFVector<double> *sD_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
void HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
std::vector<DOFVector<double> *> origVel_DOF_,
std::vector<DOFVector<double> *> vel_DOF_,
const DOFVector<double> *lS_DOF_,
DOFVector<double> *sD_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
{
TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
TEST_EXIT(velExt)("velExt not defined !\n");
......@@ -162,27 +156,26 @@ HL_SignedDist::calcVelocityExt(
if (calcSDFct || sD_DOF_ != NULL) {
calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF_, elFct);
}
else {
} else {
TEST_EXIT(vel_DOF_[0])("illegal level set function lS_DOF_ !\n");
sD_DOF = NEW DOFVector<double>(vel_DOF_[0]->getFESpace(), "sD_DOF");
sD_DOF = new DOFVector<double>(vel_DOF_[0]->getFESpace(), "sD_DOF");
calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
DELETE sD_DOF;
delete sD_DOF;
}
velExt->printVelDOF(adaptInfo_);
}
void
HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
std::vector<DOFVector<double> *> vel_DOF_,
DOFVector<double> *lS_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
void HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
std::vector<DOFVector<double> *> vel_DOF_,
DOFVector<double> *lS_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
{
TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
TEST_EXIT(velExt)("velExt not defined !\n");
......@@ -193,30 +186,27 @@ HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");
std::vector<DOFVector<double> *> newVel_DOF_(nVelDOFs);
for (int nV=0; nV<nVelDOFs; ++nV) {
TEST_EXIT(vel_DOF_[nV])("velocity vector vel_DOF_ not defined !\n");
newVel_DOF_[nV] =
NEW DOFVector<double>((vel_DOF_[nV])->getFESpace(), "vel_DOF_");
for (int nV = 0; nV < nVelDOFs; ++nV) {
TEST_EXIT(vel_DOF_[nV])("velocity vector vel_DOF_ not defined !\n");
newVel_DOF_[nV] = new DOFVector<double>((vel_DOF_[nV])->getFESpace(), "vel_DOF_");
}
velExt->setVelocity(vel_DOF_, newVel_DOF_);
velExt->printOrigVelDOF(adaptInfo_);
if (calcSDFct && elFct == NULL)
if (calcSDFct && elFct == NULL) {
calcSignedDistFct(adaptInfo_, lS_DOF_);
else {
sD_DOF = NEW DOFVector<double>(vel_DOF_[0]->getFESpace(), "sD_DOF");
} else {
sD_DOF = new DOFVector<double>(vel_DOF_[0]->getFESpace(), "sD_DOF");
calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
if (calcSDFct)
lS_DOF_->copy(*sD_DOF);
DELETE sD_DOF;
delete sD_DOF;
}
for (int nV=0; nV<nVelDOFs; ++nV)
......@@ -224,19 +214,18 @@ HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
velExt->printVelDOF(adaptInfo_);
for (int nV=0; nV<nVelDOFs; ++nV)
DELETE newVel_DOF_[nV];
for (int nV = 0; nV < nVelDOFs; ++nV)
delete newVel_DOF_[nV];
}
void
HL_SignedDist::calcVelocityExtFromVelocityField(
AdaptInfo *adaptInfo_,
std::vector<DOFVector<double> *> &velField_,
DOFVector<double> *vel_DOF_,
const DOFVector<double> *lS_DOF_,
DOFVector<double> *sD_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
void HL_SignedDist::calcVelocityExtFromVelocityField(AdaptInfo *adaptInfo_,
std::vector<DOFVector<double> *> &velField_,
DOFVector<double> *vel_DOF_,
const DOFVector<double> *lS_DOF_,
DOFVector<double> *sD_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
{
TEST_EXIT(velExtType.isSet(VEL_EXT_FROM_VEL_FIELD))
("illegal velocity extension type !\n");
......@@ -247,29 +236,27 @@ HL_SignedDist::calcVelocityExtFromVelocityField(
if (calcSDFct) {
calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF_, elFct);
}
else {
} else {
TEST_EXIT(vel_DOF_)("illegal velocity vector vel_DOF_ !\n");
sD_DOF = NEW DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
sD_DOF = new DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
DELETE sD_DOF;
delete sD_DOF;
}
velExt->printVelDOF(adaptInfo_);
}
void
HL_SignedDist::calcVelocityExtFromVelocityField(
AdaptInfo *adaptInfo_,
std::vector<DOFVector<double> *> &velField_,
DOFVector<double> *vel_DOF_,
DOFVector<double> *lS_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
void HL_SignedDist::calcVelocityExtFromVelocityField(AdaptInfo *adaptInfo_,
std::vector<DOFVector<double> *> &velField_,
DOFVector<double> *vel_DOF_,
DOFVector<double> *lS_DOF_,
bool calcSDFct,
ElementFunction<double> *elFct)
{
TEST_EXIT(velExtType.isSet(VEL_EXT_FROM_VEL_FIELD))
("illegal velocity extension type !\n");
......@@ -283,44 +270,42 @@ HL_SignedDist::calcVelocityExtFromVelocityField(
if (calcSDFct && elFct == NULL) {
calcSignedDistFct(adaptInfo_, lS_DOF_);
}
else {
} else {
TEST_EXIT(vel_DOF_)("illegal velocity vector vel_DOF_ !\n");
sD_DOF = NEW DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
sD_DOF = new DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
if (calcSDFct)
lS_DOF_->copy(*sD_DOF);
DELETE sD_DOF;
delete sD_DOF;
}
velExt->printVelDOF(adaptInfo_);
}
void
HL_SignedDist::initialize(ElementFunction<double> *elFct)
void HL_SignedDist::initialize(ElementFunction<double> *elFct)
{
if (elUpdate)
DELETE elUpdate;
delete elUpdate;
if (bndElDist)
DELETE bndElDist;
delete bndElDist;
if (elLS)
DELETE elLS;
delete elLS;
if (phi)
DELETE phi;
delete phi;
// ===== Create ElementLevelSet. =====
if (elFct == NULL) {
phi = NEW ElementFunctionDOFVec<double>(lS_DOF);
elLS = NEW ElementLevelSet("ElementLevelSet", phi, feSpace->getMesh());
}
else {
elLS = NEW ElementLevelSet("ElementLevelSet", elFct, feSpace->getMesh());
phi = new ElementFunctionDOFVec<double>(lS_DOF);
elLS = new ElementLevelSet("ElementLevelSet", phi, feSpace->getMesh());
} else {
elLS = new ElementLevelSet("ElementLevelSet", elFct, feSpace->getMesh());
}
......@@ -329,13 +314,13 @@ HL_SignedDist::initialize(ElementFunction<double> *elFct)
GET_PARAMETER(0, name + "->boundary initialization", "%d",
&boundInitFlag);
switch(boundInitFlag) {
case 0: bndElDist = NEW BoundaryElementLevelSetDist(elLS, dim);
case 0: bndElDist = new BoundaryElementLevelSetDist(elLS, dim);
break;
case 1: bndElDist = NEW BoundaryElementTopDist(elLS, dim, velExt);
case 1: bndElDist = new BoundaryElementTopDist(elLS, dim, velExt);
break;
case 2: bndElDist = NEW BoundaryElementEdgeDist(elLS, dim);
case 2: bndElDist = new BoundaryElementEdgeDist(elLS, dim);
break;
case 3: bndElDist = NEW BoundaryElementNormalDist(elLS, dim);
case 3: bndElDist = new BoundaryElementNormalDist(elLS, dim);
break;
default: ERROR_EXIT("illegal boundary initialization !\n");
break;
......@@ -346,30 +331,23 @@ HL_SignedDist::initialize(ElementFunction<double> *elFct)
("velocity extension only works with topological boundary element initialization method !\n");
switch(dim) {
case 2: elUpdate = NEW ElementUpdate_2d(velExt);
case 2: elUpdate = new ElementUpdate_2d(velExt);
break;
case 3: elUpdate = NEW ElementUpdate_3d(velExt);
case 3: elUpdate = new ElementUpdate_3d(velExt);
break;
default: ERROR_EXIT("illegal dimension !\n");
break;
}
}
void
HL_SignedDist::setSign()
void HL_SignedDist::setSign()
{
DOFVector<double>::Iterator it_sD(const_cast<DOFVector<double> *>(sD_DOF),
USED_DOFS);
DOFVector<double>::Iterator it_lS(const_cast<DOFVector<double> *>(lS_DOF),
USED_DOFS);
for(it_sD.reset(), it_lS.reset();
!it_sD.end();
++it_sD, ++it_lS) {
// Vertex lies inside the zero level set ?
if ((*it_lS) < 0) {
(*it_sD) *= -1;
}
}
DOFVector<double>::Iterator it_sD(const_cast<DOFVector<double> *>(sD_DOF), USED_DOFS);
DOFVector<double>::Iterator it_lS(const_cast<DOFVector<double> *>(lS_DOF), USED_DOFS);
// Vertex lies inside the zero level set ?
for (it_sD.reset(), it_lS.reset(); !it_sD.end(); ++it_sD, ++it_lS)
if ((*it_lS) < 0)
(*it_sD) *= -1;
}
......@@ -31,13 +31,13 @@ public:
TEST_EXIT(i < velDOF.size())("illegal index !\n");
FileWriter *fileWriter = NEW FileWriter(
FileWriter *fileWriter = new FileWriter(
"VelocityExt->velocity output",
(velDOF[i])->getFESpace()->getMesh(),
const_cast<DOFVector<double> *>(velDOF[i]));
fileWriter->writeFiles(adaptInfo, false);
DELETE fileWriter;
delete fileWriter;
}
/// Print origVelDOF to file.
......@@ -47,13 +47,13 @@ public:
TEST_EXIT(i < origVelDOF.size())("illegal index !\n");
FileWriter *fileWriter = NEW FileWriter(
FileWriter *fileWriter = new FileWriter(
"VelocityExt->interface velocity output",
(origVelDOF[i])->getFESpace()->getMesh(),
const_cast<DOFVector<double> *>(origVelDOF[i]));
fileWriter->writeFiles(adaptInfo, false);
DELETE fileWriter;
delete fileWriter;
};
/**
......
......@@ -113,15 +113,25 @@ namespace AMDiS {
secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
ElementMatrix &m =
smallElInfo->getSubElemCoordsMat_so(rowFESpace->getBasisFcts()->getDegree());
smallElInfo->getSubElemGradCoordsMat(rowFESpace->getBasisFcts()->getDegree());
tmpMat = m * mat;
mat = tmpMat;
}
if (firstOrderAssemblerGrdPsi) {
ERROR_EXIT("Not yet implemented for first order assembler!\n");
firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
if (largeElInfo == colElInfo) {
ERROR_EXIT("123 Not yet implemented for first order assembler!\n");
} else {
firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
ElementMatrix &m =
smallElInfo->getSubElemCoordsMat(rowFESpace->getBasisFcts()->getDegree());
tmpMat = mat * trans(m);
mat = tmpMat;
}
}
if (firstOrderAssemblerGrdPhi) {
......
......@@ -16,6 +16,8 @@ namespace AMDiS {
std::vector<std::map<unsigned long, mtl::dense2D<double> > > ElInfo::subElemMatrices(4);
std::vector<std::map<unsigned long, mtl::dense2D<double> > > ElInfo::subElemGradMatrices(4);
ElInfo::ElInfo(Mesh *aMesh)
: mesh(aMesh),
element(NULL),
......
......@@ -244,9 +244,9 @@ namespace AMDiS {
return subElemMatrices[degree][refinementPath];
}
virtual mtl::dense2D<double>& getSubElemCoordsMat_so(int degree) const
virtual mtl::dense2D<double>& getSubElemGradCoordsMat(int degree) const
{
return subElemMatrices[degree][refinementPath];
return subElemGradMatrices[degree][refinementPath];
}
/** \} */
......@@ -544,10 +544,12 @@ namespace AMDiS {
*/
mtl::dense2D<double> subElemCoordsMat;
mtl::dense2D<double> subElemCoordsMat_so;
mtl::dense2D<double> subElemGradCoordsMat;
public:
static std::vector<std::map<unsigned long, mtl::dense2D<double> > > subElemMatrices;
static std::vector<std::map<unsigned long, mtl::dense2D<double> > > subElemGradMatrices;
/** \brief
* child_vertex[el_type][child][i] = father's local vertex index of new
......
......@@ -40,11 +40,7 @@ namespace AMDiS {
{1.0, 0.75, 0.0},
{0.0, 0.375, 1.0}};
mtl::dense2D<double> ElInfo1d::mat_d2_right(mat_d2_right_val);
double ElInfo1d::test1_val[2][2] = {{1.0, 0.0},
{0.0, 1.0}};
mtl::dense2D<double> ElInfo1d::test2_val(test1_val);
void ElInfo1d::fillMacroInfo(const MacroElement * mel)
{
......@@ -349,27 +345,24 @@ namespace AMDiS {
}
mtl::dense2D<double>& ElInfo1d::getSubElemCoordsMat_so(int degree) const
mtl::dense2D<double>& ElInfo1d::getSubElemGradCoordsMat(int degree) const
{
FUNCNAME("ElInfo1d::getSubElemCoordsMat_so()");
FUNCNAME("ElInfo1d::getSubElemGradCoordsMat()");
ERROR_EXIT("Reimplement!\n");
TEST_EXIT(degree == 1)("Not supported for basis functions with degree > 1!\n");
return subElemMatrices[degree][refinementPath];
using namespace mtl;
#if 0
switch (degree) {
case 1:
mat = test2_val;
if (subElemMatrices[degree].count(refinementPath) == 0) {
dense2D<double> mat(mat_d1);
for (int i = 0; i < refinementPathLength; i++)
mat *= 0.5;
break;
default:
ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
subElemGradMatrices[1][refinementPath] = mat;
}
#endif
return subElemGradMatrices[degree][refinementPath];
}
......
......@@ -64,7 +64,7 @@ namespace AMDiS {
mtl::dense2D<double>& getSubElemCoordsMat(int degree) const;
mtl::dense2D<double>& getSubElemCoordsMat_so(int degree) const;
mtl::dense2D<double>& getSubElemGradCoordsMat(int degree) const;
protected:
static double mat_d1_val[2][2];
......@@ -90,11 +90,6 @@ namespace AMDiS {
static double mat_d2_right_val[3][3];
static mtl::dense2D<double> mat_d2_right;
static double test1_val[2][2];
static mtl::dense2D<double> test2_val;
};
}
......
......@@ -723,12 +723,24 @@ namespace AMDiS {
}
mtl::dense2D<double>& ElInfo2d::getSubElemCoordsMat_so(int degree) const
mtl::dense2D<double>& ElInfo2d::getSubElemGradCoordsMat(int degree) const
{
FUNCNAME("ElInfo2d::getSubElemCoordsMat_so()");
FUNCNAME(