Commit 7ae9f1a2 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Code refactoring.

parent efd6a256
...@@ -78,10 +78,9 @@ namespace AMDiS { ...@@ -78,10 +78,9 @@ namespace AMDiS {
FUNCNAME("AdaptInstationary::explicitTimeStrategy()"); FUNCNAME("AdaptInstationary::explicitTimeStrategy()");
// estimate before first adaption // estimate before first adaption
if (adaptInfo->getTime() <= adaptInfo->getStartTime()) { if (adaptInfo->getTime() <= adaptInfo->getStartTime())
problemIteration_->oneIteration(adaptInfo, ESTIMATE); problemIteration_->oneIteration(adaptInfo, ESTIMATE);
}
// increment time // increment time
adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep()); adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
...@@ -136,21 +135,20 @@ namespace AMDiS { ...@@ -136,21 +135,20 @@ namespace AMDiS {
if (problemIteration_->oneIteration(adaptInfo, FULL_ITERATION)) { if (problemIteration_->oneIteration(adaptInfo, FULL_ITERATION)) {
if (!fixedTimestep_ && if (!fixedTimestep_ &&
!adaptInfo->timeToleranceReached() && !adaptInfo->timeToleranceReached() &&
!(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep())) !(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep())) {
{ adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep());
adaptInfo->setTime(adaptInfo->getTime() - adaptInfo->getTimestep()); adaptInfo->setTimestep(adaptInfo->getTimestep() * time_delta_1);
adaptInfo->setTimestep(adaptInfo->getTimestep() * time_delta_1); problemIteration_->endIteration(adaptInfo);
problemIteration_->endIteration(adaptInfo); adaptInfo->incSpaceIteration();
adaptInfo->incSpaceIteration(); break;
break; }
}
} }
adaptInfo->incSpaceIteration(); adaptInfo->incSpaceIteration();
problemIteration_->endIteration(adaptInfo); problemIteration_->endIteration(adaptInfo);
} while(!adaptInfo->spaceToleranceReached() && } while (!adaptInfo->spaceToleranceReached() &&
adaptInfo->getSpaceIteration() <= adaptInfo->getMaxSpaceIteration()); adaptInfo->getSpaceIteration() <= adaptInfo->getMaxSpaceIteration());
} else { } else {
problemIteration_->endIteration(adaptInfo); problemIteration_->endIteration(adaptInfo);
...@@ -162,7 +160,7 @@ namespace AMDiS { ...@@ -162,7 +160,7 @@ namespace AMDiS {
adaptInfo->getTimestepIteration() <= adaptInfo->getMaxTimestepIteration()); adaptInfo->getTimestepIteration() <= adaptInfo->getMaxTimestepIteration());
if (!fixedTimestep_ && adaptInfo->timeErrorLow()) { if (!fixedTimestep_ && adaptInfo->timeErrorLow()) {
adaptInfo->setTimestep(adaptInfo->getTimestep() *time_delta_2); adaptInfo->setTimestep(adaptInfo->getTimestep() * time_delta_2);
if (dbgMode) { if (dbgMode) {
// print information about timestep increase // print information about timestep increase
} }
...@@ -170,12 +168,10 @@ namespace AMDiS { ...@@ -170,12 +168,10 @@ namespace AMDiS {
if (dbgMode) { if (dbgMode) {
std::cout << "=== ADAPT INFO DEBUG MODE ===\n"; std::cout << "=== ADAPT INFO DEBUG MODE ===\n";
std::cout << " Do not increase timestep: \n"; std::cout << " Do not increase timestep: \n";
if (fixedTimestep_) { if (fixedTimestep_)
std::cout << " fixedTimestep = true\n"; std::cout << " fixedTimestep = true\n";
} if (!adaptInfo->timeErrorLow())
if (!adaptInfo->timeErrorLow()) {
adaptInfo->printTimeErrorLowInfo(); adaptInfo->printTimeErrorLowInfo();
}
} }
} }
} }
......
...@@ -84,8 +84,8 @@ namespace AMDiS { ...@@ -84,8 +84,8 @@ namespace AMDiS {
D2Phis.resize(omp_get_overall_max_threads()); D2Phis.resize(omp_get_overall_max_threads());
for (int i = 0; i < omp_get_overall_max_threads(); i++) { for (int i = 0; i < omp_get_overall_max_threads(); i++) {
localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts); localIndices[i] = new DegreeOfFreedom[this->nBasFcts];
localVectors[i] = GET_MEMORY(T, this->nBasFcts); localVectors[i] = new T[this->nBasFcts];
grdPhis[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0); grdPhis[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
grdTmp[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0); grdTmp[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
D2Phis[i] = new DimMat<double>(dim, NO_INIT); D2Phis[i] = new DimMat<double>(dim, NO_INIT);
...@@ -96,8 +96,8 @@ namespace AMDiS { ...@@ -96,8 +96,8 @@ namespace AMDiS {
DOFVectorBase<T>::~DOFVectorBase() DOFVectorBase<T>::~DOFVectorBase()
{ {
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); delete [] localIndices[i];
FREE_MEMORY(localVectors[i], T, this->nBasFcts); delete [] localVectors[i];
delete grdPhis[i]; delete grdPhis[i];
delete grdTmp[i]; delete grdTmp[i];
delete D2Phis[i]; delete D2Phis[i];
......
...@@ -94,9 +94,8 @@ namespace AMDiS { ...@@ -94,9 +94,8 @@ namespace AMDiS {
GeoIndex position = INDEX_OF_DIM(pos, dim); GeoIndex position = INDEX_OF_DIM(pos, dim);
int ndof = 0; int ndof = 0;
for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) { for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position); ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);
}
if (ndof > 0) { if (ndof > 0) {
for (int i = 0; i < mesh->getGeo(position); i++) { for (int i = 0; i < mesh->getGeo(position); i++) {
...@@ -111,7 +110,7 @@ namespace AMDiS { ...@@ -111,7 +110,7 @@ namespace AMDiS {
} }
} }
FREE_MEMORY(dof, DegreeOfFreedom*, mesh->getNumberOfNodes()); delete [] dof;
if (child[0]) if (child[0])
child[0]->deleteElementDOFs(); child[0]->deleteElementDOFs();
......
...@@ -286,19 +286,18 @@ namespace AMDiS { ...@@ -286,19 +286,18 @@ namespace AMDiS {
: rows(r), columns(c) : rows(r), columns(c)
{ {
TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n"); TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
vec = GET_MEMORY(VectorOfFixVecs<FixVecType>*, rows); vec = new VectorOfFixVecs<FixVecType>*[rows];
for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++) { for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++)
*i = new VectorOfFixVecs<FixVecType>(dim, columns, NO_INIT); *i = new VectorOfFixVecs<FixVecType>(dim, columns, NO_INIT);
}
} }
/// destructor /// destructor
virtual ~MatrixOfFixVecs() virtual ~MatrixOfFixVecs()
{ {
for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++) { for (VectorOfFixVecs<FixVecType>** i = &vec[0]; i < &vec[rows]; i++)
delete *i; delete *i;
}
FREE_MEMORY(vec, VectorOfFixVecs<FixVecType>*, rows); delete [] vec;
} }
/// assignment operator /// assignment operator
...@@ -412,52 +411,38 @@ namespace AMDiS { ...@@ -412,52 +411,38 @@ namespace AMDiS {
class WorldVector : public FixVec<T, WORLD> class WorldVector : public FixVec<T, WORLD>
{ {
public: public:
/** \brief /// Calls the corresponding constructor of AlgoVec
* Calls the corresponding constructor of AlgoVec
*/
WorldVector() WorldVector()
: FixVec<T, WORLD>(Global::getGeo(WORLD), NO_INIT) : FixVec<T, WORLD>(Global::getGeo(WORLD), NO_INIT)
{} {}
/** \brief /// Calls the corresponding constructor of AlgoVec
* Calls the corresponding constructor of AlgoVec
*/
WorldVector(InitType initType, T* ini) WorldVector(InitType initType, T* ini)
: FixVec<T, WORLD>(Global::getGeo(WORLD), initType, ini) : FixVec<T, WORLD>(Global::getGeo(WORLD), initType, ini)
{} {}
/** \brief /// Calls the corresponding constructor of AlgoVec
* Calls the corresponding constructor of AlgoVec
*/
WorldVector(InitType initType, const T& ini) WorldVector(InitType initType, const T& ini)
: FixVec<T, WORLD>(Global::getGeo(WORLD), initType, ini) : FixVec<T, WORLD>(Global::getGeo(WORLD), initType, ini)
{} {}
/** \brief /// Sets all entries to d
* Sets all entries to d
*/
inline const WorldVector<T>& operator=(const T& d) inline const WorldVector<T>& operator=(const T& d)
{ {
this->set(d); this->set(d);
return (*this); return (*this);
} }
/** \brief /// Sets the arrays value to the geometric midpoint of the points p1 and p2.
* Sets the arrays value to the geometric midpoint of the points
* p1 and p2;
*/
inline void setMidpoint(const WorldVector<T> &p1, const WorldVector<T> &p2) inline void setMidpoint(const WorldVector<T> &p1, const WorldVector<T> &p2)
{ {
int dow = Global::getGeo(WORLD); int dow = Global::getGeo(WORLD);
for (int i = 0; i < dow; i++) { for (int i = 0; i < dow; i++)
this->valArray[i] = 0.5 * (p1[i] + p2[i]); this->valArray[i] = 0.5 * (p1[i] + p2[i]);
}
} }
/** \brief /// Multplication of a matrix with a vector.
* Multplication of a matrix with a vector.
*/
void multMatrixVec(WorldMatrix<T> &m, WorldVector<T> &v); void multMatrixVec(WorldMatrix<T> &m, WorldVector<T> &v);
}; };
......
...@@ -26,9 +26,9 @@ namespace AMDiS { ...@@ -26,9 +26,9 @@ namespace AMDiS {
("not for DIM != DIM_OF_WORLD\n"); ("not for DIM != DIM_OF_WORLD\n");
WorldVector<double> *basis = new WorldVector<double>[dim]; WorldVector<double> *basis = new WorldVector<double>[dim];
double *lengthBasis = GET_MEMORY(double, dim); double *lengthBasis = new double[dim];
WorldVector<double> *step = new WorldVector<double>[3]; WorldVector<double> *step = new WorldVector<double>[3];
for (int i = 0; i < dim; i++) { for (int i = 0; i < dim; i++) {
TEST_EXIT(numPoints[i] > 0)("numPoints < 1\n"); TEST_EXIT(numPoints[i] > 0)("numPoints < 1\n");
...@@ -121,7 +121,7 @@ namespace AMDiS { ...@@ -121,7 +121,7 @@ namespace AMDiS {
} }
delete [] basis; delete [] basis;
FREE_MEMORY(lengthBasis, double, dim); delete [] lengthBasis;
delete [] step; delete [] step;
delete outFile; delete outFile;
} }
......
...@@ -227,7 +227,7 @@ namespace AMDiS { ...@@ -227,7 +227,7 @@ namespace AMDiS {
TEST_EXIT_DBG((*vertices) == NULL)("vertices != NULL\n"); TEST_EXIT_DBG((*vertices) == NULL)("vertices != NULL\n");
int dimOfPosition = DIM_OF_INDEX(position, dim); int dimOfPosition = DIM_OF_INDEX(position, dim);
*vertices = GET_MEMORY(int, dimOfPosition + 1); *vertices = new int[dimOfPosition + 1];
if ((degree == 4) && (dimOfPosition==1)) { if ((degree == 4) && (dimOfPosition==1)) {
(*vertices)[(nodeIndex != 2) ? 0 : 1] = (*vertices)[(nodeIndex != 2) ? 0 : 1] =
...@@ -667,14 +667,15 @@ namespace AMDiS { ...@@ -667,14 +667,15 @@ namespace AMDiS {
for (int i = 0; i <= dim; i++) { // for all positions for (int i = 0; i <= dim; i++) { // for all positions
int partsAtPos = Global::getGeo(INDEX_OF_DIM(i, dim), dim); int partsAtPos = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
for (int j = 0; j < partsAtPos; j++) { // for all vertices/edges/faces/... for (int j = 0; j < partsAtPos; j++) { // for all vertices/edges/faces/...
int *coordInd = GET_MEMORY(int, i + 1); // indices of relevant coords int *coordInd = new int[i + 1]; // indices of relevant coords
for (int k = 0; k < i + 1; k++) { for (int k = 0; k < i + 1; k++) {
coordInd[k] = Global::getReferenceElement(dim)-> coordInd[k] = Global::getReferenceElement(dim)->
getVertexOfPosition(INDEX_OF_DIM(i,dim), j, k); getVertexOfPosition(INDEX_OF_DIM(i,dim), j, k);
} }
createCoords(coordInd, i + 1, 0, degree); createCoords(coordInd, i + 1, 0, degree);
FREE_MEMORY(coordInd, int, i + 1); delete [] coordInd;
if(static_cast<int>(bary->size()) == nBasFcts) return; if (static_cast<int>(bary->size()) == nBasFcts)
return;
} }
} }
} }
...@@ -683,9 +684,9 @@ namespace AMDiS { ...@@ -683,9 +684,9 @@ namespace AMDiS {
int Lagrange::getNumberOfDOFs(int dim, int degree) int Lagrange::getNumberOfDOFs(int dim, int degree)
{ {
int result = 0; int result = 0;
for (int i = 0; i <= degree; i++) { for (int i = 0; i <= degree; i++)
result += fac(dim - 1 + i) / (fac(i) * fac(dim - 1)); result += fac(dim - 1 + i) / (fac(i) * fac(dim - 1));
}
return result; return result;
} }
...@@ -828,8 +829,9 @@ namespace AMDiS { ...@@ -828,8 +829,9 @@ namespace AMDiS {
if (vec) { if (vec) {
rvec = vec; rvec = vec;
} else { } else {
if(inter) FREE_MEMORY(inter, double, nBasFcts); if (inter)
inter = GET_MEMORY(double, nBasFcts); delete [] inter;
inter = new double[nBasFcts];
inter_size = nBasFcts; inter_size = nBasFcts;
rvec = inter; rvec = inter;
} }
...@@ -928,8 +930,8 @@ namespace AMDiS { ...@@ -928,8 +930,8 @@ namespace AMDiS {
result = indices; result = indices;
} else { } else {
if (localVec) if (localVec)
FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize); delete [] localVec;
localVec = GET_MEMORY(DegreeOfFreedom, nBasFcts); localVec = new DegreeOfFreedom[nBasFcts];
localVecSize = nBasFcts; localVecSize = nBasFcts;
result = localVec; result = localVec;
} }
...@@ -946,9 +948,8 @@ namespace AMDiS { ...@@ -946,9 +948,8 @@ namespace AMDiS {
for (int i = 0; i < num; node0++, i++) { for (int i = 0; i < num; node0++, i++) {
const int *indi = orderOfPositionIndices(el, posIndex, i); const int *indi = orderOfPositionIndices(el, posIndex, i);
for (int k = 0; k < nrDOFs; k++) { for (int k = 0; k < nrDOFs; k++)
result[j++] = dof[node0][n0 + indi[k]]; result[j++] = dof[node0][n0 + indi[k]];
}
} }
} }
} }
...@@ -1005,7 +1006,7 @@ namespace AMDiS { ...@@ -1005,7 +1006,7 @@ namespace AMDiS {
if ((parametric = mesh->getParametric())) { if ((parametric = mesh->getParametric())) {
ERROR_EXIT("not yet implemented\n"); ERROR_EXIT("not yet implemented\n");
} else { } else {
wdetf_qp = GET_MEMORY(double, q->getNumPoints()); wdetf_qp = new double[q->getNumPoints()];
} }
ElInfo *el_info = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); ElInfo *el_info = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
...@@ -1030,7 +1031,7 @@ namespace AMDiS { ...@@ -1030,7 +1031,7 @@ namespace AMDiS {
el_info = stack.traverseNext(el_info); el_info = stack.traverseNext(el_info);
} }
FREE_MEMORY(wdetf_qp, double, q->getNumPoints()); delete [] wdetf_qp;
} }
...@@ -1272,7 +1273,7 @@ namespace AMDiS { ...@@ -1272,7 +1273,7 @@ namespace AMDiS {
Element *el = list->getElement(0); Element *el = list->getElement(0);
const DOFAdmin *admin = drv->getFESpace()->getAdmin(); const DOFAdmin *admin = drv->getFESpace()->getAdmin();
DegreeOfFreedom *pdof = GET_MEMORY(DegreeOfFreedom, basFct->getNumber()); DegreeOfFreedom *pdof = new DegreeOfFreedom[basFct->getNumber()];
basFct->getLocalIndices(el, admin, pdof); basFct->getLocalIndices(el, admin, pdof);
/****************************************************************************/ /****************************************************************************/
...@@ -1314,7 +1315,7 @@ namespace AMDiS { ...@@ -1314,7 +1315,7 @@ namespace AMDiS {
+ 0.75 * (*drv)[pdof[9]]); + 0.75 * (*drv)[pdof[9]]);
if (n <= 1) { if (n <= 1) {
FREE_MEMORY(pdof, DegreeOfFreedom, basFct->getNumber()); delete [] pdof;
return; return;
} }
/****************************************************************************/ /****************************************************************************/
...@@ -1352,7 +1353,7 @@ namespace AMDiS { ...@@ -1352,7 +1353,7 @@ namespace AMDiS {
- 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]]) - 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
+ 0.75 *(*drv)[pdof[9]]); + 0.75 *(*drv)[pdof[9]]);
FREE_MEMORY(pdof, DegreeOfFreedom, basFct->getNumber()); delete [] pdof;
} }
void Lagrange::refineInter3_2d(DOFIndexed<double> *drv, void Lagrange::refineInter3_2d(DOFIndexed<double> *drv,
...@@ -1681,14 +1682,9 @@ namespace AMDiS { ...@@ -1681,14 +1682,9 @@ namespace AMDiS {
if (n < 1) if (n < 1)
return; return;
Element *el; DegreeOfFreedom *pdof = new DegreeOfFreedom[basFct->getNumber()];
DegreeOfFreedom *pdof = GET_MEMORY(DegreeOfFreedom, basFct->getNumber()); const DOFAdmin *admin = drv->getFESpace()->getAdmin();
const DegreeOfFreedom *cdof; Element *el = list->getElement(0);
const DOFAdmin *admin;
el = list->getElement(0);
admin = drv->getFESpace()->getAdmin();
basFct->getLocalIndices(el, admin, pdof); basFct->getLocalIndices(el, admin, pdof);
...@@ -1696,7 +1692,8 @@ namespace AMDiS { ...@@ -1696,7 +1692,8 @@ namespace AMDiS {
/* values on child[0] */ /* values on child[0] */
/****************************************************************************/ /****************************************************************************/
cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL); const DegreeOfFreedom *cdof =
basFct->getLocalIndices(el->getChild(0), admin, NULL);
(*drv)[cdof[2]] = (*drv)[pdof[10]]; (*drv)[cdof[2]] = (*drv)[pdof[10]];
(*drv)[cdof[3]] = (*drv)[cdof[3]] =
...@@ -1763,7 +1760,7 @@ namespace AMDiS { ...@@ -1763,7 +1760,7 @@ namespace AMDiS {
(*drv)[cdof[14]] = (*drv)[pdof[13]]; (*drv)[cdof[14]] = (*drv)[pdof[13]];
if (n <= 1) { if (n <= 1) {
FREE_MEMORY(pdof, DegreeOfFreedom, basFct->getNumber()); delete [] pdof;
return; return;
} }
/****************************************************************************/ /****************************************************************************/
...@@ -1825,7 +1822,7 @@ namespace AMDiS { ...@@ -1825,7 +1822,7 @@ namespace AMDiS {
- 0.03125*(*drv)[pdof[11]] + 0.75*(*drv)[pdof[14]]); - 0.03125*(*drv)[pdof[11]] + 0.75*(*drv)[pdof[14]]);
(*drv)[cdof[14]] = (*drv)[pdof[13]]; (*drv)[cdof[14]] = (*drv)[pdof[13]];
FREE_MEMORY(pdof, DegreeOfFreedom, basFct->getNumber()); delete [] pdof;
} }
void Lagrange::refineInter4_2d(DOFIndexed<double> *drv, void Lagrange::refineInter4_2d(DOFIndexed<double> *drv,
......
...@@ -49,8 +49,8 @@ namespace AMDiS { ...@@ -49,8 +49,8 @@ namespace AMDiS {
int dim = mesh->getDim(); int dim = mesh->getDim();
int el1, el2; int el1, el2;
int *verticesEl1 = GET_MEMORY(int, dim); int *verticesEl1 = new int[dim];
int *verticesEl2 = GET_MEMORY(int, dim); int *verticesEl2 = new int[dim];
int mode = -1; // 0: drop dofs, 1: associate dofs int mode = -1; // 0: drop dofs, 1: associate dofs
int result; int result;
BoundaryType boundaryType; BoundaryType boundaryType;
...@@ -172,8 +172,8 @@ namespace AMDiS { ...@@ -172,8 +172,8 @@ namespace AMDiS {
} }
} }
FREE_MEMORY(verticesEl1, int, dim); delete [] verticesEl1;
FREE_MEMORY(verticesEl2, int, dim); delete [] verticesEl2;
// change periodic vertex dofs // change periodic vertex dofs
for (int i = 0; i < mesh->getNumberOfVertices(); i++) { for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
...@@ -362,25 +362,22 @@ namespace AMDiS { ...@@ -362,25 +362,22 @@ namespace AMDiS {
mesh->addMacroElement(mel[i]); mesh->addMacroElement(mel[i]);
} }
dof = GET_MEMORY(DegreeOfFreedom*, nv); dof = new DegreeOfFreedom*[nv];
coords = GET_MEMORY(WorldVector<double>, nv); coords = new WorldVector<double>[nv];
mel_vertex = GET_MEMORY(int*, ne); mel_vertex = new int*[ne];
for (int i = 0; i < ne; i++) { for (int i = 0; i < ne; i++)
mel_vertex[i] = GET_MEMORY(int, mesh->getGeo(VERTEX)); mel_vertex[i] = new int[mesh->getGeo(VERTEX)];
}
for (int i = 0; i < nv; i++) { for (int i = 0; i < nv; i++)
dof[i] = mesh->getDOF(VERTEX); dof[i] = mesh->getDOF(VERTEX);
}
for (int i = 0; i < ne; i++) { for (int i = 0; i < ne; i++) {
mel[i]->element = mesh->createNewElement(); mel[i]->element = mesh->createNewElement();
(mel)[i]->index = i; (mel)[i]->index = i;
if (dim == 3) {