diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index 9f99f2c553aa5ee2fca5236d6ec95385709bfd19..db16edd4a19a1f5afc38a03180775cc8dc7c71b3 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -389,16 +389,148 @@ namespace AMDiS { std::vector<DOFAdmin*>::iterator dai = std::find(admin.begin(),admin.end(),localAdmin); - if (dai!= admin.end()) { - ERROR("admin %s is already associated to mesh %s\n", - localAdmin->getName().c_str(), this->getName().c_str()); - } + TEST_EXIT(dai == admin.end()) + ("admin %s is already associated to mesh %s\n", + localAdmin->getName().c_str(), this->getName().c_str()); // ===== adding dofs to already existing elements ============================ - // If adding DOFAdmins to already initilized meshes is required, see older - // AMDiS version (revision < 244) at the same code position. - TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work!\n"); + TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work yet!\n"); + + if (initialized) { + static bool pnd_1d_0[2] = {true, true}; + static bool pnd_1d_1[1] = {false}; + static bool pnd_2d_0[3] = {true, true, true}; + static bool pnd_2d_1[3] = {true, true, false}; + static bool pnd_2d_2[1] = {false}; + static bool pnd_3d_0[4] = {true, true, true, true}; + static bool pnd_3d_1[6] = {false, true, true, true, true, true}; + static bool pnd_3d_2[4] = {true, true, false, false}; + static bool pnd_3d_3[1] = {false}; + static bool *pnd_1d[2] = {pnd_1d_0, pnd_1d_1}; + static bool *pnd_2d[3] = {pnd_2d_0, pnd_2d_1, pnd_2d_2}; + static bool *pnd_3d[4] = {pnd_3d_0, pnd_3d_1, pnd_3d_2, pnd_3d_3}; + static bool **parentNeedsDOF[4] = {NULL, pnd_1d, pnd_2d, pnd_3d}; + + + std::list<struct delmem> delList; + std::map< std::set<DegreeOfFreedom>, DegreeOfFreedom*> dofPtrMap; + const DOFAdmin *vertexAdmin = getVertexAdmin(); + int vertexAdminPreDOFs = vertexAdmin->getNumberOfPreDOFs(VERTEX); + + // finding necessary node number for new admin + + int newNNode = 0; + GeoIndex geoIndex; + + for (int d = 0; d < dim + 1; d++) { + geoIndex = INDEX_OF_DIM(d, dim); + + if (localAdmin->getNumberOfDOFs(geoIndex)>0||nDOF[geoIndex]>0) + newNNode+=getGeo(geoIndex); + } + + bool extendNodes = (newNNode>nNodeEl); + nNodeEl = newNNode; + + TraverseStack stack; + ElInfo *elInfo = NULL; + + elInfo = stack.traverseFirst(this, -1, CALL_EVERY_EL_PREORDER); + while(elInfo) { + Element *element = elInfo->getElement(); + DegreeOfFreedom *newDOF, **oldDOF, **dof = + const_cast<DegreeOfFreedom**>(element->getDOF()); + + int index = 0; + + if (extendNodes) { + oldDOF=dof; + element->setDOFPtrs(); + dof=const_cast<DegreeOfFreedom**>(element->getDOF()); + int index=0,oldIndex=0; + for (int d = 0; d < dim+1; d++) { + geoIndex = INDEX_OF_DIM(d, dim); + if (nDOF[geoIndex]>0) { + for(int i=0;i<getGeo(geoIndex);++i) + dof[index++]=oldDOF[oldIndex++]; + } + else { + if (localAdmin->getNumberOfDOFs(geoIndex)>0) + index+=getGeo(geoIndex); + } + } + + FREE_MEMORY(oldDOF, DegreeOfFreedom*, oldNNodes); + + TEST_EXIT_DBG(index == nNodeEl)("ERROR: Number of entered nodes %f != number of nodes %f\n",index,nNodeEl); + + } + + + index=0; + + // allocate new memory at elements + for(int d = 0; d < dim+1; d++) { + geoIndex = INDEX_OF_DIM(d, dim); + + int numberOfDOFs = localAdmin->getNumberOfDOFs(geoIndex); + int numberOfPreDOFs = nDOF[geoIndex]; + + if (numberOfDOFs>0||numberOfPreDOFs>0) { + + // for all vertices/edges/... + for(int i = 0; i < getGeo(geoIndex); i++, index++) { + std::set<DegreeOfFreedom> dofSet; + for(int j = 0; j < d+1; j++) { + dofSet.insert(dof[element->getVertexOfPosition(geoIndex, i, j)][vertexAdminPreDOFs]); + } + + if(element->isLeaf() || parentNeedsDOF[dim][d][i]) { + if(dofPtrMap[dofSet] == NULL) { + if(localAdmin->getNumberOfDOFs(geoIndex)) { + newDOF = GET_MEMORY(DegreeOfFreedom, numberOfPreDOFs + numberOfDOFs); + // copy old dofs to new memory and free old memory + if(dof[index]) { + for(int j = 0; j < numberOfPreDOFs; j++) { + newDOF[j] = dof[index][j]; + } + // FREE_MEMORY(dof[index], DegreeOfFreedom, numberOfPreDOFs); + // Do not free memory. The information has to be used to identify the part in other elements. + // The memory is only marked for freeing. + struct delmem fm; + fm.ptr=dof[index]; + fm.len=numberOfPreDOFs; + delList.push_back(fm); + } + for(int j = 0; j < numberOfDOFs; j++) { + newDOF[numberOfPreDOFs + j] = localAdmin->getDOFIndex(); + } + dof[index] = newDOF; + } + dofPtrMap[dofSet] = dof[index]; + } else { + dof[index] = dofPtrMap[dofSet]; + } + } + } + } + } + elInfo = stack.traverseNext(elInfo); + } + + // now free the old dof memory: + + std::list<struct delmem>::iterator it=delList.begin(); + + while(it!=delList.end()) { + FREE_MEMORY((*it).ptr, DegreeOfFreedom, (*it).len); + it++; + } + + delList.clear(); + + } admin.push_back(localAdmin); @@ -540,14 +672,14 @@ namespace AMDiS { } - const DOFAdmin *Mesh::createDOFAdmin(const std::string& lname,DimVec<int> lnDOF) + const DOFAdmin *Mesh::createDOFAdmin(const std::string& lname, DimVec<int> lnDOF) { FUNCNAME("Mesh::createDOFAdmin()"); DOFAdmin *localAdmin = new DOFAdmin(this, lname); - for (int i = 0; i < dim+1; i++) - localAdmin->setNumberOfDOFs(i,lnDOF[i]); + for (int i = 0; i < dim + 1; i++) + localAdmin->setNumberOfDOFs(i, lnDOF[i]); addDOFAdmin(localAdmin); diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc index 728525639e26f2e19d8fde8a0abdbf7f92f7baae..4993ea4718262897d247c7f37a9fe8cefafde8a0 100644 --- a/AMDiS/src/ProblemVec.cc +++ b/AMDiS/src/ProblemVec.cc @@ -38,7 +38,7 @@ namespace AMDiS { WARNING("meshes already created\n"); } else { if (initFlag.isSet(CREATE_MESH) || - ((!adoptFlag.isSet(INIT_MESH))&& + (!adoptFlag.isSet(INIT_MESH) && (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) { createMesh(); } @@ -72,7 +72,7 @@ namespace AMDiS { WARNING("feSpaces already created\n"); } else { if (initFlag.isSet(INIT_FE_SPACE) || - (initFlag.isSet(INIT_SYSTEM)&&!adoptFlag.isSet(INIT_FE_SPACE))) { + (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) { createFESpace(); } if (adoptProblem && @@ -87,9 +87,8 @@ namespace AMDiS { TEST_EXIT(feSpaces.size() == 1)("Daran muss ich noch arbeiten!\n"); componentSpaces.resize(nComponents); - for (int i = adoptProblem->getNumComponents(); i < nComponents; i++) { + for (int i = adoptProblem->getNumComponents(); i < nComponents; i++) componentSpaces[i] = componentSpaces[0]; - } } } @@ -125,27 +124,21 @@ namespace AMDiS { WARNING("no solver created\n"); // === create estimator === - if (initFlag.isSet(INIT_ESTIMATOR)) { + if (initFlag.isSet(INIT_ESTIMATOR)) createEstimator(); - } - if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR)) { + + if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR)) estimator = adoptProblem->getEstimator(); - } // === create marker === - if (initFlag.isSet(INIT_MARKER)) { + if (initFlag.isSet(INIT_MARKER)) createMarker(); - } - if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) { + if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) marker = adoptProblem->getMarker(); - } - // === create file writer === - if (initFlag.isSet(INIT_FILEWRITER)) { + if (initFlag.isSet(INIT_FILEWRITER)) createFileWriter(); - } - // === read serialization and init mesh ===