diff --git a/AMDiS/src/DOFAdmin.h b/AMDiS/src/DOFAdmin.h
index 7338108726052a13070882496b20e758910393c6..1b99e91da1df3eb3336c1237d3f9ab0e566b2bc5 100644
--- a/AMDiS/src/DOFAdmin.h
+++ b/AMDiS/src/DOFAdmin.h
@@ -225,7 +225,13 @@ namespace AMDiS {
*/
/// Sets \ref nrDOF[i] = v
- void setNumberOfDOFs(int i, int v);
+ void setNumberOfDOFs(int i, int v);
+
+ /// Sets all values of \ref nrDOF
+ void setNumberOfDOFs(DimVec<int> v)
+ {
+ nrDOF = v;
+ }
/// Sets \ref nr0DOF[i] = v
void setNumberOfPreDOFs(int i, int v);
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index bdc099ae547abc46b8b0732ffee411bc89913c5f..edc2e1252fa2e7b914bc77e9d323fd5a8d74b4b4 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -749,9 +749,10 @@ namespace AMDiS {
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
- if (quad && quadFast)
+ if (quad && quadFast) {
TEST_EXIT_DBG(quad == quadFast->getQuadrature())
("quad != quadFast->quadrature\n");
+ }
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
diff --git a/AMDiS/src/DataCollector.cc b/AMDiS/src/DataCollector.cc
index 72e965f6cacfa37055a8003433a10c7fa38f3ad5..c4643feb8f4bd7708d5ecc6bda347815b735ab36 100644
--- a/AMDiS/src/DataCollector.cc
+++ b/AMDiS/src/DataCollector.cc
@@ -16,12 +16,18 @@ namespace AMDiS {
level(l),
traverseFlag(flag),
feSpace(fe),
+ nVertices(0),
+ nElements(0),
+ nInterpPoints(0),
+ nConnection(0),
+ elementDataCollected(false),
+ valueDataCollected(false),
+ periodicDataCollected(false),
writeElem(writeElemFct)
{
FUNCNAME("DataCollector::DataCollector()");
- TEST_EXIT(feSpace)("no feSpace\n");
- TEST_EXIT(values)("no value Vector\n");
+ TEST_EXIT(feSpace)("No finite elem space defined!\n");
// === get mesh ===
mesh = feSpace->getMesh();
@@ -29,36 +35,26 @@ namespace AMDiS {
// === get admin ===
localAdmin = const_cast<DOFAdmin*>(feSpace->getAdmin());
// === create vertex info vector ===
- vertexInfos = new DOFVector< std::list<VertexInfo> >(feSpace, "vertex infos");
+ vertexInfos = new DOFVector<std::list<VertexInfo> >(feSpace, "vertex infos");
interpPointInd = new DOFVector<int>(feSpace, "interpolation point indices");
interpPointCoords = new DOFVector< std::list<WorldVector<double> > >(feSpace, "interpolation point coordinates");
dofCoord = new DOFVector< std::list<WorldVector<double> > >(feSpace, "dof coords");
- dim = mesh->getDim();
-
+ dim = mesh->getDim();
nPreDofs = localAdmin->getNumberOfPreDOFs(VERTEX);
- nVertices = 0;
- nElements = 0;
- nInterpPoints = 0;
- nConnection = 0;
-
- elementDataCollected = false;
- valueDataCollected = false;
- periodicDataCollected = false;
-
- vertexCoords = new WorldVector<double>;
}
+
DataCollector::~DataCollector()
{
delete vertexInfos;
delete interpPointInd;
delete interpPointCoords;
delete dofCoord;
- delete vertexCoords;
}
+
void DataCollector::fillAllData()
{
if (!elementDataCollected)
@@ -71,6 +67,7 @@ namespace AMDiS {
startCollectingValueData();
}
+
void DataCollector::startCollectingElementData()
{
FUNCNAME("DataCollector::startCollectingElementData()");
@@ -107,6 +104,7 @@ namespace AMDiS {
elementDataCollected = true;
}
+
void DataCollector::startCollectingValueData()
{
FUNCNAME("DataCollector::startCollectingValueData()");
@@ -156,6 +154,7 @@ namespace AMDiS {
valueDataCollected = true;
}
+
void DataCollector::startCollectingPeriodicData()
{
FUNCNAME("DataCollector::startCollectingPeriodicData()");
@@ -200,6 +199,7 @@ namespace AMDiS {
periodicDataCollected = true;
}
+
void DataCollector::addElementData(ElInfo* elInfo)
{
FUNCNAME("DataCollector::addElementData()");
@@ -272,10 +272,12 @@ namespace AMDiS {
elements.push_back(elementInfo);
}
+
void DataCollector::addValueData(ElInfo *elInfo)
{
FUNCNAME("DataCollector::addValueData()");
+ WorldVector<double> vertexCoords;
basisFcts->getLocalIndices(elInfo->getElement(), localAdmin, localDOFs);
// First, traverse all DOFs at the vertices of the element, determine
@@ -286,18 +288,18 @@ namespace AMDiS {
(*interpPointInd)[dofi] = -2; // mark as vertex
// get coords of this vertex
- *vertexCoords = elInfo->getCoord(i);
+ vertexCoords = elInfo->getCoord(i);
// search for coords at this dof
std::list<WorldVector<double> >::iterator it =
find((*dofCoord)[dofi].begin(),
(*dofCoord)[dofi].end(),
- *vertexCoords);
+ vertexCoords);
// coords not yet in list?
if (it == (*dofCoord)[dofi].end()) {
// add new coords to list
- (*dofCoord)[dofi].push_back(*vertexCoords);
+ (*dofCoord)[dofi].push_back(vertexCoords);
}
}
@@ -307,7 +309,7 @@ namespace AMDiS {
for (int i = mesh->getGeo(VERTEX); i < nBasisFcts; i++) {
DegreeOfFreedom dofi = localDOFs[i];
- elInfo->coordToWorld(*basisFcts->getCoords(i), *vertexCoords);
+ elInfo->coordToWorld(*basisFcts->getCoords(i), vertexCoords);
if ((*interpPointInd)[dofi] == -1) {
// mark as interpolation point
@@ -318,16 +320,17 @@ namespace AMDiS {
std::list<WorldVector<double> >::iterator it =
find((*interpPointCoords)[dofi].begin(),
(*interpPointCoords)[dofi].end(),
- *vertexCoords);
+ vertexCoords);
if (it == (*interpPointCoords)[dofi].end()) {
- (*interpPointCoords)[dofi].push_back(*vertexCoords);
+ (*interpPointCoords)[dofi].push_back(vertexCoords);
nInterpPoints++;
}
}
}
}
+
void DataCollector::addInterpData(ElInfo *elInfo)
{
FUNCNAME("DataCollector::addInterpData()");
@@ -341,6 +344,7 @@ namespace AMDiS {
interpPoints.push_back(elemInterpPoints);
}
+
void DataCollector::addPeriodicData(ElInfo *elInfo)
{
FUNCNAME("DataCollector::addPeriodicData");
@@ -401,6 +405,7 @@ namespace AMDiS {
}
}
+
std::list<ElementInfo>* DataCollector::getElementInfos()
{
if (!elementDataCollected)
@@ -409,6 +414,7 @@ namespace AMDiS {
return &elements;
}
+
DOFVector< std::list<VertexInfo> >* DataCollector::getVertexInfos()
{
if (!elementDataCollected)
@@ -417,6 +423,7 @@ namespace AMDiS {
return vertexInfos;
}
+
std::list<PeriodicInfo>* DataCollector::getPeriodicInfos()
{
if (!periodicDataCollected)
@@ -425,6 +432,7 @@ namespace AMDiS {
return &periodicInfos;
}
+
int DataCollector::getNumberVertices()
{
if (!elementDataCollected)
@@ -433,6 +441,7 @@ namespace AMDiS {
return nVertices;
}
+
int DataCollector::getNumberElements()
{
if (!elementDataCollected)
@@ -441,6 +450,7 @@ namespace AMDiS {
return nElements;
}
+
int DataCollector::getNumberInterpPoints()
{
if (!valueDataCollected)
@@ -448,6 +458,7 @@ namespace AMDiS {
return nInterpPoints;
}
+
int DataCollector::getNumberConnections()
{
@@ -457,15 +468,6 @@ namespace AMDiS {
return nConnection;
}
- const FiniteElemSpace* DataCollector::getFeSpace()
- {
- return feSpace;
- }
-
- Mesh* DataCollector::getMesh()
- {
- return mesh;
- }
DOFVector<double>* DataCollector::getValues()
{
@@ -475,6 +477,7 @@ namespace AMDiS {
return values;
}
+
DOFVector< std::list<WorldVector<double> > >* DataCollector::getDofCoords()
{
if (!valueDataCollected)
@@ -483,6 +486,7 @@ namespace AMDiS {
return dofCoord;
}
+
DOFVector<int>* DataCollector::getInterpPointInd()
{
if (!valueDataCollected)
@@ -491,6 +495,7 @@ namespace AMDiS {
return interpPointInd;
}
+
DOFVector< std::list<WorldVector<double> > >* DataCollector::getInterpPointCoords()
{
if (!valueDataCollected)
@@ -499,6 +504,7 @@ namespace AMDiS {
return interpPointCoords;
}
+
std::vector< std::vector<int> >* DataCollector::getInterpPoints()
{
if (!valueDataCollected)
diff --git a/AMDiS/src/DataCollector.h b/AMDiS/src/DataCollector.h
index 19073164a92677efba88e7ba9bab05a35018ac4b..35826204a3d3447378c8baaccc6d893c0d899ab2 100644
--- a/AMDiS/src/DataCollector.h
+++ b/AMDiS/src/DataCollector.h
@@ -43,7 +43,7 @@ namespace AMDiS {
public:
/// Constructor
DataCollector(const FiniteElemSpace *feSpace,
- DOFVector<double> *values,
+ DOFVector<double> *values = NULL,
int level = -1,
Flag traverseFlag = Mesh::CALL_LEAF_EL,
bool (*writeElem)(ElInfo*) = NULL);
@@ -60,7 +60,10 @@ namespace AMDiS {
DOFVector< std::list<VertexInfo> >* getVertexInfos();
/// Returns the finite element space of the problem.
- const FiniteElemSpace* getFeSpace();
+ const FiniteElemSpace* getFeSpace()
+ {
+ return feSpace;
+ }
/// Returns vector with value information.
DOFVector<double>* getValues();
@@ -93,7 +96,10 @@ namespace AMDiS {
int getNumberConnections();
/// Returns the mesh of the problem.
- Mesh* getMesh();
+ Mesh* getMesh()
+ {
+ return mesh;
+ }
void setMesh(Mesh *m)
{
@@ -165,16 +171,16 @@ namespace AMDiS {
DOFVector<int> *interpPointInd;
/// Stores for each simplex the interpolation points.
- std::vector< std::vector<int> > interpPoints;
+ std::vector<std::vector<int> > interpPoints;
/** \brief
* Stores for each DOF a list of its coordinates. If there are now periodic
* boundaries than there is also only one coordinate per DOF.
*/
- DOFVector< std::list<WorldVector<double> > > *interpPointCoords;
+ DOFVector<std::list<WorldVector<double> > > *interpPointCoords;
/// list of coords for each dof
- DOFVector< std::list<WorldVector<double> > > *dofCoord;
+ DOFVector<std::list<WorldVector<double> > > *dofCoord;
/// List that stores an ElementInfo for each element.
std::list<ElementInfo> elements;
@@ -183,7 +189,7 @@ namespace AMDiS {
std::list<PeriodicInfo> periodicInfos;
/// Stores a list of vertex infos for each dof.
- DOFVector< std::list<VertexInfo> > *vertexInfos;
+ DOFVector<std::list<VertexInfo> > *vertexInfos;
/** \brief
* periodicConnections[i][j] stores whether the connection at side j of
diff --git a/AMDiS/src/FileWriter.cc b/AMDiS/src/FileWriter.cc
index 00cb4fbb9568bdcf8217c232a937a05364b3732c..13f5425461e27e423f4a785cccbe960895797a1f 100644
--- a/AMDiS/src/FileWriter.cc
+++ b/AMDiS/src/FileWriter.cc
@@ -158,6 +158,7 @@ namespace AMDiS {
}
}
+
void FileWriter::writeFiles(AdaptInfo *adaptInfo,
bool force,
int level,
@@ -170,7 +171,7 @@ namespace AMDiS {
return;
// Containers, which store the data to be written;
- std::vector< DataCollector* > dataCollectors(solutionVecs.size());
+ std::vector<DataCollector* > dataCollectors(solutionVecs.size());
if (writeElem) {
for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++)
@@ -219,10 +220,10 @@ namespace AMDiS {
if (writeAMDiSFormat) {
MacroWriter::writeMacro(dataCollectors[0],
- const_cast<char*>((fn + amdisMeshExt).c_str()),
- adaptInfo ? adaptInfo->getTime() : 0.0);
+ const_cast<char*>((fn + amdisMeshExt).c_str()),
+ adaptInfo ? adaptInfo->getTime() : 0.0);
MSG("macro file written to %s\n", (fn + amdisMeshExt).c_str());
-
+
ValueWriter::writeValues(dataCollectors[0],
(fn + amdisDataExt).c_str(),
adaptInfo ? adaptInfo->getTime() : 0.0);
diff --git a/AMDiS/src/FiniteElemSpace.cc b/AMDiS/src/FiniteElemSpace.cc
index 50e202cabc0435607b4e16b9f3a2b562b75eefa0..cc75ccf07b9926e33baf4a334302aa802dbf2a23 100644
--- a/AMDiS/src/FiniteElemSpace.cc
+++ b/AMDiS/src/FiniteElemSpace.cc
@@ -20,7 +20,8 @@ namespace AMDiS {
{
FUNCNAME("FiniteElemSpace::FiniteElemSpace()");
- TEST_EXIT(mesh)("no Mesh\n");
+ TEST_EXIT(mesh)("No mesh!\n");
+ TEST_EXIT(basFcts)("No basis functions!\n");
if (!admin) {
const DOFAdmin *admin_local = NULL;
diff --git a/AMDiS/src/FirstOrderTerm.cc b/AMDiS/src/FirstOrderTerm.cc
index c39a63f2fe639a0b6f21d0dad59ceee4068e0ef5..8380049e17c0654f3485de106eb6fe44784fdbee 100644
--- a/AMDiS/src/FirstOrderTerm.cc
+++ b/AMDiS/src/FirstOrderTerm.cc
@@ -62,7 +62,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int dow = Global::getGeo(WORLD);
@@ -101,7 +101,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double f)
{
int dow = Global::getGeo(WORLD);
@@ -112,6 +112,7 @@ namespace AMDiS {
double resultQP = 0.0;
for (int i = 0; i < dow; i++)
resultQP += grdUhAtQP[iq][i];
+
result[iq] += f * factor * resultQP;
}
}
@@ -139,7 +140,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double f)
{
int dow = Global::getGeo(WORLD);
@@ -192,7 +193,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
if (grdUhAtQP) {
@@ -239,7 +240,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
if (grdUhAtQP) {
@@ -273,7 +274,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double f)
{
int dow = Global::getGeo(WORLD);
@@ -331,7 +332,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
if (grdUhAtQP) {
@@ -365,7 +366,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
if (grdUhAtQP) {
@@ -415,7 +416,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
if (grdUhAtQP)
@@ -461,7 +462,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
if (grdUhAtQP)
@@ -527,7 +528,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
if (grdUhAtQP)
@@ -602,7 +603,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
if (grdUhAtQP)
@@ -680,7 +681,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -782,7 +783,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
diff --git a/AMDiS/src/FirstOrderTerm.h b/AMDiS/src/FirstOrderTerm.h
index 457fefba7dfc4358283ab24205f8ae92a74f72f5..fc2b126c3f85b3e9e95a077075ad779c72f89580 100644
--- a/AMDiS/src/FirstOrderTerm.h
+++ b/AMDiS/src/FirstOrderTerm.h
@@ -56,7 +56,7 @@ namespace AMDiS {
const ElementVector&,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -230,7 +230,7 @@ namespace AMDiS {
const ElementVector&,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *,
- double *result,
+ ElementVector& result,
double factor)
{
if (grdUhAtQP)
@@ -281,7 +281,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -324,7 +324,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -363,7 +363,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -403,7 +403,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector &result,
double factor);
protected:
@@ -442,7 +442,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -483,7 +483,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -527,7 +527,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -573,7 +573,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -613,7 +613,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -646,7 +646,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -679,7 +679,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -714,7 +714,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -753,7 +753,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
@@ -798,7 +798,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
protected:
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index 1d97fa00c9c6e8fabe964d6163e4ef4aec555fd3..1ea009075998c14db74617d0f6cf9932eedf1636 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -21,15 +21,15 @@
namespace AMDiS {
- MacroInfo* MacroReader::readMacro(const char *filename,
+ MacroInfo* MacroReader::readMacro(std::string filename,
Mesh* mesh,
- const char *periodicFile,
+ std::string periodicFile,
int check)
{
FUNCNAME("Mesh::readMacro()");
- TEST_EXIT(filename)("no file specified; filename NULL pointer\n");
-
+ TEST_EXIT(filename != "")("no file specified; filename NULL pointer\n");
+
MacroInfo *macroInfo = new MacroInfo();
macroInfo->readAMDiSMacro(filename, mesh);
@@ -39,11 +39,11 @@ namespace AMDiS {
DegreeOfFreedom **dof = macroInfo->dof;
// === read periodic data =================================
- if (periodicFile && (strcmp(periodicFile, "") != 0)) {
+ if (periodicFile != "") {
WARNING("Periodic boundaries may lead to errors in small meshes if element neighbours are not set!\n");
- FILE *file = fopen(periodicFile, "r");
- TEST_EXIT(file)("can't open file %s\n", periodicFile);
+ FILE *file = fopen(periodicFile.c_str(), "r");
+ TEST_EXIT(file)("can't open file %s\n", periodicFile.c_str());
int n;
int dim = mesh->getDim();
@@ -219,7 +219,7 @@ namespace AMDiS {
}
if (!macroInfo->neigh_set) {
- TEST_EXIT(!periodicFile)
+ TEST_EXIT(periodicFile == "")
("periodic boundary condition => element neighbours must be set\n");
computeNeighbours(mesh);
} else {
@@ -310,7 +310,7 @@ namespace AMDiS {
if (mesh->getDim() > 1) {
char filenew[128];
- strncpy(filenew, filename, 128);
+ strncpy(filenew, filename.c_str(), 128);
filenew[127] = 0;
strncat(filenew, ".new", 128);
filenew[127] = 0;
@@ -454,7 +454,7 @@ namespace AMDiS {
/****************************************************************************/
- void MacroInfo::readAMDiSMacro(const char *filename, Mesh* mesh)
+ void MacroInfo::readAMDiSMacro(std::string filename, Mesh* mesh)
{
FUNCNAME("MacroInfo::readAMDiSMacro()");
@@ -467,12 +467,12 @@ namespace AMDiS {
const char *key;
DimVec<int> *ind = NULL;
- TEST_EXIT(filename)("no file specified; filename NULL pointer\n");
- TEST_EXIT(strlen(filename) < static_cast<unsigned int>(127))
+ TEST_EXIT(filename != "")("No filename specified!\n");
+ TEST_EXIT(filename.length() < static_cast<unsigned int>(127))
("can only handle filenames up to 127 characters\n");
- FILE *file = fopen(filename, "r");
- TEST_EXIT(file)("cannot open file %s\n", filename);
+ FILE *file = fopen(filename.c_str(), "r");
+ TEST_EXIT(file)("cannot open file %s\n", filename.c_str());
/****************************************************************************/
/* looking for all keys in the macro file ... */
@@ -487,7 +487,7 @@ namespace AMDiS {
int i_key = get_key_no(key);
TEST_EXIT(i_key >= 0)
("macro file %s must not contain key %s on line %d\n",
- filename, key, line_no);
+ filename.c_str(), key, line_no);
TEST_EXIT(!key_def[i_key])("key %s defined second time on line %d in file %s\n");
sort_key[n_keys++] = i_key;
@@ -504,7 +504,7 @@ namespace AMDiS {
for (j = 0; j < n_keys; j++)
if (sort_key[j] == i_key) break;
TEST_EXIT(j < n_keys)("You do not have specified data for %s in %s\n",
- keys[i_key], filename);
+ keys[i_key], filename.c_str());
for (j = 0; j < n_keys; j++)
if (sort_key[j] == 2) break;
@@ -522,18 +522,18 @@ namespace AMDiS {
case 4:
TEST_EXIT(nv_key < i_key)
("Before reading data for %s, you have to specify the %s in file\n",
- keys[4], keys[2], filename);
+ keys[4], keys[2], filename.c_str());
break;
case 5:
TEST_EXIT(nv_key < i_key && ne_key < i_key)
("Before reading data for %s, you have to specify the %s and %s in file %s\n",
- keys[5], keys[3], keys[2], filename);
+ keys[5], keys[3], keys[2], filename.c_str());
case 6:
case 7:
case 8:
TEST_EXIT(ne_key < i_key)
("Before reading data for %s, you have to specify the %s in file %s\n",
- keys[i_key], keys[3], filename);
+ keys[i_key], keys[3], filename.c_str());
}
}
@@ -544,8 +544,8 @@ namespace AMDiS {
/* and now, reading data ... */
/****************************************************************************/
- file = fopen(filename, "r");
- TEST_EXIT(file)("cannot open file %s\n", filename);
+ file = fopen(filename.c_str(), "r");
+ TEST_EXIT(file)("cannot open file %s\n", filename.c_str());
int result;
@@ -556,7 +556,7 @@ namespace AMDiS {
case 0:
// line "DIM"
result = fscanf(file, "%*s %d", &dim);
- TEST_EXIT(result == 1)("cannot read DIM correctly in file %s\n", filename);
+ TEST_EXIT(result == 1)("cannot read DIM correctly in file %s\n", filename.c_str());
ind = new DimVec<int>(dim, NO_INIT);
@@ -567,7 +567,7 @@ namespace AMDiS {
// line "DIM_OF_WORLD"
result = fscanf(file, "%*s %d", &dow);
TEST_EXIT(result == 1)
- ("cannot read Global::getGeo(WORLD) correctly in file %s\n", filename);
+ ("cannot read Global::getGeo(WORLD) correctly in file %s\n", filename.c_str());
TEST_EXIT(dow == Global::getGeo(WORLD))
("dimension of world = %d != Global::getGeo(WORLD) = %d\n",
dow, Global::getGeo(WORLD));
@@ -579,7 +579,7 @@ namespace AMDiS {
// line "number of vertices"
result = fscanf(file, "%*s %*s %*s %d", &nv);
TEST_EXIT(result == 1)
- ("cannot read number of vertices correctly in file %s\n", filename);
+ ("cannot read number of vertices correctly in file %s\n", filename.c_str());
TEST_EXIT(nv > 0)
("number of vertices = %d must be bigger than 0\n", nv);
@@ -592,7 +592,7 @@ namespace AMDiS {
// line "number of elements"
result = fscanf(file, "%*s %*s %*s %d", &ne);
TEST_EXIT(result == 1)
- ("cannot read number of elements correctly in file %s\n", filename);
+ ("cannot read number of elements correctly in file %s\n", filename.c_str());
TEST_EXIT(ne > 0)
("number of elements = %d must be bigger than 0\n", ne);
@@ -608,7 +608,7 @@ namespace AMDiS {
for (j = 0; j <Global::getGeo(WORLD) ; j++) {
result = fscanf(file, "%lf", &dbl);
TEST_EXIT(result == 1)
- ("error while reading coordinates, check file %s\n", filename);
+ ("error while reading coordinates, check file %s\n", filename.c_str());
coords[i][j] = dbl;
}
}
@@ -626,7 +626,7 @@ namespace AMDiS {
for (int i = 0; i < ne; i++) {
result = read_indices(file, *ind);
TEST_EXIT(result)
- ("cannot read vertex indices of element %d in file %s\n", i, filename);
+ ("cannot read vertex indices of element %d in file %s\n", i, filename.c_str());
for (k = 0; k < mesh->getGeo(VERTEX); k++)
mel_vertex[i][k] = (*ind)[k];
@@ -647,7 +647,7 @@ namespace AMDiS {
result = read_indices(file, *ind);
TEST_EXIT(result)
- ("cannot read boundary type of element %d in file %s\n", i, filename);
+ ("cannot read boundary type of element %d in file %s\n", i, filename.c_str());
// fill boundary of macro-element
MacroReader::fillMelBoundary(mesh,
@@ -700,7 +700,7 @@ namespace AMDiS {
for (int i = 0; i < ne; i++) {
result = fscanf(file, "%d", &j);
TEST_EXIT(result == 1)
- ("cannot read elType of element %d in file %s\n", i, filename);
+ ("cannot read elType of element %d in file %s\n", i, filename.c_str());
if (dim == 3)
(mel)[i]->elType = j;
}
@@ -722,7 +722,7 @@ namespace AMDiS {
for (int i = 0; i < ne; i++) {
result = read_indices(file, *ind);
TEST_EXIT(result)
- ("cannot read boundary projector of element %d in file %s\n", i, filename);
+ ("cannot read boundary projector of element %d in file %s\n", i, filename.c_str());
Projection *projector = Projection::getProjection((*ind)[0]);
@@ -757,7 +757,7 @@ namespace AMDiS {
for (int i = 0; i < ne; i++) {
result = fscanf(file, "%d", &j);
TEST_EXIT(result == 1)
- ("cannot read region of element %d in file %s\n", i, filename);
+ ("cannot read region of element %d in file %s\n", i, filename.c_str());
if (j >= 0) {
Element *el = mel[i]->getElement();
ElementRegion_ED *elementRegion =
@@ -775,7 +775,7 @@ namespace AMDiS {
for (int i = 0; i < ne; i++) {
result = read_indices(file, *ind);
TEST_EXIT(result)
- ("cannot read surface regions of element %d in file %s\n", i, filename);
+ ("cannot read surface regions of element %d in file %s\n", i, filename.c_str());
Element *el = mel[i]->getElement();
@@ -810,25 +810,6 @@ namespace AMDiS {
fclose(file);
}
-
- int macro_type(const char *filename, const char *type)
- {
- const char *fn, *t;
-
- if (strlen(filename) <= strlen(type))
- return(false);
-
- fn = filename;
- while (*fn) fn++;
- t = type;
- while (*t) t++;
-
- while (t != type && *t == *fn) t--;
-
- return(t == type);
- }
-
-
/****************************************************************************/
/* sets the boundary of all edges/faces with no neigbour to a straight */
/* line/face with Dirichlet boundary type */
@@ -836,7 +817,7 @@ namespace AMDiS {
void MacroInfo::dirichletBoundary()
{
- for (int i = 0; i < static_cast<int>( mel.size()); i++) {
+ for (int i = 0; i < static_cast<int>(mel.size()); i++) {
for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
if (mel[i]->neighbour[k])
mel[i]->boundary[k] = INTERIOR;
@@ -849,6 +830,8 @@ namespace AMDiS {
void MacroInfo::fillBoundaryInfo(Mesh *mesh)
{
+ FUNCNAME("MacroInfo::fillBoundaryInfo()");
+
int i, nv = mesh->getNumberOfVertices();
std::deque<MacroElement*>::iterator melIt;
std::vector<BoundaryType> bound(nv);
diff --git a/AMDiS/src/MacroReader.h b/AMDiS/src/MacroReader.h
index 413a28fc9aaa35243ce0daedf836d4f2d06e4ccb..10acb48ad21b36b3f4d188a1bba2a8662bd73cd3 100644
--- a/AMDiS/src/MacroReader.h
+++ b/AMDiS/src/MacroReader.h
@@ -42,9 +42,9 @@ namespace AMDiS {
{
public:
/// Creates a Mesh by reading the macro file with the given filename.
- static MacroInfo* readMacro(const char *filename,
+ static MacroInfo* readMacro(std::string filename,
Mesh* mesh,
- const char *periodicFile,
+ std::string periodicFile,
int check);
protected:
@@ -130,7 +130,7 @@ namespace AMDiS {
* Fills MacroInfo structure.
* Called by Mesh::readMacro(), fills missing information
*/
- void readAMDiSMacro(const char *filename, Mesh* mesh);
+ void readAMDiSMacro(std::string filename, Mesh* mesh);
/// Frees memory of MacroInfo
void clear();
diff --git a/AMDiS/src/MacroWriter.cc b/AMDiS/src/MacroWriter.cc
index 1bd2e2c73b1c419d47e698fa050dcea928e6fdba..ce27abdb36c99803a2a88fee4e1336647bfe0785 100644
--- a/AMDiS/src/MacroWriter.cc
+++ b/AMDiS/src/MacroWriter.cc
@@ -48,7 +48,7 @@ namespace AMDiS {
std::ofstream file;
std::list<ElementInfo> *elements = dc->getElementInfos();
- DOFVector< std::list<VertexInfo> > *vertexInfos = dc->getVertexInfos();
+ DOFVector<std::list<VertexInfo> > *vertexInfos = dc->getVertexInfos();
int dow = Global::getGeo(WORLD);
int dim = dc->getMesh()->getDim();
@@ -58,11 +58,7 @@ namespace AMDiS {
file.open(name);
- // === print file header ===
- // comment out annoying usless info
- //file << "mesh name: " << dc->getMesh()->getName() << std::endl << std::endl;
- //file << "time: " << std::scientific << time << std::endl << std::endl;
-
+ // === print file header ===
file << "DIM: " << dim << std::endl;
file << "DIM_OF_WORLD: " << dow << std::endl << std::endl;
@@ -73,7 +69,7 @@ namespace AMDiS {
file << "vertex coordinates:" << std::endl;
DOFVector< std::list<VertexInfo> >::Iterator it(vertexInfos, USED_DOFS);
- int i, counter = 0;
+ int counter = 0;
// for all DOFs
for (it.reset(); !it.end(); ++it) {
@@ -81,7 +77,7 @@ namespace AMDiS {
std::list<VertexInfo>::iterator it2;
for (it2 = it->begin(); it2 != it->end(); ++it2) {
it2->outputIndex = counter++;
- for (i = 0; i < dow; i++)
+ for (int i = 0; i < dow; i++)
file << std::scientific << it2->coords[i] << " ";
file << std::endl;
}
@@ -95,7 +91,7 @@ namespace AMDiS {
for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
// for all vertices
- for (i = 0; i < vertices; i++)
+ for (int i = 0; i < vertices; i++)
file << elementIt->vertexInfo[i]->outputIndex << " ";
file << "\n";
}
@@ -105,7 +101,7 @@ namespace AMDiS {
for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
// for all vertices
- for (i = 0; i < vertices; i++)
+ for (int i = 0; i < vertices; i++)
file << elementIt->boundary[i] << " ";
file << "\n";
}
@@ -115,7 +111,7 @@ namespace AMDiS {
for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
// for all vertices
- for (i = 0; i < vertices; i++)
+ for (int i = 0; i < vertices; i++)
file << elementIt->neighbour[i] << " ";
file << "\n";
}
@@ -125,8 +121,8 @@ namespace AMDiS {
for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
// for all vertices
- for (i = 0; i < vertices; i++) {
- if(elementIt->projection[i])
+ for (int i = 0; i < vertices; i++) {
+ if (elementIt->projection[i])
file << elementIt->projection[i]->getID() << " ";
else
file << "0 ";
@@ -144,7 +140,7 @@ namespace AMDiS {
file << std::endl << "surface region:" << std::endl;
for (elementIt = elements->begin(); elementIt != elements->end(); ++elementIt) {
- for (i = 0; i < vertices; i++)
+ for (int i = 0; i < vertices; i++)
file << elementIt->surfaceRegions[i] << " ";
file << std::endl;
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 5d2d8eb77384c14a3050f36361d7fd3cfe27f323..027fd427bb5658ee557ce2e80f6bf529008ac929 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -24,6 +24,7 @@
#include "Projection.h"
#include "ElInfoStack.h"
#include "Serializer.h"
+#include "Lagrange.h"
namespace AMDiS {
@@ -1129,32 +1130,192 @@ namespace AMDiS {
void Mesh::initialize()
{
+ FUNCNAME("Mesh::initialize()");
+
std::string macroFilename("");
std::string valueFilename("");
- std::string periodicFile("");
+ std::string periodicFilename("");
int check = 1;
GET_PARAMETER(0, name + "->macro file name", ¯oFilename);
GET_PARAMETER(0, name + "->value file name", &valueFilename);
- GET_PARAMETER(0, name + "->periodic file", &periodicFile);
+ GET_PARAMETER(0, name + "->periodic file", &periodicFilename);
GET_PARAMETER(0, name + "->check", "%d", &check);
GET_PARAMETER(0, name + "->preserve coarse dofs", "%d", &preserveCoarseDOFs);
if (macroFilename.length()) {
- macroFileInfo = MacroReader::readMacro(macroFilename.c_str(), this,
- periodicFile == "" ? NULL : periodicFile.c_str(),
- check);
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ if (checkParallelMacroFile(macroFilename, periodicFilename, check))
+ macroFileInfo =
+ MacroReader::readMacro(macroFilename + ".tmp", this,
+ periodicFilename + ".tmp", check);
+ else
+ macroFileInfo =
+ MacroReader::readMacro(macroFilename, this, periodicFilename, check);
+
+#else
+ // In sequentiel computations just read the macro file to the mesh.
+ macroFileInfo =
+ MacroReader::readMacro(macroFilename, this, periodicFilename, check);
+#endif
// If there is no value file which should be written, we can delete
// the information of the macro file.
if (!valueFilename.length())
clearMacroFileInfo();
+
}
initialized = true;
}
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ bool Mesh::checkParallelMacroFile(std::string macroFilename,
+ std::string periodicFilename,
+ int check)
+ {
+ FUNCNAME("Mesh::checkParallelMacroFile()");
+
+ // In parallel computations, first the mesh must be checked if it is refined
+ // in an apropriate way.
+
+ TEST_EXIT(admin.size() == 1)("Not yet implemented!\n");
+
+
+ // === Create a temporary mesh and load the macro file to it. ===
+
+ Mesh testMesh(name, dim);
+ DOFAdmin *localAdmin = new DOFAdmin(&testMesh, admin[0]->getName());
+ localAdmin->setNumberOfDOFs(admin[0]->getNumberOfDOFs());
+ testMesh.addDOFAdmin(localAdmin);
+
+ MacroInfo *testMacroInfo =
+ MacroReader::readMacro(macroFilename, &testMesh, periodicFilename, check);
+
+
+ // === Check the mesh structure. ===
+
+ int nMacroElements = 0;
+ int elType = -1;
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(&testMesh, -1, Mesh::CALL_LEAF_EL);
+ while (elInfo) {
+ if (elType == -1) {
+ elType = elInfo->getType();
+ } else {
+ TEST_EXIT(elType == elInfo->getType())
+ ("All elements in mesh must have the same element type!\n");
+ }
+
+ nMacroElements++;
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+
+ // === Calculate the number of global refinements such that the new mesh ===
+ // === would fulfill all requirements. ===
+
+ // There should be at least 10 macro Elements per processor, therefore:
+ // nMacroElements * 2^gr >= nProcs * 10
+ // => gr = log_2(nProcs * 10 / nMacroElements)
+
+ double tmp = 10.0 * MPI::COMM_WORLD.Get_size() / nMacroElements;
+ int nrRefine = ceil(log(tmp) / log(2));
+ if (nrRefine < 0)
+ nrRefine = 0;
+
+ if (dim == 3) {
+ int newElType = (elType + nrRefine) % 3;
+ switch (newElType) {
+ case 1:
+ if (nrRefine > 0)
+ nrRefine--;
+ else
+ nrRefine = 2;
+ break;
+ case 2:
+ nrRefine++;
+ break;
+ }
+
+ TEST_EXIT((elType + nrRefine) % 3 == 0)("This should not happen!\n");
+ }
+
+
+ // === If we do not need to refine the mesh, return back. ===
+
+ if (nrRefine == 0)
+ return false;
+
+
+ // === If macro weights are explicitly given, we cannot change the mesh. ===
+
+ std::string macroWeightsFilename = "";
+ GET_PARAMETER(0, name + "->macro weights", ¯oWeightsFilename);
+ if (macroWeightsFilename != "") {
+ ERROR_EXIT("Should not happen!\n");
+ }
+
+
+ // === Rank 0 creates a new mesh file. ===
+
+ if (MPI::COMM_WORLD.Get_rank() == 0) {
+ RefinementManager *refManager;
+ if (dim == 2)
+ refManager = new RefinementManager2d();
+ else
+ refManager = new RefinementManager3d();
+
+ refManager->globalRefine(&testMesh, nrRefine);
+ delete refManager;
+
+ Lagrange* basFcts = Lagrange::getLagrange(dim, 1);
+ FiniteElemSpace *feSpace =
+ FiniteElemSpace::provideFeSpace(localAdmin, basFcts, &testMesh, "tmp");
+
+ DataCollector dc(feSpace);
+ MacroWriter::writeMacro(&dc, (macroFilename + ".tmp").c_str());
+
+ if (periodicFilename != "")
+ MacroWriter::writePeriodicFile(&dc, (periodicFilename + ".tmp").c_str());
+ }
+
+
+ // === All ranks must wait until rank 0 has created the new macro mesh file. ===
+
+ MPI::COMM_WORLD.Barrier();
+
+
+ // === We have refined the mesh, so reduce the number of global refinements. ===
+
+ int globalRefinements = 0;
+ GET_PARAMETER(0, name + "->global refinements", "%d", &globalRefinements);
+
+ if (globalRefinements < nrRefine)
+ globalRefinements = 0;
+ else
+ globalRefinements -= nrRefine;
+
+ std::stringstream oss;
+ oss << globalRefinements;
+
+ ADD_PARAMETER(0, name + "->global refinements", oss.str().c_str());
+
+
+ // === Print a note to the screen that another mesh file will be used. ===
+
+ MSG("The macro mesh file \"%s\" was refined %d times and stored to file \"%s\".\n",
+ macroFilename.c_str(), nrRefine, (macroFilename + ".tmp").c_str());
+
+
+ return true;
+ }
+#endif
+
+
bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2)
{
std::map<BoundaryType, VertexVector*>::iterator it;
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index fe1e3f39f281a7da95f9849ce513cb46f2d718fb..1afe8995b8f58022b1e3aa5b813bfdb30b439cb2 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -673,6 +673,32 @@ namespace AMDiS {
int outside,
ElInfo *final_el_info);
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+ /** \brief
+ * This functions is called in parallel computations by the function \ref
+ * Mesh::initialize(). It checks that the macro file has enough macro elements
+ * for the number of used processors and that all macro elements are of type 0.
+ * If this is not the case, that macro mesh is globally refined in an
+ * apropriate way and is written to a new macro file.
+ *
+ * The function returns true, if a new macro and parallel file were created for
+ * the current parallel usage. In this case, a new macro file with the same name
+ * plus ".tmp", and if required, a new periodic file with the same name plus
+ * ".tmp" are created.
+ *
+ * \param[in] macroFilename Name of the macro mesh file.
+ * \param[in] periodicFilename If periodic boundaries are used, name of the
+ * periodicity file. Otherwise, the string must
+ * be empty.
+ * \param[in] check If the mesh should be checked to be a correct
+ * AMDiS macro mesh, the value must be 1 and 0
+ * otherwise.
+ */
+ bool checkParallelMacroFile(std::string macroFilename,
+ std::string periodicFilename,
+ int check);
+#endif
+
protected:
/// maximal number of DOFs at one position
static const int MAX_DOF;
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index ebafa21a28b233dc187a69d13d8eee2ddb312866..6fb4d43a6513c8d5bc4ce7d2c833783de0b9ad4e 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -196,7 +196,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor) const
{
int myRank = omp_get_thread_num();
@@ -213,7 +213,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor) const
{
int myRank = omp_get_thread_num();
@@ -230,7 +230,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor) const
{
int myRank = omp_get_thread_num();
@@ -247,7 +247,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor) const
{
int myRank = omp_get_thread_num();
@@ -304,7 +304,7 @@ namespace AMDiS {
}
/// Calls getC() for each term in \ref zeroOrder and adds the results to c.
- void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &c)
+ void getC(const ElInfo *elInfo, int nPoints, ElementVector& c)
{
int myRank = omp_get_thread_num();
diff --git a/AMDiS/src/OperatorTerm.h b/AMDiS/src/OperatorTerm.h
index 044b890b2dc57192b7635e3115d23f08d6335801..d96e0919a3cefc8f3515a51f709140f7c30882cc 100644
--- a/AMDiS/src/OperatorTerm.h
+++ b/AMDiS/src/OperatorTerm.h
@@ -103,7 +103,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor) = 0;
/** \brief
diff --git a/AMDiS/src/ParMetisPartitioner.cc b/AMDiS/src/ParMetisPartitioner.cc
index 21c473aa5a934aca44bd0588ace5f0bd691b540c..be345a3433c7d28bc9eaaf89934ec7f416f6ba33 100644
--- a/AMDiS/src/ParMetisPartitioner.cc
+++ b/AMDiS/src/ParMetisPartitioner.cc
@@ -40,7 +40,7 @@ namespace AMDiS {
nElements = elementCounter;
- TEST_EXIT(nElements > 0)("no elements in ParMETIS mesh\n");
+ TEST_EXIT(nElements > 0)("No elements in ParMETIS mesh!\n");
// allocate memory
eptr = new int[nElements + 1];
@@ -380,14 +380,14 @@ namespace AMDiS {
// update ParMETIS mesh to new partitioning
if (!parMetisMesh)
- parMetisMesh = new ParMetisMesh(mesh, mpiComm);
+ parMetisMesh = new ParMetisMesh(mesh, mpiComm);
int mpiRank = mpiComm->Get_rank();
int mpiSize = mpiComm->Get_size();
- int *nPartitionElements = new int[mpiSize];
+ std::vector<int> nPartitionElements(mpiSize);
int *elmdist = parMetisMesh->getElementDist();
- for (int i = 0; i < mpiSize; i++)
+ for (int i = 0; i < mpiSize; i++)
nPartitionElements[i] = elmdist[i + 1] - elmdist[i];
// === count number of elements ===
@@ -395,14 +395,14 @@ namespace AMDiS {
int localElements = parMetisMesh->getNumElements();
mpiComm->Allreduce(&localElements, &nElements, 1, MPI_INT, MPI_SUM);
- int *partitionElements = new int[nElements];
+ std::vector<int> partitionElements(nElements);
// distribute partition elements
mpiComm->Allgatherv(parMetisMesh->getAMDiSIndices(),
nPartitionElements[mpiRank],
MPI_INT,
- partitionElements,
- nPartitionElements,
+ &(partitionElements[0]),
+ &(nPartitionElements[0]),
elmdist,
MPI_INT);
@@ -410,9 +410,6 @@ namespace AMDiS {
for (int i = 0; i < mpiSize; i++)
for (int j = 0; j < nPartitionElements[i]; j++)
(*partitionVec)[partitionElements[elmdist[i] + j]] = i;
-
- delete [] partitionElements;
- delete [] nPartitionElements;
}
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 8a7d582a01d655e30152f50c507c69620373aefa..1d4f8866068ef64490f65a20daeb562a37b9d0f9 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -260,7 +260,8 @@ namespace AMDiS {
}
componentMeshes[i] = meshForRefinementSet[refSet];
}
- switch(dim) {
+
+ switch (dim) {
case 1:
coarseningManager = new CoarseningManager1d();
refinementManager = new RefinementManager1d();
diff --git a/AMDiS/src/ResidualEstimator.cc b/AMDiS/src/ResidualEstimator.cc
index 6ce78d7a352cd290e3e70397a984c4d00827d22c..1ebd55a460687d92710827289f5e786eeaf9a225 100644
--- a/AMDiS/src/ResidualEstimator.cc
+++ b/AMDiS/src/ResidualEstimator.cc
@@ -88,7 +88,8 @@ namespace AMDiS {
uhQP.change_dim(nPoints);
uhOldQP.change_dim(nPoints);
}
- riq = new double[nPoints];
+
+ riq.change_dim(nPoints);
grdUh_qp = NULL;
D2uhqp = NULL;
@@ -168,7 +169,6 @@ namespace AMDiS {
MSG("time estimate for component %d = %.8e\n", row, est_t_sum);
}
- delete [] riq;
delete [] basFcts;
delete [] quadFast;
@@ -260,9 +260,7 @@ namespace AMDiS {
std::vector<double*>::iterator itfac;
double det = elInfo->getDet();
double h2 = h2_from_det(det, dim);
-
- for (int iq = 0; iq < nPoints; iq++)
- riq[iq] = 0.0;
+ riq = 0.0;
for (int system = 0; system < nSystems; system++) {
if (matrix[system] == NULL)
@@ -512,7 +510,7 @@ namespace AMDiS {
DOFMatrix *A,
DOFVector<double> *fh,
Quadrature *quad,
- double *result)
+ ElementVector& result)
{
std::vector<Operator*>::iterator it;
std::vector<double*>::iterator fac;
@@ -556,7 +554,7 @@ namespace AMDiS {
if (num_rows(uhOldIq) > 0)
(*it)->evalZeroOrder(nPoints, uhOldIq, grdUhOldIq, D2UhOldIq, result, -factor);
} else {
- std::vector<double> fx(nPoints, 0.0);
+ ElementVector fx(nPoints, 0.0);
(*it)->getC(elInfo, nPoints, fx);
for (int iq = 0; iq < nPoints; iq++)
diff --git a/AMDiS/src/ResidualEstimator.h b/AMDiS/src/ResidualEstimator.h
index a4a85b92bc5faadae3fae68028b937f3403dc308..ddba9399f6a50566f6fb65e461d1db62166b37ac 100644
--- a/AMDiS/src/ResidualEstimator.h
+++ b/AMDiS/src/ResidualEstimator.h
@@ -46,7 +46,7 @@ namespace AMDiS {
DOFMatrix *A,
DOFVector<double> *fh,
Quadrature *quad,
- double *result);
+ ElementVector& result);
/// Returns pow(det,2.0/dim). Not Member of Estimator to avoid multiple instantiation.
inline double h2_from_det(double det, int dim)
@@ -148,7 +148,7 @@ namespace AMDiS {
ElementVector uhOldQP;
/// Stores the element residual computed at the quadrature points of the element.
- double *riq;
+ ElementVector riq;
WorldVector<double> *grdUh_qp;
diff --git a/AMDiS/src/RobinBC.cc b/AMDiS/src/RobinBC.cc
index 5577315edf71dfb2b4fa8327a6568ed747a5f873..ed6871bda45fe6584b48e59e548b6f5742201322 100644
--- a/AMDiS/src/RobinBC.cc
+++ b/AMDiS/src/RobinBC.cc
@@ -242,11 +242,11 @@ namespace AMDiS {
int nPoints = quadrature->getNumPoints();
mtl::dense_vector<double> uhAtQp(nPoints);
dv->getVecAtQPs(elInfo, quadrature, NULL, uhAtQp);
- std::vector<double> f(nPoints, 0.0);
+ ElementVector f(nPoints);
+ f = 0.0;
if (robinOperators)
- (*robinOperators)[face]->evalZeroOrder(nPoints, uhAtQp,
- NULL, NULL, &f[0], 1.0);
+ (*robinOperators)[face]->evalZeroOrder(nPoints, uhAtQp, NULL, NULL, f, 1.0);
std::vector<WorldVector<double> > grdUh(nPoints);
std::vector<WorldVector<double> > A_grdUh(nPoints);
@@ -265,8 +265,8 @@ namespace AMDiS {
val = 0.0;
for (int iq = 0; iq < nPoints; iq++) {
- n_A_grdUh = (normal*A_grdUh[iq]) - f[iq];
- val += quadrature->getWeight(iq)*sqr(n_A_grdUh);
+ n_A_grdUh = (normal * A_grdUh[iq]) - f[iq];
+ val += quadrature->getWeight(iq) * sqr(n_A_grdUh);
}
}
}
diff --git a/AMDiS/src/SecondOrderTerm.cc b/AMDiS/src/SecondOrderTerm.cc
index d5a52a17d1d840ade489fbb3dd0cea14543c5a9d..85ce9e7c9397bfbe00a3a4672970c7fcf54c4630 100644
--- a/AMDiS/src/SecondOrderTerm.cc
+++ b/AMDiS/src/SecondOrderTerm.cc
@@ -95,7 +95,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -135,7 +135,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -201,7 +201,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int dow = Global::getGeo(WORLD);
@@ -263,7 +263,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int dow = Global::getGeo(WORLD);
@@ -311,7 +311,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double f)
{
int dow = Global::getGeo(WORLD);
@@ -385,7 +385,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -454,7 +454,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int dow = Global::getGeo(WORLD);
@@ -517,7 +517,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -584,7 +584,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int dow = Global::getGeo(WORLD);
@@ -644,7 +644,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int dow = Global::getGeo(WORLD);
@@ -708,7 +708,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -786,7 +786,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -888,7 +888,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -1020,7 +1020,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -1108,7 +1108,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int dow = Global::getGeo(WORLD);
@@ -1158,7 +1158,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int dow = Global::getGeo(WORLD);
@@ -1207,7 +1207,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double f)
{
if (D2UhAtQP) {
@@ -1264,7 +1264,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
if (D2UhAtQP) {
diff --git a/AMDiS/src/SecondOrderTerm.h b/AMDiS/src/SecondOrderTerm.h
index 0f6b15c5c961fe44ab158ac1c673f47109b52914..deaf11b14a94590e7b182d436ea9ac4261dc3f18 100644
--- a/AMDiS/src/SecondOrderTerm.h
+++ b/AMDiS/src/SecondOrderTerm.h
@@ -131,7 +131,7 @@ namespace AMDiS {
const ElementVector&,
const WorldVector<double> *,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor)
{
int dow = Global::getGeo(WORLD);
@@ -198,7 +198,7 @@ namespace AMDiS {
const ElementVector&,
const WorldVector<double> *,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double f)
{
int dow = Global::getGeo(WORLD);
@@ -258,7 +258,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -313,7 +313,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -370,7 +370,7 @@ namespace AMDiS {
const ElementVector&,
const WorldVector<double> *,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
if (D2UhAtQP)
@@ -427,7 +427,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -471,7 +471,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -518,7 +518,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -570,7 +570,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -627,7 +627,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -677,7 +677,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -733,7 +733,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -784,7 +784,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
@@ -836,7 +836,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -886,7 +886,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -936,7 +936,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -994,7 +994,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implenetation of SecondOrderTerm::weakEval().
@@ -1054,7 +1054,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implements SecondOrderTerm::weakEval().
@@ -1104,7 +1104,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implements SecondOrderTerm::weakEval().
@@ -1158,7 +1158,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implements SecondOrderTerm::weakEval().
@@ -1204,7 +1204,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double factor);
/// Implements SecondOrderTerm::weakEval().
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index 90fc7c0c194b94f03d1e03af5040b9320c08c958..f8032eb7e36433a3579a28f1d55b7da782bf518c 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -91,7 +91,7 @@ namespace AMDiS {
const BasisFunction *phi = colFeSpace->getBasisFcts();
int nPoints = quadrature->getNumPoints();
- std::vector<double> c(nPoints, 0.0);
+ ElementVector c(nPoints, 0.0);
std::vector<double> phival(nCol);
int myRank = omp_get_thread_num();
@@ -140,7 +140,7 @@ namespace AMDiS {
void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{
int nPoints = quadrature->getNumPoints();
- std::vector<double> c(nPoints, 0.0);
+ ElementVector c(nPoints, 0.0);
int myRank = omp_get_thread_num();
std::vector<OperatorTerm*>::iterator termIt;
@@ -185,11 +185,10 @@ namespace AMDiS {
}
}
- std::vector<double> &c = tmpC[myRank];
- if (static_cast<int>(c.size()) < nPoints)
- c.resize(nPoints);
- for (int i = 0; i < nPoints; i++)
- c[i] = 0.0;
+ ElementVector &c = tmpC[myRank];
+ if (num_rows(c) < nPoints)
+ c.change_dim(nPoints);
+ c = 0.0;
for (std::vector<OperatorTerm*>::iterator termIt = terms[myRank].begin();
termIt != terms[myRank].end(); ++termIt)
@@ -244,19 +243,17 @@ namespace AMDiS {
}
}
- std::vector<double> c(nPoints, 0.0);
+ ElementVector c(nPoints, 0.0);
std::vector<OperatorTerm*>::iterator termIt;
- for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
+ for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
- }
for (int iq = 0; iq < nPoints; iq++) {
c[iq] *= elInfo->getDet();
const double *psi = psiFast->getPhi(iq);
- for (int i = 0; i < nRow; i++) {
- vec[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
- }
+ for (int i = 0; i < nRow; i++)
+ vec[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
}
}
@@ -283,7 +280,7 @@ namespace AMDiS {
}
}
- std::vector<double> c(1, 0.0);
+ ElementVector c(1, 0.0);
int myRank = omp_get_thread_num();
int size = static_cast<int>(terms[myRank].size());
@@ -329,7 +326,7 @@ namespace AMDiS {
std::vector<OperatorTerm*>::iterator termIt;
int myRank = omp_get_thread_num();
- std::vector<double> c(1, 0.0);
+ ElementVector c(1, 0.0);
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
(static_cast<ZeroOrderTerm*>(*termIt))->getC(elInfo, 1, c);
diff --git a/AMDiS/src/ZeroOrderAssembler.h b/AMDiS/src/ZeroOrderAssembler.h
index 8b349c1acbb7ad1278922dbc057fc65a7f2db297..49005137e47a65e0df2d723214b2e286e3304755 100644
--- a/AMDiS/src/ZeroOrderAssembler.h
+++ b/AMDiS/src/ZeroOrderAssembler.h
@@ -87,7 +87,7 @@ namespace AMDiS {
void calculateElementVector(const ElInfo *elInfo, ElementVector& vec);
protected:
- std::vector<std::vector<double> > tmpC;
+ std::vector<ElementVector> tmpC;
};
diff --git a/AMDiS/src/ZeroOrderTerm.cc b/AMDiS/src/ZeroOrderTerm.cc
index c8b1fd123cbd434c33246fd0d825fa3bd030e973..67e485aa63484afe54b4a613d940d84dc307bccf 100644
--- a/AMDiS/src/ZeroOrderTerm.cc
+++ b/AMDiS/src/ZeroOrderTerm.cc
@@ -32,7 +32,7 @@ namespace AMDiS {
}
- void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
if (f) {
for (int iq = 0; iq < nPoints; iq++)
@@ -48,7 +48,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
if (f) {
@@ -87,7 +87,7 @@ namespace AMDiS {
}
- void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]);
@@ -98,7 +98,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -140,7 +140,7 @@ namespace AMDiS {
getVectorAtQPs(vec2, smallElInfo, largeElInfo, subAssembler, quad, vecAtQPs2);
}
- void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
@@ -150,7 +150,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -184,7 +184,7 @@ namespace AMDiS {
getVectorAtQPs(vec3, elInfo, subAssembler, quad, vecAtQPs3);
}
- void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);
@@ -194,7 +194,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -222,7 +222,7 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void FctGradientCoords_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void FctGradientCoords_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
@@ -232,7 +232,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -261,7 +261,7 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
@@ -271,7 +271,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -301,7 +301,7 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
@@ -311,7 +311,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -345,7 +345,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -355,7 +355,7 @@ namespace AMDiS {
uhAtQP[iq];
}
- void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]);
@@ -380,7 +380,7 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
}
- void FctGradient_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void FctGradient_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(gradAtQPs[iq]);
@@ -390,7 +390,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -417,7 +417,7 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
}
- void VecAndGradAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void VecAndGradAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
@@ -427,7 +427,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -457,7 +457,7 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vecGrd, elInfo, subAssembler, quad);
}
- void VecAndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void VecAndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
@@ -467,7 +467,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -501,7 +501,7 @@ namespace AMDiS {
grad2AtQPs = getGradientsAtQPs(vecGrd2, elInfo, subAssembler, quad);
}
- void VecAndGradVec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void VecAndGradVec2AtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]);
@@ -511,7 +511,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -546,7 +546,7 @@ namespace AMDiS {
}
- void VecOfDOFVecsAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void VecOfDOFVecsAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
unsigned int size = vecs.size();
std::vector<double> arg(size);
@@ -564,7 +564,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int size = static_cast<int>(vecs.size());
@@ -603,7 +603,7 @@ namespace AMDiS {
gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);
}
- void VecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void VecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
int size = static_cast<int>(vecs.size());
std::vector<WorldVector<double>*> arg(size);
@@ -620,7 +620,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int size = static_cast<int>(vecs.size());
@@ -666,7 +666,7 @@ namespace AMDiS {
}
- void VecDivergence_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void VecDivergence_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
int size = static_cast<int>(vecs.size());
@@ -679,7 +679,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int size = static_cast<int>(vecs.size());
@@ -725,7 +725,7 @@ namespace AMDiS {
}
void VecAndVecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints,
- std::vector<double> &C)
+ ElementVector& C)
{
int size = static_cast<int>(vecs.size());
std::vector<WorldVector<double>*> arg(size);
@@ -742,7 +742,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int size = static_cast<int>(vecs.size());
@@ -786,7 +786,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -795,7 +795,7 @@ namespace AMDiS {
uhAtQP[iq];
}
- void Vec2AndGrad2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void Vec2AndGrad2AtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], gradAtQPs1[iq], gradAtQPs2[iq]);
@@ -828,7 +828,7 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vecGrd, elInfo, subAssembler, quad);
}
- void Vec2AndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void Vec2AndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs[iq]);
@@ -838,7 +838,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -887,7 +887,7 @@ namespace AMDiS {
gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
}
- void General_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void General_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
unsigned int nVecs = vecs_.size();
unsigned int nGrads = grads_.size();
@@ -906,7 +906,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
unsigned int nVecs = vecs_.size();
@@ -965,7 +965,7 @@ namespace AMDiS {
gradsAtQPs[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
}
- void GeneralParametric_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ void GeneralParametric_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
{
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
@@ -987,7 +987,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
int nVecs = static_cast<int>(vecs_.size());
@@ -1016,7 +1016,7 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void CoordsAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C)
+ void CoordsAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += (*g)(coordsAtQPs[iq]);
@@ -1026,7 +1026,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
diff --git a/AMDiS/src/ZeroOrderTerm.h b/AMDiS/src/ZeroOrderTerm.h
index cbcc557fb321eba46a3a2490f937a6dc0881874e..0a6054a7e55beed819fa2e366050f8fd9250c9a5 100644
--- a/AMDiS/src/ZeroOrderTerm.h
+++ b/AMDiS/src/ZeroOrderTerm.h
@@ -44,9 +44,7 @@ namespace AMDiS {
virtual ~ZeroOrderTerm() {}
/// Evaluates \f$ c \f$
- virtual void getC(const ElInfo *elInfo, int nPoints,
- std::vector<double> &C) = 0;
-
+ virtual void getC(const ElInfo *elInfo, int nPoints, ElementVector& C) = 0;
};
@@ -63,7 +61,7 @@ namespace AMDiS {
Simple_ZOT(double f = 1.0) : ZeroOrderTerm(0), factor(f) {}
/// Implements ZeroOrderTerm::getC().
- inline void getC(const ElInfo *, int nPoints, std::vector<double> &C)
+ inline void getC(const ElInfo *, int nPoints, ElementVector& C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += factor;
@@ -74,7 +72,7 @@ namespace AMDiS {
const ElementVector& uhAtQP,
const WorldVector<double> *,
const WorldMatrix<double> *,
- double *result,
+ ElementVector& result,
double fac)
{
for (int iq = 0; iq < nPoints; iq++)
@@ -112,14 +110,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -154,14 +152,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -204,14 +202,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -250,14 +248,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -290,14 +288,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getC().
- void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *elInfo, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -333,14 +331,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -376,14 +374,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -418,14 +416,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -462,14 +460,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements SecondOrderTerm::getC().
- void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *elInfo, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -504,14 +502,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -542,14 +540,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -584,14 +582,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -630,14 +628,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -664,14 +662,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -700,14 +698,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -732,14 +730,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -770,13 +768,13 @@ namespace AMDiS {
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -800,13 +798,13 @@ namespace AMDiS {
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -835,14 +833,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -884,14 +882,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
@@ -935,14 +933,14 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
- void getC(const ElInfo *, int nPoints, std::vector<double> &C);
+ void getC(const ElInfo *, int nPoints, ElementVector& C);
/// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const ElementVector& uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
- double *result,
+ ElementVector& result,
double fac);
protected:
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index d3e3707f5a57b1e5aa6d1a9776e35762bd0c0705..d442e6905db25ba6982f34d5a08b8d1b6742f1bb 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -71,7 +71,8 @@ namespace AMDiS {
TEST_EXIT(mpiSize > 1)
("Parallelization does not work with only one process!\n");
- TEST_EXIT(mesh)("No mesh has been defined for mesh distribution!\n");
+ TEST_EXIT(feSpace)("No FE space has been defined for the mesh distributor!\n");
+ TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
// If the problem has been already read from a file, we need only to set
// isRankDofs to all matrices and rhs vector and to remove periodic
@@ -82,14 +83,13 @@ namespace AMDiS {
return;
}
-
+
// Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
// only for macro meshes, so it will not work yet if the mesh is already refined
// in some way.
testForMacroMesh();
-
// create an initial partitioning of the mesh
partitioner->createPartitionData();
// set the element weights, which are 1 at the very first begin
@@ -97,7 +97,6 @@ namespace AMDiS {
// and now partition the mesh
partitionMesh(adaptInfo);
-
#if (DEBUG != 0)
debug::ElementIdxToDofs elMap;
debug::createSortedDofs(mesh, elMap);
@@ -112,7 +111,6 @@ namespace AMDiS {
MSG("Skip write part mesh!\n");
}
}
-
ParallelDebug::testAllElements(*this);
#endif
@@ -323,6 +321,9 @@ namespace AMDiS {
while (elInfo) {
TEST_EXIT(elInfo->getLevel() == 0)
("Mesh is already refined! This does not work with parallelization!\n");
+
+ TEST_EXIT(elInfo->getType() == 0)
+ ("Only macro elements with level 0 are supported!\n");
nMacroElements++;
@@ -476,7 +477,7 @@ namespace AMDiS {
void MeshDistributor::checkMeshChange()
{
- FUNCNAME("MeshDistributor::checkMeshChange()");
+ FUNCNAME("MeshDistributor::checkMeshChange()");
#if (DEBUG != 0)
debug::writeMesh(feSpace, -1, "before_check_mesh");