load "PETSc" load "SLEPc" macro dimension()2//EOM include "macro_ddm.idp" // macro to get the current mesh size per DOF (used in ReMeshIndicator) // parameters // input: // - Th: the mesh // - Vh: P1 fespace on Th // output: // - h: the Vh finite element finite set to the current mesh size macro MeshSizeComputationDOF(Th, Vh, h) { /*Th mesh * Vh P1 finite element space * h the P1 mesh size value */ /* initialising vector of length of number of vertices (nv) to store number of adjacent edges for each vertex */ Vh NoAdjacentEdges =0 ; /* sum of lengths of edges for edges adjacent to DOF entry_j = sum_(edges adjacent to DOF j) integral_(edge) basis_j */ varf SumMeshSize(unused, v) = intalledges(Th, qfnbpE=1)(v); /* number of edges per DOF (lenEdge = int_e 1) entry_j = sum_(edges adjacent to DOF j) integral_(edge) basis_j / length(edge) */ varf EdgeCount(unused, v) = intalledges(Th, qfnbpE=1)(v/lenEdge); /* number of adjacent edges for each DOF in Vh */ NoAdjacentEdges[] = EdgeCount(0, Vh); /* entry_j = sum_(edges adjacent to DOF j) length(edge) / (number of adjacent edges to DOF j) = mean value of length of edges around DOF j */ h[] = SumMeshSize(0, Vh); h[] = h[]./NoAdjacentEdges[]; } // end of macro MeshSizeComputationDOF // macro to remesh according to the residual indicator // input: // - Th: the mesh // - Ph: P0 fespace on Th // - Vh: fespace on Th // - ErrorIndicator: the varf to evaluate the indicator squared // - alpha: coef to determine DOFs to be refined macro ReMeshIndicator(Th, Ph, Vh, ErrorIndicator, alpha, MaxEta,j) { /* evalutate the mesh size per DOF */ Vh h ; MeshSizeComputationDOF(Th, Vh, h); /* initialising error indicator as a piecewise constant function per triangle (here all 0) and filling it with values from ErrorIndicator squared and taking the square root -> gives us the error on each triangle */ Ph Eta ; Eta[] = ErrorIndicator(0, Ph) ; Eta[] = sqrt(Eta[]) ; /* printing maximal and minimal error over all triangles cout << "Bound for Eta, min = " << Eta[].min << ", max = " << Eta[].max << endl; */ MaxEta[j] = Eta[].max ; /* adapting the mesh to have the required new mesh size */ int NoBigEta = 0 ; /* integer to know how many triangles have big eta*/ int[int] ListOfTriIndices(NoBigEta) ; /* list of indices of the triangles with big eta */ int[int] ListOfGlobalIndices(NoBigEta * Vh.ndofK) ; /* list of global indices of the DOFs in such triangles */ int counter = 0 ; /* parameter for index purposes */ for(int k=0; k alpha*Eta[].max){ NoBigEta++ ; ListOfTriIndices.resize(NoBigEta) ; ListOfGlobalIndices.resize(NoBigEta * Vh.ndofK) ; ListOfTriIndices[counter] = k ; /* global index of the k-th triangle */ for(int l=0; l