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

* Some more optimizations

parent c7cbb878
...@@ -239,45 +239,23 @@ namespace AMDiS { ...@@ -239,45 +239,23 @@ namespace AMDiS {
for (int i = 0; i < n_row; i++) { // for all rows of element matrix for (int i = 0; i < n_row; i++) { // for all rows of element matrix
row = elMat.rowIndices[i]; row = elMat.rowIndices[i];
// if (omp_get_thread_num() == 0) {
// ::std::cout << "bound[i] = " << bound[i] << ::std::endl;
// }
BoundaryCondition *condition = BoundaryCondition *condition =
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL; bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
if (condition && condition->isDirichlet()) { if (condition && condition->isDirichlet()) {
// if (omp_get_thread_num() == 0) {
// ::std::cout << "IN BC" << ::std::endl;
// }
MatrixRow *matrixRow = &(matrix[row]); MatrixRow *matrixRow = &(matrix[row]);
if (coupleMatrix) { if (coupleMatrix) {
// if (omp_get_thread_num() == 0) {
// ::std::cout << "RESIZE 0" << ::std::endl;
// }
matrixRow->resize(0); matrixRow->resize(0);
} else { } else {
// if (omp_get_thread_num() == 0) {
// ::std::cout << "RESIZE 1 with col = " << row << ::std::endl;
// }
matrixRow->resize(1); matrixRow->resize(1);
((*matrixRow)[0]).col = row; ((*matrixRow)[0]).col = row;
((*matrixRow)[0]).entry = 1.0; ((*matrixRow)[0]).entry = 1.0;
} }
} else { } else {
// if (omp_get_thread_num() == 0) {
// ::std::cout << "OUT BC" << ::std::endl;
// }
for (int j = 0; j < n_col; j++) { // for all columns for (int j = 0; j < n_col; j++) { // for all columns
col = elMat.colIndices[j]; col = elMat.colIndices[j];
entry = elMat[i][j]; entry = elMat[i][j];
addSparseDOFEntry(sign, row, col, entry, add); addSparseDOFEntry(sign, row, col, entry, add);
// if (omp_get_thread_num() == 0) {
// ::std::cout << "addSparseDOFEntry(" << sign << "," << row << "," << col << "," << entry << "," << add << ")" << ::std::endl;
// }
} }
} }
} }
...@@ -340,30 +318,27 @@ namespace AMDiS { ...@@ -340,30 +318,27 @@ namespace AMDiS {
("Column index %d out of range 0-%d\n", jcol, colFESpace->getAdmin()->getUsedSize() - 1); ("Column index %d out of range 0-%d\n", jcol, colFESpace->getAdmin()->getUsedSize() - 1);
// first entry is diagonal entry // first entry is diagonal entry
if(rowFESpace==colFESpace) if (rowFESpace==colFESpace)
if (rowSize == 0) { if (rowSize == 0) {
MatEntry newEntry = {irow, 0.0}; MatEntry newEntry = {irow, 0.0};
row->push_back(newEntry); row->push_back(newEntry);
rowSize = 1; rowSize = 1;
// if (omp_get_thread_num() == 0) {
// ::std::cout << "PB: " << irow << " " << 0.0 << ::std::endl;
// }
} }
// search jcol // search jcol
for(i=0; i < rowSize; i++) { for (i = 0; i < rowSize; i++) {
// jcol found ? // jcol found ?
if((*row)[i].col == jcol) { if ((*row)[i].col == jcol) {
break; break;
} }
// remember free entry // remember free entry
if((*row)[i].col == UNUSED_ENTRY) { if ((*row)[i].col == UNUSED_ENTRY) {
freeCol = i; freeCol = i;
} }
// no more entries // no more entries
if((*row)[i].col == NO_MORE_ENTRIES) { if ((*row)[i].col == NO_MORE_ENTRIES) {
freeCol = i; freeCol = i;
if(rowSize > i+1) { if (rowSize > i+1) {
(*row)[i+1].entry = NO_MORE_ENTRIES; (*row)[i+1].entry = NO_MORE_ENTRIES;
} }
break; break;
...@@ -371,29 +346,22 @@ namespace AMDiS { ...@@ -371,29 +346,22 @@ namespace AMDiS {
} }
// jcol found? // jcol found?
if(i < rowSize) { if (i < rowSize) {
if(!add) (*row)[i].entry = 0.0; if (!add)
(*row)[i].entry = 0.0;
(*row)[i].entry += sign * entry; (*row)[i].entry += sign * entry;
result = &((*row)[i].entry); result = &((*row)[i].entry);
if (omp_get_thread_num() == 0) {
// ::std::cout << "ADD: " << i << " " << sign * entry << ::std::endl;
}
} else { } else {
if(freeCol == -1) { if(freeCol == -1) {
MatEntry newEntry = {jcol, sign * entry}; MatEntry newEntry = {jcol, sign * entry};
row->push_back(newEntry); row->push_back(newEntry);
result = &((*row)[row->size() - 1].entry); result = &((*row)[row->size() - 1].entry);
if (omp_get_thread_num() == 0) {
// ::std::cout << "PB: " << jcol << " " << sign * entry << ::std::endl;
}
} else { } else {
(*row)[freeCol].col = jcol; (*row)[freeCol].col = jcol;
if(!add) (*row)[freeCol].entry = 0.0; if (!add)
(*row)[freeCol].entry = 0.0;
(*row)[freeCol].entry += sign * entry; (*row)[freeCol].entry += sign * entry;
result = &((*row)[freeCol].entry); result = &((*row)[freeCol].entry);
if (omp_get_thread_num() == 0) {
// ::std::cout << "ADD: " << jcol << " " << sign * entry << ::std::endl;
}
} }
} }
......
...@@ -38,6 +38,8 @@ namespace AMDiS { ...@@ -38,6 +38,8 @@ namespace AMDiS {
for (int i = 0; i < neighbourCoord_.getSize(); i++) { for (int i = 0; i < neighbourCoord_.getSize(); i++) {
neighbourCoord_[i].init(mesh_->getDim()); neighbourCoord_[i].init(mesh_->getDim());
} }
dimOfWorld = Global::getGeo(WORLD);
} }
ElInfo::~ElInfo() ElInfo::~ElInfo()
...@@ -53,33 +55,22 @@ namespace AMDiS { ...@@ -53,33 +55,22 @@ namespace AMDiS {
const WorldVector<double> *ElInfo::coordToWorld(const DimVec<double>& l, const WorldVector<double> *ElInfo::coordToWorld(const DimVec<double>& l,
WorldVector<double>* w) const WorldVector<double>* w) const
{
return coordToWorld(l, &coord_, w);
}
const WorldVector<double> *ElInfo::coordToWorld(const DimVec<double>& l,
const FixVec<WorldVector<double>, VERTEX> *coords,
WorldVector<double> *w)
{ {
int dim = l.getSize() - 1; int dim = l.getSize() - 1;
int dimOfWorld = Global::getGeo(WORLD);
static WorldVector<double> world[2]; static WorldVector<double> world[2];
WorldVector<double> *ret = w ? w : &world[omp_get_thread_num()]; WorldVector<double> *ret = w ? w : &world[omp_get_thread_num()];
WorldVector<double> v = (*coords)[0];
double c = l[0]; double c = l[0];
for (int j = 0; j < dimOfWorld; j++) for (int j = 0; j < dimOfWorld; j++)
(*ret)[j] = c * v[j]; (*ret)[j] = c * coord_[0][j];
int vertices = Global::getGeo(VERTEX, dim); int vertices = Global::getGeo(VERTEX, dim);
for (int i = 1; i < vertices; i++) { for (int i = 1; i < vertices; i++) {
v = (*coords)[i];
c = l[i]; c = l[i];
for (int j = 0; j < dimOfWorld; j++) for (int j = 0; j < dimOfWorld; j++)
(*ret)[j] += c * v[j]; (*ret)[j] += c * coord_[i][j];
} }
return ret; return ret;
} }
...@@ -121,7 +112,7 @@ namespace AMDiS { ...@@ -121,7 +112,7 @@ namespace AMDiS {
e3[i] = coords[3][i] - v0[i]; e3[i] = coords[3][i] - v0[i];
} }
switch(dow) { switch (dow) {
case 1: case 1:
det = coords[1][0] - coords[0][0]; det = coords[1][0] - coords[0][0];
break; break;
...@@ -137,9 +128,9 @@ namespace AMDiS { ...@@ -137,9 +128,9 @@ namespace AMDiS {
WorldVector<double> n; WorldVector<double> n;
if (dim > 1) { if (dim > 1) {
n[0] = e1[1]*e2[2] - e1[2]*e2[1]; n[0] = e1[1] * e2[2] - e1[2] * e2[1];
n[1] = e1[2]*e2[0] - e1[0]*e2[2]; n[1] = e1[2] * e2[0] - e1[0] * e2[2];
n[2] = e1[0]*e2[1] - e1[1]*e2[0]; n[2] = e1[0] * e2[1] - e1[1] * e2[0];
} }
if (dim == 1) { if (dim == 1) {
...@@ -149,11 +140,11 @@ namespace AMDiS { ...@@ -149,11 +140,11 @@ namespace AMDiS {
} else if (dim == 3) { } else if (dim == 3) {
det = n[0] * e3[0] + n[1] * e3[1] + n[2] * e3[2]; det = n[0] * e3[0] + n[1] * e3[1] + n[2] * e3[2];
} else } else
ERROR_EXIT("not yet for problem dimension = %d",dim); ERROR_EXIT("not yet for problem dimension = %d", dim);
break; break;
} }
default: default:
ERROR_EXIT("not yet for Global::getGeo(WORLD) = %d",dow); ERROR_EXIT("not yet for Global::getGeo(WORLD) = %d", dow);
} }
return abs(det); return abs(det);
......
...@@ -441,24 +441,6 @@ namespace AMDiS { ...@@ -441,24 +441,6 @@ namespace AMDiS {
}; };
// ===== protected methods ====================================================
protected:
/** \brief
* Used by ElInfo's coordToWorld().
* \param lambda barycentric coordinates which should be converted.
* \param coords coords of element vertices.
* \param world world coordinates corresponding to lambda
*/
static const WorldVector<double>
*coordToWorld(const DimVec<double>& lambda,
const FixVec<WorldVector<double>, VERTEX> *coords,
WorldVector<double> *world);
// ===== protected members ====================================================
protected: protected:
/** \brief /** \brief
...@@ -555,6 +537,11 @@ namespace AMDiS { ...@@ -555,6 +537,11 @@ namespace AMDiS {
*/ */
bool parametric_; bool parametric_;
/** \brief
* Stores the world dimension.
*/
int dimOfWorld;
// ===== static public members ================================================ // ===== static public members ================================================
public: public:
/** \brief /** \brief
......
...@@ -124,7 +124,7 @@ namespace AMDiS { ...@@ -124,7 +124,7 @@ namespace AMDiS {
TEST_EXIT_DBG(lambda)("lambda must not be NULL\n"); TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
TEST_EXIT_DBG(dim == 1)("dim!=1\n"); TEST_EXIT_DBG(dim == 1)("dim!=1\n");
TEST_EXIT_DBG(Global::getGeo(WORLD) == dim)("not yet for DIM != DIM_OF_WORLD\n"); TEST_EXIT_DBG(dimOfWorld == dim)("not yet for DIM != DIM_OF_WORLD\n");
if (abs(length) < DBL_TOL) { if (abs(length) < DBL_TOL) {
ERROR_EXIT("length = %le; abort\n", length); ERROR_EXIT("length = %le; abort\n", length);
...@@ -173,15 +173,13 @@ namespace AMDiS { ...@@ -173,15 +173,13 @@ namespace AMDiS {
{ {
FUNCNAME("ElInfo::getElementNormal()"); FUNCNAME("ElInfo::getElementNormal()");
double det = 0.0; TEST_EXIT_DBG(dimOfWorld == 2)
int dow = Global::getGeo(WORLD); (" element normal only well defined for DIM_OF_WORLD = DIM + 1 !!");
TEST_EXIT_DBG(dow == 2)(" element normal only well defined for DIM_OF_WORLD = DIM + 1 !!");
elementNormal[0] = coord_[1][1] - coord_[0][1]; elementNormal[0] = coord_[1][1] - coord_[0][1];
elementNormal[1] = coord_[0][0] - coord_[1][0]; elementNormal[1] = coord_[0][0] - coord_[1][0];
det = norm(&elementNormal); double det = norm(&elementNormal);
TEST_EXIT_DBG(det > 1.e-30)("det = 0"); TEST_EXIT_DBG(det > 1.e-30)("det = 0");
......
...@@ -148,8 +148,6 @@ namespace AMDiS { ...@@ -148,8 +148,6 @@ namespace AMDiS {
parent_ = elem; parent_ = elem;
level_ = elinfo_old->level_ + 1; level_ = elinfo_old->level_ + 1;
int dow = Global::getGeo(WORLD);
if (fillFlag_.isSet(Mesh::FILL_COORDS) || if (fillFlag_.isSet(Mesh::FILL_COORDS) ||
fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_DET) ||
fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) { fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
...@@ -402,22 +400,21 @@ namespace AMDiS { ...@@ -402,22 +400,21 @@ namespace AMDiS {
const WorldVector<double> v0 = coord_[0]; const WorldVector<double> v0 = coord_[0];
testFlag(Mesh::FILL_COORDS); testFlag(Mesh::FILL_COORDS);
int dow = Global::getGeo(WORLD);
int dim = mesh_->getDim(); int dim = mesh_->getDim();
for (int i = 0; i < dow; i++) { for (int i = 0; i < dimOfWorld; i++) {
e1[i] = coord_[1][i] - v0[i]; e1[i] = coord_[1][i] - v0[i];
e2[i] = coord_[2][i] - v0[i]; e2[i] = coord_[2][i] - v0[i];
} }
if (dow == 2) { if (dimOfWorld == 2) {
double sdet = e1[0] * e2[1] - e1[1] * e2[0]; double sdet = e1[0] * e2[1] - e1[1] * e2[0];
adet = abs(sdet); adet = abs(sdet);
if (adet < 1.0E-25) { if (adet < 1.0E-25) {
MSG("abs(det) = %f\n", adet); MSG("abs(det) = %f\n", adet);
for (int i = 0; i <= dim; i++) for (int i = 0; i <= dim; i++)
for (int j = 0; j < dow; j++) for (int j = 0; j < dimOfWorld; j++)
grd_lam[i][j] = 0.0; grd_lam[i][j] = 0.0;
} else { } else {
double det1 = 1.0 / sdet; double det1 = 1.0 / sdet;
...@@ -443,7 +440,7 @@ namespace AMDiS { ...@@ -443,7 +440,7 @@ namespace AMDiS {
if (adet < 1.0E-15) { if (adet < 1.0E-15) {
MSG("abs(det) = %lf\n", adet); MSG("abs(det) = %lf\n", adet);
for (int i = 0; i <= dim; i++) for (int i = 0; i <= dim; i++)
for (int j = 0; j < dow; j++) for (int j = 0; j < dimOfWorld; j++)
grd_lam[i][j] = 0.0; grd_lam[i][j] = 0.0;
} else { } else {
vectorProduct(e2, normal, grd_lam[1]); vectorProduct(e2, normal, grd_lam[1]);
...@@ -451,7 +448,7 @@ namespace AMDiS { ...@@ -451,7 +448,7 @@ namespace AMDiS {
double adet2 = 1.0 / (adet * adet); double adet2 = 1.0 / (adet * adet);
for (int i = 0; i < dow; i++) { for (int i = 0; i < dimOfWorld; i++) {
grd_lam[1][i] *= adet2; grd_lam[1][i] *= adet2;
grd_lam[2][i] *= adet2; grd_lam[2][i] *= adet2;
} }
...@@ -477,7 +474,6 @@ namespace AMDiS { ...@@ -477,7 +474,6 @@ namespace AMDiS {
static DimVec<double> vec(mesh_->getDim(), NO_INIT); static DimVec<double> vec(mesh_->getDim(), NO_INIT);
int dim = mesh_->getDim(); int dim = mesh_->getDim();
int dimOfWorld = Global::getGeo(WORLD);
for (int j = 0; j < dimOfWorld; j++) { for (int j = 0; j < dimOfWorld; j++) {
double x0 = coord_[dim][j]; double x0 = coord_[dim][j];
...@@ -523,7 +519,7 @@ namespace AMDiS { ...@@ -523,7 +519,7 @@ namespace AMDiS {
int i0 = (side + 1) % 3; int i0 = (side + 1) % 3;
int i1 = (side + 2) % 3; int i1 = (side + 2) % 3;
if (Global::getGeo(WORLD) == 2){ if (dimOfWorld == 2){
normal[0] = coord_[i1][1] - coord_[i0][1]; normal[0] = coord_[i1][1] - coord_[i0][1];
normal[1] = coord_[i0][0] - coord_[i1][0]; normal[1] = coord_[i0][0] - coord_[i1][0];
} else { // dow == 3 } else { // dow == 3
...@@ -559,7 +555,8 @@ namespace AMDiS { ...@@ -559,7 +555,8 @@ namespace AMDiS {
{ {
FUNCNAME("ElInfo::getElementNormal()"); FUNCNAME("ElInfo::getElementNormal()");
TEST_EXIT_DBG(Global::getGeo(WORLD) == 3)(" element normal only well defined for DIM_OF_WORLD = DIM + 1 !!"); TEST_EXIT_DBG(dimOfWorld == 3)
(" element normal only well defined for DIM_OF_WORLD = DIM + 1 !!");
WorldVector<double> e0 = coord_[1] - coord_[0]; WorldVector<double> e0 = coord_[1] - coord_[0];
WorldVector<double> e1 = coord_[2] - coord_[0]; WorldVector<double> e1 = coord_[2] - coord_[0];
......
...@@ -111,7 +111,7 @@ namespace AMDiS { ...@@ -111,7 +111,7 @@ namespace AMDiS {
{ {
FUNCNAME("ElInfo3d::calcGrdLambda()"); FUNCNAME("ElInfo3d::calcGrdLambda()");
TEST_EXIT_DBG(Global::getGeo(WORLD) == 3) TEST_EXIT_DBG(dimOfWorld == 3)
("dim != dim_of_world ! use parametric elements!\n"); ("dim != dim_of_world ! use parametric elements!\n");
WorldVector<double> e1, e2, e3; WorldVector<double> e1, e2, e3;
...@@ -183,7 +183,6 @@ namespace AMDiS { ...@@ -183,7 +183,6 @@ namespace AMDiS {
TEST_EXIT_DBG(lambda)("lambda must not be NULL\n"); TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");
int dim = mesh_->getDim(); int dim = mesh_->getDim();
int dimOfWorld = Global::getGeo(WORLD);
TEST_EXIT_DBG(dim == dimOfWorld)("dim!=dimOfWorld not yet implemented\n"); TEST_EXIT_DBG(dim == dimOfWorld)("dim!=dimOfWorld not yet implemented\n");
...@@ -267,7 +266,6 @@ namespace AMDiS { ...@@ -267,7 +266,6 @@ namespace AMDiS {
int neighbours = mesh_->getGeo(NEIGH); int neighbours = mesh_->getGeo(NEIGH);
int vertices = mesh_->getGeo(VERTEX); int vertices = mesh_->getGeo(VERTEX);
int dow = Global::getGeo(WORLD);
if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
Tetrahedron *nb; Tetrahedron *nb;
...@@ -295,7 +293,7 @@ namespace AMDiS { ...@@ -295,7 +293,7 @@ namespace AMDiS {
if (nb->isNewCoordSet()) if (nb->isNewCoordSet())
oppCoord_[ineigh] = *(nb->getNewCoord()); oppCoord_[ineigh] = *(nb->getNewCoord());
else else
for (int j = 0; j < dow; j++) for (int j = 0; j < dimOfWorld; j++)
oppCoord_[ineigh][j] = (oppCoord_[ineigh][j] + coord_[k][j]) / 2; oppCoord_[ineigh][j] = (oppCoord_[ineigh][j] + coord_[k][j]) / 2;
} }
neighbour_[ineigh] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(nb->getChild(1-ov))); neighbour_[ineigh] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(nb->getChild(1-ov)));
...@@ -313,14 +311,12 @@ namespace AMDiS { ...@@ -313,14 +311,12 @@ namespace AMDiS {
double det = 0.0; double det = 0.0;
WorldVector<double> e0, e1, e2; WorldVector<double> e0, e1, e2;
int dow = Global::getGeo(WORLD); if (dimOfWorld == 3) {
if (dow == 3) {
int i0 = (face + 1) % 4; int i0 = (face + 1) % 4;
int i1 = (face + 2) % 4; int i1 = (face + 2) % 4;
int i2 = (face + 3) % 4; int i2 = (face + 3) % 4;
for (int i = 0; i < dow; i++) { for (int i = 0; i < dimOfWorld; i++) {
e0[i] = coord_[i1][i] - coord_[i0][i]; e0[i] = coord_[i1][i] - coord_[i0][i];
e1[i] = coord_[i2][i] - coord_[i0][i]; e1[i] = coord_[i2][i] - coord_[i0][i];
e2[i] = coord_[face][i] - coord_[i0][i]; e2[i] = coord_[face][i] - coord_[i0][i];
...@@ -329,7 +325,7 @@ namespace AMDiS { ...@@ -329,7 +325,7 @@ namespace AMDiS {
vectorProduct(e0, e1, normal); vectorProduct(e0, e1, normal);
if ((e2*normal) < 0.0) if ((e2*normal) < 0.0)
for (int i = 0; i < dow; i++) for (int i = 0; i < dimOfWorld; i++)
normal[i] = -normal[i]; normal[i] = -normal[i];
det = norm(&normal); det = norm(&normal);
...@@ -339,7 +335,7 @@ namespace AMDiS { ...@@ -339,7 +335,7 @@ namespace AMDiS {
normal[1] /= det; normal[1] /= det;
normal[2] /= det; normal[2] /= det;
} else { } else {
MSG("not implemented for DIM_OF_WORLD = %d in 3d\n", dow); MSG("not implemented for DIM_OF_WORLD = %d in 3d\n", dimOfWorld);
} }
return(det); return(det);
...@@ -367,7 +363,7 @@ namespace AMDiS { ...@@ -367,7 +363,7 @@ namespace AMDiS {
Flag fill_opp_coords; Flag fill_opp_coords;
Mesh *mesh = elinfo_old->getMesh(); Mesh *mesh = elinfo_old->getMesh();
TEST_EXIT_DBG(el_old->getChild(0))("missing child?\n"); /* Kuni 22.08.96 */ TEST_EXIT_DBG(el_old->getChild(0))("missing child?\n");
element_ = const_cast<Element*>( el_old->getChild(ichild)); element_ = const_cast<Element*>( el_old->getChild(ichild));
macroElement_ = elinfo_old->macroElement_; macroElement_ = elinfo_old->macroElement_;
...@@ -385,20 +381,18 @@ namespace AMDiS { ...@@ -385,20 +381,18 @@ namespace AMDiS {