mfem: Implement AMR refiner/derefiner for block systems

Hi, @tangqi , The trouble here is that Mesh::GetNode() expects a “node” index but Element::GetVerticies() returns “vertex” indices. Where “Vertices” are topological entities which tie the mesh together at the corners of elements. On the other hand “Nodes” are geometric entities that describe the shape of an element in physical space.

In your mesh the nodes are described using a discontinuous scalar field which has four entries per element which are arranged in a tensor product ordering. Specifically, the first element contains nodes 0, 1, 2, and 3 where 0 and 1 lie on vertices 0 and 1 in your diagram but node 2 is at vertex 3 and node 3 is at vertex 4. Your code prints out nodes 0 and 1 correctly because these vertices and nodes happen to be numbered the same. However, node 4 is actually the first node of the second element which happens to lie on top of the “1” in your diagram. Similarly, node 3 lies at the location of vertex “4” in your diagram. I believe this explains your results.

The harder question is how do you actually find the node indices that you want? I’m not sure about that yet. I’ll try to look into it and post again unless someone else beats me to it.

Best wishes, Mark

I got AMR working well for the serial version of my code. I will close this one. I may open another one when I encounter more issues in the parallel code. Thanks again for your help!

Oh, we could actually shorten this a little by using Mesh::GetNode() as @jakubcerveny suggested. Rather than asking for vdofs we could just get the dofs which are the nodes indices needed by GetNode(). Something like this…

   double vert[3];

   FiniteElementSpace * fes = mesh.GetNodes()->FESpace();
   
   Array<int> dofs;
   for (int i = 0; i < mesh.GetNE(); i++)
   {
      fes->GetElementDofs(i, dofs);
      for (int j=0; j<dofs.Size(); j++)
      {
        mesh.GetNode(dofs[j], vert);
      }
   }

I think Mesh::GetNode(i, coord) should do what you need.

Thanks for your help, Mark. That was really helpful.

Hi @jakubcerveny I did not find any option to give minimal grid spacing in ThresholdRefiner. For a time dependent problem, it would be nice to have it so that I have a good idea of what kind of time step I should use. Is there a better way to use ThresholdRefiner to achieve this? Note that I just call refiner once per time step.

Qi

Hi, @tangqi ,

I just made a few quick changes to ex6.cpp and ex6p.cpp and confirmed that the size of a BilinearForm object changes after a call to FormLinearSystem() but not after a call to Assemble(). In the parallel version the size of the ParBilinearForm object does not change. This is a confusing inconsistency. I wonder if @v-dobrev or @jakubcerveny would like to offer their opinions on this matter. I suppose it could be argued that BilinearForm is trying to conserve memory by only keeping one version of the matrix but this same argument would apply in the parallel version where two versions of the matrix are stored.

To resolve your current issue I see two options. First you could build two copies of M. This means more memory and computational effort but remember that applying a BC using ess_tdof_list alters select rows of the operator and this would impact the values of z on these rows which may not be desirable. The second approach would be to compute M->Mult(w,z) before the call to FormLinearSystem. This may not be possible but consider it.

Best wishes, Mark

HI, @tangqi ,

Actually DRe->Assemble() should always produce an operator of size fespace.GetVSize(). However, if the mesh happens to be conforming and it is not partitioned across multiple processors then fespace.GetVSize() and fespace.GetTrueVSize() return the same number. If the size of this operator does actually change after calling Assemble() when using non-conforming meshes this would be very confusing. I’m afraid I rarely use the serial versions of these classes so I’m not entirely sure how they function in this situation.

Some of what we’ve talked about is discussed on the page https://mfem.org/pri-dual-vec/.

Best wishes, Mark

Hi, @tangqi ,

OK, I think that makes sense. The calls to MakeTRef only update an optional vector within a GridFunction called t_vec. However, t_vec isn’t actually used for anything so the important part of GridFunction is not being updated.

H1Update->Mult will certainly work for either conforming or non-conforming meshes but its size will be different in these two cases so you need to be careful to use vectors with compatible sizes. Just to be clear, GridFunctions have a length that we call VSize which accounts for all of the degrees of freedom needed to build the finite element approximations inside each element. On non-conforming meshes neighboring elements may not share the same degrees of freedom so the GridFunction needs to be a little longer to store two different sets of degrees of freedom for a given face. However, these two different sets of degrees of freedom are related to each other by a set of constraint equations (to ensure the approximation has the correct continuity across the face). The “true” degree of freedom vector, with size TrueVSize, uses these constraint equations to eliminate dependent degrees of freedom so TrueVSize <= VSize.

Back to your issue. If S will be based on TrueVSize then it isn’t large enough, in general, to store the GridFunction data so you’ll need to store this elsewhere. This could be another, longer, BlockVector or you could just use phi, psi, w, and j. The good news is that I believe we can accomplish this without using S_tmp. Consider the following:

/// Let S contain the true dofs associated with phi, psi, w, and j
/// The GridFunctions will be overwritten with data from S
void AMRUpdate(BlockVector &S,
               Array<int> &true_offset,
               GridFunction &phi,
               GridFunction &psi,
               GridFunction &w,
               GridFunction &j)
{
   // Update the GridFunctions so that they match S
   // (this step can be skipped if you can be sure that the GridFunctions already agree with S)
   phi.SetFromTrueDofs(S.GetBlock(0));
   psi.SetFromTrueDofs(S.GetBlock(1));
   w.SetFromTrueDofs(S.GetBlock(2));
   j.SetFromTrueDofs(S.GetBlock(3));

   FiniteElementSpace* H1FESpace = phi.FESpace();

   //update fem space
   H1FESpace->Update();

   // Compute new dofs on the new mesh
   phi.Update();
   psi.Update();
   w.Update();
   j.Update();

   int fe_size = H1FESpace->GetTrueVSize();

   //update offset vector
   true_offset[0] = 0;
   true_offset[1] = fe_size;
   true_offset[2] = 2*fe_size;
   true_offset[3] = 3*fe_size;
   true_offset[4] = 4*fe_size;

   // Resize S
   S.Update(true_offset);

   // Compute "true" dofs and store them in S
   phi.GetTrueDofs(S.GetBlock(0));
   psi.GetTrueDofs(S.GetBlock(1));
     w.GetTrueDofs(S.GetBlock(2));
     j.GetTrueDofs(S.GetBlock(3));

   H1FESpace->UpdatesFinished();
}

There are many ways to accomplish these tasks and they each have pros and cons.

Best wishes, Mark

BTW, is anisotropic in ZienkiewiczZhuEstimator a flag for anisotropic AMR? I assume I could do isotropic AMR for now.

Yes to both questions.

I agree with @mlstowell that you most likely need to merge the results of the two estimators. If you are not able to use the Refiner interface for that, you can always just create a list of elements and call Mesh::GeneralRefinement, that’s what the Refiner eventually does too.

Derefinement tends to be more tricky than refinement. One needs to derefine carefully so as not to damage the solution (in contrast, excess refinement is usually not a problem). Also, you are only allowed to derefine certain groups of elements (those that were refined previously). The low-level function for derefinement is Mesh::DerefineByError.