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

* Global refinement bug fix (for different meshes of different components)

parent 23e3a357
...@@ -279,7 +279,8 @@ namespace AMDiS { ...@@ -279,7 +279,8 @@ namespace AMDiS {
} }
} }
} else { } else {
DimVec<double> grdPhi(dim, DEFAULT_VALUE, 0.0); // DimVec<double> grdPhi(dim, DEFAULT_VALUE, 0.0);
DimVec<double>* grdPhi = grdPhis[omp_get_thread_num()];
for (int i = 0; i < numPoints; i++) { for (int i = 0; i < numPoints; i++) {
for (int j = 0; j < dim + 1; j++) { for (int j = 0; j < dim + 1; j++) {
...@@ -287,9 +288,9 @@ namespace AMDiS { ...@@ -287,9 +288,9 @@ namespace AMDiS {
} }
for (int j = 0; j < nBasFcts; j++) { for (int j = 0; j < nBasFcts; j++) {
(*(basFcts->getGrdPhi(j)))(quad->getLambda(i), grdPhi); (*(basFcts->getGrdPhi(j)))(quad->getLambda(i), *grdPhi);
for (int k = 0; k < parts; k++) { for (int k = 0; k < parts; k++) {
grd1[k] += grdPhi[k] * localVec[j]; grd1[k] += (*grdPhi)[k] * localVec[j];
} }
} }
...@@ -379,7 +380,9 @@ namespace AMDiS { ...@@ -379,7 +380,9 @@ namespace AMDiS {
} }
} }
} else { } else {
DimMat<double> D2Phi(dim, NO_INIT); // DimMat<double> D2Phi(dim, NO_INIT);
DimMat<double>* D2Phi = D2Phis[omp_get_thread_num()];
for (iq = 0; iq < numPoints; iq++) { for (iq = 0; iq < numPoints; iq++) {
for (k = 0; k < parts; k++) for (k = 0; k < parts; k++)
for (l = 0; l < parts; l++) for (l = 0; l < parts; l++)
...@@ -387,11 +390,11 @@ namespace AMDiS { ...@@ -387,11 +390,11 @@ namespace AMDiS {
for (i = 0; i < nBasFcts; i++) { for (i = 0; i < nBasFcts; i++) {
WARNING("not tested after index correction\n"); WARNING("not tested after index correction\n");
(*(basFcts->getD2Phi(i)))(quad->getLambda(iq), D2Phi); (*(basFcts->getD2Phi(i)))(quad->getLambda(iq), *D2Phi);
for (k = 0; k < parts; k++) for (k = 0; k < parts; k++)
for (l = 0; l < parts; l++) for (l = 0; l < parts; l++)
D2Tmp[k][l] += localVec[i] * D2Phi[k][l]; D2Tmp[k][l] += localVec[i] * (*D2Phi)[k][l];
} }
for (i = 0; i < dow; i++) for (i = 0; i < dow; i++)
......
...@@ -228,7 +228,17 @@ namespace AMDiS { ...@@ -228,7 +228,17 @@ namespace AMDiS {
* Are used to store temporary local dofs of an element. * Are used to store temporary local dofs of an element.
* Thread safe. * Thread safe.
*/ */
std::vector<DegreeOfFreedom* > localIndices; std::vector<DegreeOfFreedom*> localIndices;
/** \brief
* Temporarly used in \ref getGrdAtQPs. Thread safe.
*/
std::vector<DimVec<double>*> grdPhis;
/** \brief
* Temporarly used in \ref getD2AtQPs. Thread safe.
*/
std::vector<DimMat<double>*> D2Phis;
}; };
......
...@@ -25,11 +25,16 @@ namespace AMDiS { ...@@ -25,11 +25,16 @@ namespace AMDiS {
boundaryManager(NULL) boundaryManager(NULL)
{ {
nBasFcts = feSpace->getBasisFcts()->getNumber(); nBasFcts = feSpace->getBasisFcts()->getNumber();
int dim = feSpace->getMesh()->getDim();
localIndices.resize(omp_get_max_threads()); localIndices.resize(omp_get_max_threads());
grdPhis.resize(omp_get_max_threads());
D2Phis.resize(omp_get_max_threads());
for (int i = 0; i < omp_get_max_threads(); i++) { for (int i = 0; i < omp_get_max_threads(); i++) {
localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts); localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts);
grdPhis[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
D2Phis[i] = NEW DimMat<double>(dim, NO_INIT);
} }
} }
...@@ -38,6 +43,8 @@ namespace AMDiS { ...@@ -38,6 +43,8 @@ namespace AMDiS {
{ {
for (int i = 0; i < static_cast<int>(localIndices.size()); i++) { for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
FREE_MEMORY(localIndices[i], DegreeOfFreedom, this->nBasFcts); FREE_MEMORY(localIndices[i], DegreeOfFreedom, this->nBasFcts);
DELETE grdPhis[i];
DELETE D2Phis[i];
} }
} }
......
...@@ -29,9 +29,7 @@ namespace AMDiS { ...@@ -29,9 +29,7 @@ namespace AMDiS {
neighbour_(mesh_->getDim(), NO_INIT), neighbour_(mesh_->getDim(), NO_INIT),
neighbourCoord_(mesh_->getDim(), NO_INIT), neighbourCoord_(mesh_->getDim(), NO_INIT),
oppVertex_(mesh_->getDim(), NO_INIT), oppVertex_(mesh_->getDim(), NO_INIT),
grdLambda_(mesh_->getDim(), NO_INIT) grdLambda_(mesh_->getDim(), NO_INIT)
//parametric_(false)
{ {
projection_.set(NULL); projection_.set(NULL);
...@@ -58,7 +56,7 @@ namespace AMDiS { ...@@ -58,7 +56,7 @@ namespace AMDiS {
{ {
int dim = l.getSize() - 1; int dim = l.getSize() - 1;
static WorldVector<double> world[2]; static WorldVector<double> world[4];
WorldVector<double> *ret = w ? w : &world[omp_get_thread_num()]; WorldVector<double> *ret = w ? w : &world[omp_get_thread_num()];
double c = l[0]; double c = l[0];
...@@ -106,9 +104,9 @@ namespace AMDiS { ...@@ -106,9 +104,9 @@ namespace AMDiS {
for (int i = 0; i < dow; i++) { for (int i = 0; i < dow; i++) {
e1[i] = coords[1][i] - v0[i]; e1[i] = coords[1][i] - v0[i];
if(dim > 1) if (dim > 1)
e2[i] = coords[2][i] - v0[i]; e2[i] = coords[2][i] - v0[i];
if(dim > 2) if (dim > 2)
e3[i] = coords[3][i] - v0[i]; e3[i] = coords[3][i] - v0[i];
} }
...@@ -117,22 +115,22 @@ namespace AMDiS { ...@@ -117,22 +115,22 @@ namespace AMDiS {
det = coords[1][0] - coords[0][0]; det = coords[1][0] - coords[0][0];
break; break;
case 2: case 2:
if(dim == 1) { if (dim == 1) {
det = norm(&e1); det = norm(&e1);
} else { } else {
det = e1[0]*e2[1] - e1[1]*e2[0]; det = e1[0] * e2[1] - e1[1] * e2[0];
} }
break; break;
case 3: case 3:
{ {
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) {
det = norm(&e1); det = norm(&e1);
} else if (dim == 2) { } else if (dim == 2) {
...@@ -141,12 +139,12 @@ namespace AMDiS { ...@@ -141,12 +139,12 @@ namespace AMDiS {
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);
} }
......
...@@ -159,7 +159,7 @@ namespace AMDiS { ...@@ -159,7 +159,7 @@ namespace AMDiS {
nonLinSolver_ = nonLinSolverCreator->create(); nonLinSolver_ = nonLinSolverCreator->create();
nonLinSolver_->setVectorCreator(NEW SystemVector::Creator("temp", nonLinSolver_->setVectorCreator(NEW SystemVector::Creator("temp",
componentSpaces_, componentSpaces_,
numComponents_)); nComponents));
} }
...@@ -171,13 +171,13 @@ namespace AMDiS { ...@@ -171,13 +171,13 @@ namespace AMDiS {
void ProblemNonLinVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag) void ProblemNonLinVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag)
{ {
FUNCNAME("ProblemNonLinVec::buildAfterCoarsen()"); FUNCNAME("ProblemNonLinVec::buildAfterCoarsen()");
int i, numMeshes = static_cast<int>(meshes_.size()); int nMeshes = static_cast<int>(meshes_.size());
for(i = 0; i < numMeshes; i++) { for (int i = 0; i < nMeshes; i++) {
meshes_[i]->dofCompress(); meshes_[i]->dofCompress();
} }
for(i = 0; i < numComponents_; i++) { for (int i = 0; i < nComponents; i++) {
MSG("%d DOFs for %s\n", MSG("%d DOFs for %s\n",
componentSpaces_[i]->getAdmin()->getUsedSize(), componentSpaces_[i]->getAdmin()->getUsedSize(),
componentSpaces_[i]->getName().c_str()); componentSpaces_[i]->getName().c_str());
...@@ -187,15 +187,15 @@ namespace AMDiS { ...@@ -187,15 +187,15 @@ namespace AMDiS {
ElInfo *elInfo; ElInfo *elInfo;
// for all elements ... // for all elements ...
for(i = 0; i < numComponents_; i++) { for (int i = 0; i < nComponents; i++) {
elInfo = stack.traverseFirst(componentMeshes_[i], -1, elInfo = stack.traverseFirst(componentMeshes_[i], -1,
Mesh::CALL_LEAF_EL | Mesh::CALL_LEAF_EL |
Mesh::FILL_BOUND | Mesh::FILL_BOUND |
Mesh::FILL_COORDS | Mesh::FILL_COORDS |
Mesh::FILL_DET | Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA); Mesh::FILL_GRD_LAMBDA);
while(elInfo) { while (elInfo) {
if(solution_->getDOFVector(i)->getBoundaryManager()) { if (solution_->getDOFVector(i)->getBoundaryManager()) {
solution_->getDOFVector(i)-> solution_->getDOFVector(i)->
getBoundaryManager()->fillBoundaryConditions(elInfo, solution_->getDOFVector(i)); getBoundaryManager()->fillBoundaryConditions(elInfo, solution_->getDOFVector(i));
} }
......
...@@ -167,7 +167,7 @@ namespace AMDiS { ...@@ -167,7 +167,7 @@ namespace AMDiS {
nonLinSolver_(NULL), nonLinSolver_(NULL),
updater_(NULL) updater_(NULL)
{ {
u0_.resize(numComponents_); u0_.resize(nComponents);
u0_.set(NULL); u0_.set(NULL);
}; };
...@@ -178,7 +178,7 @@ namespace AMDiS { ...@@ -178,7 +178,7 @@ namespace AMDiS {
int index) int index)
{ {
FUNCNAME("ProblemNonLinVec::setU0()"); FUNCNAME("ProblemNonLinVec::setU0()");
TEST_EXIT(index < numComponents_)("invalid index\n"); TEST_EXIT(index < nComponents)("invalid index\n");
u0_[index] = u0Fct; u0_[index] = u0Fct;
solution_->getDOFVector(index)->interpol(u0Fct); solution_->getDOFVector(index)->interpol(u0Fct);
}; };
......
...@@ -130,7 +130,7 @@ namespace AMDiS { ...@@ -130,7 +130,7 @@ namespace AMDiS {
// before timeout of the runqueue. // before timeout of the runqueue.
int readSerialization = 0; int readSerialization = 0;
::std::string serializationFilename = ""; std::string serializationFilename = "";
GET_PARAMETER(0, "argv->rs", &serializationFilename); GET_PARAMETER(0, "argv->rs", &serializationFilename);
// If the parameter -rs is set, we do nothing here, because the problem will be // If the parameter -rs is set, we do nothing here, because the problem will be
...@@ -153,24 +153,21 @@ namespace AMDiS { ...@@ -153,24 +153,21 @@ namespace AMDiS {
TEST_EXIT(serializationFilename != "")("no serialization file\n"); TEST_EXIT(serializationFilename != "")("no serialization file\n");
MSG("Deserialization from file: %s\n", serializationFilename.c_str()); MSG("Deserialization from file: %s\n", serializationFilename.c_str());
::std::ifstream in(serializationFilename.c_str()); std::ifstream in(serializationFilename.c_str());
deserialize(in); deserialize(in);
in.close(); in.close();
} else { } else {
int globalRefinements = 0;
GET_PARAMETER(0, meshes_[0]->getName() + "->global refinements", "%d",
&globalRefinements);
// Initialize the meshes if there is no serialization file. // Initialize the meshes if there is no serialization file.
for (int i = 0; i < static_cast<int>(meshes_.size()); i++) { for (int i = 0; i < static_cast<int>(meshes_.size()); i++) {
if (initFlag.isSet(INIT_MESH) && if (initFlag.isSet(INIT_MESH) &&
meshes_[i] && meshes_[i] &&
!(meshes_[i]->isInitialized())) { !(meshes_[i]->isInitialized())) {
meshes_[i]->initialize(); meshes_[i]->initialize();
refinementManager_->globalRefine(meshes_[i], globalRefinements);
// do global refinements
int globalRefinements = 0;
GET_PARAMETER(0, meshes_[0]->getName() + "->global refinements", "%d",
&globalRefinements);
for (int gr = 0; gr < static_cast<int>(meshes_.size()); gr++) {
refinementManager_->globalRefine(meshes_[gr], globalRefinements);
}
} }
} }
} }
...@@ -183,23 +180,22 @@ namespace AMDiS { ...@@ -183,23 +180,22 @@ namespace AMDiS {
{ {
FUNCNAME("ProblemVec::createMesh()"); FUNCNAME("ProblemVec::createMesh()");
componentMeshes_.resize(numComponents_); componentMeshes_.resize(nComponents);
::std::map<int, Mesh*> meshForRefinementSet; std::map<int, Mesh*> meshForRefinementSet;
char number[3]; char number[3];
::std::string meshName(""); std::string meshName("");
GET_PARAMETER(0, name_ + "->mesh", &meshName); GET_PARAMETER(0, name_ + "->mesh", &meshName);
TEST_EXIT(meshName != "")("no mesh name spezified\n"); TEST_EXIT(meshName != "")("no mesh name spezified\n");
int dim = 0; int dim = 0;
GET_PARAMETER(0, name_ + "->dim", "%d", &dim); GET_PARAMETER(0, name_ + "->dim", "%d", &dim);
TEST_EXIT(dim)("no problem dimension spezified!\n"); TEST_EXIT(dim)("no problem dimension spezified!\n");
for (int i = 0; i < numComponents_; i++) { for (int i = 0; i < nComponents; i++) {
sprintf(number, "%d", i); sprintf(number, "%d", i);
int refSet = -1; int refSet = -1;
GET_PARAMETER(0, name_ + "->refinement set[" + number + "]", "%d", &refSet); GET_PARAMETER(0, name_ + "->refinement set[" + number + "]", "%d", &refSet);
if (refSet < 0) { if (refSet < 0) {
WARNING("refinement set for component %d not set\n", i);
refSet = 0; refSet = 0;
} }
if (meshForRefinementSet[refSet] == NULL) { if (meshForRefinementSet[refSet] == NULL) {
...@@ -234,30 +230,30 @@ namespace AMDiS { ...@@ -234,30 +230,30 @@ namespace AMDiS {
int degree = 1; int degree = 1;
char number[3]; char number[3];
::std::map< ::std::pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap; std::map< std::pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
int dim = -1; int dim = -1;
GET_PARAMETER(0, name_ + "->dim", "%d", &dim); GET_PARAMETER(0, name_ + "->dim", "%d", &dim);
TEST_EXIT(dim != -1)("no problem dimension spezified!\n"); TEST_EXIT(dim != -1)("no problem dimension spezified!\n");
componentSpaces_.resize(numComponents_, NULL); componentSpaces_.resize(nComponents, NULL);
for (int i = 0; i < numComponents_; i++) { for (int i = 0; i < nComponents; i++) {
sprintf(number, "%d", i); sprintf(number, "%d", i);
GET_PARAMETER(0, name_ + "->polynomial degree[" + number + "]","%d", &degree); GET_PARAMETER(0, name_ + "->polynomial degree[" + number + "]","%d", &degree);
TEST_EXIT(componentSpaces_[i] == NULL)("feSpace already created\n"); TEST_EXIT(componentSpaces_[i] == NULL)("feSpace already created\n");
if (feSpaceMap[::std::pair<Mesh*, int>(componentMeshes_[i], degree)] == NULL) { if (feSpaceMap[std::pair<Mesh*, int>(componentMeshes_[i], degree)] == NULL) {
FiniteElemSpace *newFESpace = FiniteElemSpace *newFESpace =
FiniteElemSpace::provideFESpace(NULL, FiniteElemSpace::provideFESpace(NULL,
Lagrange::getLagrange(dim, degree), Lagrange::getLagrange(dim, degree),
componentMeshes_[i], componentMeshes_[i],
name_ + "->feSpace"); name_ + "->feSpace");
feSpaceMap[::std::pair<Mesh*, int>(componentMeshes_[i], degree)] = newFESpace; feSpaceMap[std::pair<Mesh*, int>(componentMeshes_[i], degree)] = newFESpace;
feSpaces_.push_back(newFESpace); feSpaces_.push_back(newFESpace);
} }
componentSpaces_[i] = componentSpaces_[i] =
feSpaceMap[::std::pair<Mesh*, int>(componentMeshes_[i], degree)]; feSpaceMap[std::pair<Mesh*, int>(componentMeshes_[i], degree)];
} }
// create dof admin for vertex dofs if neccessary // create dof admin for vertex dofs if neccessary
...@@ -278,21 +274,21 @@ namespace AMDiS { ...@@ -278,21 +274,21 @@ namespace AMDiS {
// === create vectors and system matrix === // === create vectors and system matrix ===
systemMatrix_ = NEW Matrix<DOFMatrix*>(numComponents_, numComponents_); systemMatrix_ = NEW Matrix<DOFMatrix*>(nComponents, nComponents);
systemMatrix_->set(NULL); systemMatrix_->set(NULL);
rhs_ = NEW SystemVector("rhs", componentSpaces_, numComponents_); rhs_ = NEW SystemVector("rhs", componentSpaces_, nComponents);
solution_ = NEW SystemVector("solution", componentSpaces_, numComponents_); solution_ = NEW SystemVector("solution", componentSpaces_, nComponents);
char number[10]; char number[10];
::std::string numberedName; std::string numberedName;
for (i = 0; i < numComponents_; i++) { for (i = 0; i < nComponents; i++) {
(*systemMatrix_)[i][i] = NEW DOFMatrix(componentSpaces_[i], (*systemMatrix_)[i][i] = NEW DOFMatrix(componentSpaces_[i],
componentSpaces_[i], "A_ii"); componentSpaces_[i], "A_ii");
(*systemMatrix_)[i][i]->setCoupleMatrix(false); (*systemMatrix_)[i][i]->setCoupleMatrix(false);
sprintf(number, "[%d]", i); sprintf(number, "[%d]", i);
numberedName = "rhs" + ::std::string(number); numberedName = "rhs" + std::string(number);
rhs_->setDOFVector(i, NEW DOFVector<double>(componentSpaces_[i], numberedName)); rhs_->setDOFVector(i, NEW DOFVector<double>(componentSpaces_[i], numberedName));
numberedName = name_ + ::std::string(number); numberedName = name_ + std::string(number);
solution_->setDOFVector(i, NEW DOFVector<double>(componentSpaces_[i], solution_->setDOFVector(i, NEW DOFVector<double>(componentSpaces_[i],
numberedName)); numberedName));
solution_->getDOFVector(i)->refineInterpol(true); solution_->getDOFVector(i)->refineInterpol(true);
...@@ -309,7 +305,7 @@ namespace AMDiS { ...@@ -309,7 +305,7 @@ namespace AMDiS {
FUNCNAME("ProblemVec::createSolver()"); FUNCNAME("ProblemVec::createSolver()");
// === create solver === // === create solver ===
::std::string solverType("no"); std::string solverType("no");
GET_PARAMETER(0, name_ + "->solver", &solverType); GET_PARAMETER(0, name_ + "->solver", &solverType);
OEMSolverCreator<SystemVector> *solverCreator = OEMSolverCreator<SystemVector> *solverCreator =
dynamic_cast<OEMSolverCreator<SystemVector>*>( dynamic_cast<OEMSolverCreator<SystemVector>*>(
...@@ -322,10 +318,10 @@ namespace AMDiS { ...@@ -322,10 +318,10 @@ namespace AMDiS {
solver_->initParameters(); solver_->initParameters();
// === create preconditioners === // === create preconditioners ===
::std::string preconType("no"); std::string preconType("no");
PreconditionerScal *scalPrecon; PreconditionerScal *scalPrecon;
PreconditionerVec *vecPrecon = NEW PreconditionerVec(numComponents_); PreconditionerVec *vecPrecon = NEW PreconditionerVec(nComponents);
GET_PARAMETER(0, name_ + "->solver->left precon", &preconType); GET_PARAMETER(0, name_ + "->solver->left precon", &preconType);
CreatorInterface<PreconditionerScal> *preconCreator = CreatorInterface<PreconditionerScal> *preconCreator =
...@@ -337,12 +333,12 @@ namespace AMDiS { ...@@ -337,12 +333,12 @@ namespace AMDiS {
dynamic_cast<PreconditionerScalCreator*>(preconCreator)-> dynamic_cast<PreconditionerScalCreator*>(preconCreator)->
setName(name_ + "->solver->left precon"); setName(name_ + "->solver->left precon");
for(i = 0; i < numComponents_; i++) { for(i = 0; i < nComponents; i++) {
dynamic_cast<PreconditionerScalCreator*>(preconCreator)-> dynamic_cast<PreconditionerScalCreator*>(preconCreator)->
setSizeAndRow(numComponents_, i); setSizeAndRow(nComponents, i);
scalPrecon = preconCreator->create(); scalPrecon = preconCreator->create();
for(j = 0; j < numComponents_; j++) { for(j = 0; j < nComponents; j++) {
scalPrecon->setMatrix(&(*systemMatrix_)[i][j], j); scalPrecon->setMatrix(&(*systemMatrix_)[i][j], j);
} }
vecPrecon->setScalarPrecon(i, scalPrecon); vecPrecon->setScalarPrecon(i, scalPrecon);
...@@ -351,7 +347,7 @@ namespace AMDiS { ...@@ -351,7 +347,7 @@ namespace AMDiS {
} }
vecPrecon = NEW PreconditionerVec(numComponents_); vecPrecon = NEW PreconditionerVec(nComponents);
GET_PARAMETER(0, name_ + "->solver->right precon", &preconType); GET_PARAMETER(0, name_ + "->solver->right precon", &preconType);
preconCreator = preconCreator =
...@@ -362,12 +358,12 @@ namespace AMDiS { ...@@ -362,12 +358,12 @@ namespace AMDiS {
setName(name_ + "->solver->left precon"); setName(name_ + "->solver->left precon");
for(i = 0; i < numComponents_; i++