Commit 639f286f authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Global::getGeo optimized by precalculation

* added ifdef OPENMP to delete warning during compilation
parent b8cd3cef
......@@ -680,8 +680,6 @@ namespace AMDiS {
void Pre0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
double *c = GET_MEMORY(double, 1);
if (firstCall) {
q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
......@@ -691,20 +689,22 @@ namespace AMDiS {
firstCall = false;
}
c[0] = 0.0;
// c[0] = 0.0;
double c = 0.0;
int myRank = omp_get_thread_num();
int size = static_cast<int>(terms[myRank].size());
for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
(static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, c);
for (int i = 0; i < size; i++) {
(static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, &c);
}
c[0] *= elInfo->getDet();
c *= elInfo->getDet();
if (symmetric) {
for (int i = 0; i < nRow; i++) {
(*mat)[i][i] += c[0] * q00->getValue(i,i);
(*mat)[i][i] += c * q00->getValue(i,i);
for (int j = i + 1; j < nCol; j++) {
double val = c[0] * q00->getValue(i, j);
double val = c * q00->getValue(i, j);
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
......@@ -712,16 +712,12 @@ namespace AMDiS {
} else {
for (int i = 0; i < nRow; i++)
for (int j = 0; j < nCol; j++)
(*mat)[i][j] += c[0]*q00->getValue(i,j);
(*mat)[i][j] += c * q00->getValue(i, j);
}
FREE_MEMORY(c, double, 1);
}
void Pre0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
double *c = GET_MEMORY(double, 1);
if (firstCall) {
q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
......@@ -734,17 +730,15 @@ namespace AMDiS {
::std::vector<OperatorTerm*>::iterator termIt;
int myRank = omp_get_thread_num();
c[0] = 0.0;
double c = 0.0;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
(static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, c);
(static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, &c);
}
c[0] *= elInfo->getDet();
c *= elInfo->getDet();
for (int i = 0; i < nRow; i++)
(*vec)[i] += c[0] * q0->getValue(i);
FREE_MEMORY(c, double, 1);
(*vec)[i] += c * q0->getValue(i);
}
Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad)
......
......@@ -35,7 +35,7 @@ namespace AMDiS {
{
TEST_EXIT(rowFESpace)("no rowFESpace\n");
if(!colFESpace) {
if (!colFESpace) {
colFESpace = rowFESpace;
}
......@@ -57,7 +57,7 @@ namespace AMDiS {
{
FUNCNAME("DOFMatrix::~DOFMatrix()");
if(rowFESpace && rowFESpace->getAdmin()) {
if (rowFESpace && rowFESpace->getAdmin()) {
(const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
}
}
......
......@@ -35,10 +35,6 @@ namespace AMDiS {
}
}
// if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
// det = elGrdLambda(*Lambda);
// }
if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
WorldVector<double> oppC;
......@@ -221,10 +217,6 @@ namespace AMDiS {
}
}
// if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
// det = calcGrdLambda(*Lambda);
// }
if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
WorldVector<double> oppC;
......
......@@ -348,7 +348,6 @@ namespace AMDiS {
{
FUNCNAME("ElInfo3d::fillElInfo()");
int el_type_local = 0; /* el_type in {0,1,2} */
int ochild = 0; /* index of other child = 1-ichild */
int *cv = NULL; /* cv = child_vertex[el_type][ichild] */
const int (*cvg)[4] = NULL; /* cvg = child_vertex[el_type] */
......@@ -370,12 +369,12 @@ namespace AMDiS {
fillFlag_ = fillFlag__local;
parent_ = el_old;
level_ = elinfo_old->level_ + 1;
el_type = (( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->el_type + 1) % 3;
int el_type_local = ( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->getType();
el_type = (el_type_local + 1) % 3;
TEST_EXIT_DBG(element_)("missing child %d?\n", ichild);
if (fillFlag__local.isAnySet()) {
el_type_local = ( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->getType();
cvg = Tetrahedron::childVertex[el_type_local];
cv = const_cast<int*>( cvg[ichild]);
ochild = 1 - ichild;
......
......@@ -81,13 +81,6 @@ namespace AMDiS {
*/
double getNormal(int side, WorldVector<double> &normal) const;
/** \brief
* 3-dimensional realisation of ElInfo's getElementNormal method.
*/
//double getElementNormal( WorldVector<double> &normal) const;
/** \brief
* update ElInfo after refinement (of some neighbours). Only in 3d!
*/
......
......@@ -10,6 +10,7 @@
#include "Flag.h"
#include "ElInfo.h"
#include "Mesh.h"
#include "OpenMP.h"
namespace AMDiS {
......@@ -150,6 +151,10 @@ namespace AMDiS {
GET_PARAMETER(0, name + "->index decimals", "%d", &indexDecimals);
GET_PARAMETER(0, name + "->write every i-th timestep", "%d", &tsModulo);
GET_PARAMETER(0, name + "->delay", "%d", &delayWriting_);
TEST_EXIT(!delayWriting_ || amdisHaveOpenMP)
("Delayed writing only possible with OpenMP support!\n");
}
void FileWriter::writeFiles(AdaptInfo *adaptInfo,
......
......@@ -13,14 +13,18 @@ namespace AMDiS {
const char *funcName = NULL;
const char *Msg::oldFuncName = NULL;
::std::ostream* Msg::out=NULL;
::std::ostream* Msg::error=NULL;
::std::ostream* Msg::out = NULL;
::std::ostream* Msg::error = NULL;
int Global::dimOfWorld = 0;
::std::vector< ::std::vector< int > > Global::geoIndexTable;
int Msg::msgInfo = 10;
bool Msg::msgWait = true;
Element *Global::referenceElement[4] =
{ NULL, NEW Line(NULL), NEW Triangle(NULL), NEW Tetrahedron(NULL) };
{ NULL,
NEW Line(NULL),
NEW Triangle(NULL),
NEW Tetrahedron(NULL) };
void Msg::wait(bool w)
{
......@@ -55,26 +59,23 @@ namespace AMDiS {
void Msg::open_file(const char *filename, OPENMODE type)
{
FUNCNAME("Msg::open_file");
FUNCNAME("Msg::open_file()");
::std::ofstream *fp;
if (filename && (fp = new ::std::ofstream(filename, type)))
{
if (out && *out != ::std::cout && *out != ::std::cerr)
{
if (filename && (fp = new ::std::ofstream(filename, type))) {
if (out && *out != ::std::cout && *out != ::std::cerr) {
dynamic_cast< ::std::ofstream*>(out)->close();
delete out;
}
out = fp;
}
else
{
} else {
if (filename)
ERROR("can not open %s;\n", filename);
else
ERROR("no filename specified;\n");
ERROR("use previous stream for messages furthermore\n");
}
return;
......@@ -82,19 +83,16 @@ namespace AMDiS {
void Msg::change_error_out(::std::ofstream *fp)
{
FUNCNAME("Msg::change_error_out");
if (fp)
{
if (error && *error != ::std::cout && *error != ::std::cerr)
{
FUNCNAME("Msg::change_error_out()");
if (fp) {
if (error && *error != ::std::cout && *error != ::std::cerr) {
dynamic_cast< ::std::ofstream*>(error)->close();
delete error;
}
error = fp;
}
else
{
} else {
ERROR("file pointer is pointer to nil;\n");
ERROR("use previous stream for errors furthermore\n");
}
......@@ -104,27 +102,24 @@ namespace AMDiS {
void Msg::open_error_file(const char *filename, OPENMODE type)
{
FUNCNAME("Msg::open_error_file");
FUNCNAME("Msg::open_error_file()");
::std::ofstream *fp;
if (filename && (fp = new ::std::ofstream(filename, type)))
{
if (error && *error != ::std::cout && *error != ::std::cerr)
{
if (filename && (fp = new ::std::ofstream(filename, type))) {
if (error && *error != ::std::cout && *error != ::std::cerr) {
dynamic_cast< ::std::ofstream*>(error)->close();
delete error;
}
error = fp;
}
else
{
} else {
if (filename)
ERROR("can not open %s;\n", filename);
else
ERROR("no filename specified;\n");
ERROR("use previous stream for errors furthermore\n");
}
return;
}
......@@ -133,12 +128,9 @@ namespace AMDiS {
static char name[256];
char *cp;
if (path == NULL || path[0] == '\0')
{
if (path == NULL || path[0] == '\0') {
sprintf(name, "./%s", fn);
}
else
{
} else {
const char* ccp = path;
while (*ccp)
ccp++;
......@@ -159,7 +151,8 @@ namespace AMDiS {
void Msg::print_funcname(const char *funcName)
{
if (!out) out = &::std::cout;
if (!out)
out = &::std::cout;
if (funcName && oldFuncName != funcName) {
(*out)<< funcName << ":" << ::std::endl;
......@@ -176,19 +169,21 @@ namespace AMDiS {
{
static int old_line = -1;
if (!error) error = &::std::cerr;
if (!error)
error = &::std::cerr;
if (funcName && oldFuncName != funcName) {
(*error)<<funcName<< ": ";
}
else if (!funcName) {
if (line-old_line > 5) (*error)<< "*unknown function*";
if (line-old_line > 5)
(*error) << "*unknown function*";
}
if (oldFuncName != funcName) {
(*error)<<"ERROR in "<<file<<", line "<<line<<"\n";
(*error) << "ERROR in " << file << ", line " << line << "\n";
oldFuncName = funcName;
} else if (line - old_line > 5)
(*error)<< "ERROR in "<<file<<", line "<<line<<"\n";
(*error) << "ERROR in " << file << ", line " << line << "\n";
old_line = line;
}
......@@ -198,7 +193,8 @@ namespace AMDiS {
va_list arg;
char buff[255];
if (!error) error = &::std::cerr;
if (!error)
error = &::std::cerr;
va_start(arg, format);
vsprintf(buff, format, arg);
......@@ -214,7 +210,8 @@ namespace AMDiS {
char buff[255];
if (!error) error = &::std::cerr;
if (!error)
error = &::std::cerr;
va_start(arg, format);
vsprintf(buff, format, arg);
......@@ -230,19 +227,20 @@ namespace AMDiS {
{
static int old_line = -1;
if (!out) out = &::std::cout;
if (!out)
out = &::std::cout;
if (funcName && oldFuncName != funcName) {
(*out)<<funcName<<": ";
(*out) << funcName << ": ";
}
else if (!funcName){
(*out)<< "*unknown function*";
(*out) << "*unknown function*";
}
if (oldFuncName != funcName) {
(*out)<<"WARNING in "<<file<<", line "<<line<<"\n";
(*out) << "WARNING in " << file << ", line " << line << "\n";
oldFuncName = funcName;
} else if (line - old_line > 5)
(*out)<<"WARNING in "<<file<<", line "<<line<<"\n";
(*out) << "WARNING in " << file << ", line " << line << "\n";
old_line = line;
}
......@@ -252,7 +250,8 @@ namespace AMDiS {
va_list arg;
char buff[255];
if (!out) out = &::std::cout;
if (!out)
out = &::std::cout;
va_start(arg, format);
vsprintf(buff, format, arg);
......@@ -273,7 +272,7 @@ namespace AMDiS {
va_start(arg, format);
vsprintf(buff, format, arg);
(*out)<<buff;
(*out) << buff;
va_end(arg);
}
......@@ -283,54 +282,53 @@ namespace AMDiS {
// get dimension
TEST_EXIT(Parameters::initialized())("Parameters not initialized!\n");
Parameters::getGlobalParameter(0,"dimension of world","%d",&d);
TEST_EXIT(d > 0)("can not initialize dimension\n");
Parameters::getGlobalParameter(0, "dimension of world","%d",&d);
TEST_EXIT(d > 0)("Cannot initialize dimension!\n");
TEST_EXIT((d == 1) || (d == 2) || (d == 3))("Invalid world dimension %d!\n",d);
// set dimension
dimOfWorld = d;
// set msgWait
Parameters::getGlobalParameter(0,"WAIT","%d",&d);
Msg::setMsgWait(!(d==0));
};
int Global::getGeo(GeoIndex p, int dim)
{
FUNCNAME("Global::getGeo()");
initTest();
TEST_EXIT_DBG((p >= MINPART) && (p <= MAXPART))
("Calling for invalid geometry value %d\n",p);
TEST_EXIT_DBG((dim >= 0) && (dim < 4))
("invalid dim: %d\n", dim);
if (p == WORLD)
return dimOfWorld;
if (dim == 0) {
switch(p) {
case PARTS:
return 1;
case VERTEX:
return 1;
case EDGE:
return 0;
case FACE:
return 0;
default:
ERROR_EXIT("dim = 0\n");
// prepare geoIndex-Table
int geoTableSize = abs(static_cast<int>(MINPART)) + MAXPART + 1;
geoIndexTable.resize(4);
for (int i = 0; i < 4; i++) {
geoIndexTable[i].resize(geoTableSize);
for (int j = 0; j < geoTableSize; j++) {
geoIndexTable[i][j] = 0;
}
}
return referenceElement[dim]->getGeo(p);
geoIndexTable[0][PARTS - MINPART] = 1;
geoIndexTable[0][VERTEX - MINPART] = 1;
geoIndexTable[0][EDGE - MINPART] = 0;
geoIndexTable[0][FACE - MINPART] = 0;
geoIndexTable[0][WORLD - MINPART] = dimOfWorld;
for (int i = 1; i < 4; i++) {
geoIndexTable[i][CENTER - MINPART] = referenceElement[i]->getGeo(CENTER);
geoIndexTable[i][VERTEX - MINPART] = referenceElement[i]->getGeo(VERTEX);
geoIndexTable[i][EDGE - MINPART] = referenceElement[i]->getGeo(EDGE);
geoIndexTable[i][FACE - MINPART] = referenceElement[i]->getGeo(FACE);
geoIndexTable[i][DIMEN - MINPART] = referenceElement[i]->getGeo(DIMEN);
geoIndexTable[i][PARTS - MINPART] = referenceElement[i]->getGeo(PARTS);
geoIndexTable[i][NEIGH - MINPART] = referenceElement[i]->getGeo(NEIGH);
geoIndexTable[i][WORLD - MINPART] = dimOfWorld;
geoIndexTable[i][BOUNDARY - MINPART] = referenceElement[i]->getGeo(BOUNDARY);
geoIndexTable[i][PROJECTION - MINPART] = referenceElement[i]->getGeo(PROJECTION);
}
// set msgWait
Parameters::getGlobalParameter(0, "WAIT", "%d", &d);
Msg::setMsgWait(!(d == 0));
};
int fac(int i)
{
if(i<=1) return 1;
else return i*fac(i-1);
if (i <= 1)
return 1;
else
return i * fac(i - 1);
}
::std::string memSizeStr(int size)
......
......@@ -453,8 +453,6 @@ namespace AMDiS {
* returns geometrical information. Currently this is only dimOfWorld.
*/
static inline int getGeo(GeoIndex p) {
FUNCNAME("Global::getGeo()");
initTest();
if (WORLD == p)
return dimOfWorld;
......@@ -467,7 +465,17 @@ namespace AMDiS {
* returns geometrical information about elements of the dimension dim.
* getGeo(VERTEX, 3) returns 4 because a Tetrahedron has 4 vertices.
*/
static int getGeo(GeoIndex p, int dim);
static inline int getGeo(GeoIndex p, int dim) {
initTest();
TEST_EXIT_DBG((p >= MINPART) && (p <= MAXPART))
("Calling for invalid geometry value %d\n",p);
TEST_EXIT_DBG((dim >= 0) && (dim < 4))
("invalid dim: %d\n", dim);
TEST_EXIT_DBG((dim != 0) || (p == PARTS || p == VERTEX || p == EDGE || p == FACE))
("dim = 0\n");
return geoIndexTable[dim][p - MINPART];
}
private:
/** \brief
......@@ -503,6 +511,8 @@ namespace AMDiS {
* Tetrahedron wich is 4 => no switch statement is necessary.
*/
static Element *referenceElement[4];
static ::std::vector< ::std::vector< int > > geoIndexTable;
};
#define COMMA ,
......
......@@ -792,7 +792,7 @@ namespace AMDiS {
DimVec<int> parts(dim, NO_INIT);
for (int i = 0; i < dim + 1; i++)
parts[i] = Global::getGeo(INDEX_OF_DIM(i, dim),dim);
parts[i] = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
int index = 0;
......
......@@ -617,15 +617,13 @@ namespace AMDiS {
static const MacroElement *mel = NULL;
DimVec<double> lambda(dim, NO_INIT);
ElInfo *mel_info = NULL;
int i, k;
bool inside;
mel_info = createNewElInfo();
if (start_mel != NULL)
mel = start_mel;
else
if((mel == NULL)||(mel->getElement()->getMesh() != this))
if ((mel == NULL) || (mel->getElement()->getMesh() != this))
mel = *(macroElements.begin());
mel_info->setFillFlag(Mesh::FILL_COORDS);
......@@ -635,7 +633,7 @@ namespace AMDiS {
mel_info->fillMacroInfo(mel);
int k;
while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
if (mel->getNeighbour(k)) {
mel = mel->getNeighbour(k);
......@@ -646,8 +644,10 @@ namespace AMDiS {
}
/* now, descend in tree to find leaf element at point */
inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
for (i=0; i<=dim; i++) bary[i] = final_lambda[i];
bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
for (int i = 0; i <= dim; i++) {
bary[i] = final_lambda[i];
}
DELETE mel_info;
......@@ -661,12 +661,8 @@ namespace AMDiS {
const WorldVector<double> *xy0,
double *sp)
{
ElInfo *el_info = NULL;
int val;
el_info = createNewElInfo();
val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
ElInfo *el_info = createNewElInfo();
int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
*elp = el_info->getElement();
......@@ -682,32 +678,31 @@ namespace AMDiS {
int outside,
ElInfo* final_el_info)
{
FUNCNAME("Mesh::findElementAtPointRecursive");
FUNCNAME("Mesh::findElementAtPointRecursive()");
Element *el = el_info->getElement();
ElInfo *c_el_info = NULL;
DimVec<double> c_lambda(dim, NO_INIT);
int i, inside;
int inside;
int ichild, c_outside;
if (el->isLeaf()) {
*final_el_info = *el_info;
if (outside < 0) {
for (i=0; i<=dim; i++) final_lambda[i] = lambda[i];
return(true);
for (int i = 0; i <= dim; i++) {
final_lambda[i] = lambda[i];
}
else
{ /* outside */
if (g_xy0)
{ /* find boundary point of [xy0, xy] */
double s;