mfem: Problem with Example 1 for continuous beam-quad.mesh

Thank you for clearing up the Listing 5 confusion and many thanks for the example code! This has been a huge help.

I hope it’s not a bad idea to steal the coordinate data from the topology, e.g. in case it doesn’t follow a grid. In the case of above code this can be done with a minor alteration to the coordinates loop:

   for (int q = 0; q < quadCount+1; q++)
   {
      (*coords)[8*q+0] = mesh->GetVertex(q+0)[0]; 
      (*coords)[8*q+1] = mesh->GetVertex(q+0)[1]; 

      (*coords)[8*q+2] = mesh->GetVertex(q+1)[0]; 
      (*coords)[8*q+3] = mesh->GetVertex(q+1)[1]; 

      (*coords)[8*q+4] = mesh->GetVertex(q+2+quadCount)[0]; 
      (*coords)[8*q+5] = mesh->GetVertex(q+2+quadCount)[1]; 

      (*coords)[8*q+6] = mesh->GetVertex(q+3+quadCount)[0]; 
      (*coords)[8*q+7] = mesh->GetVertex(q+3+quadCount)[1]; 
   }

(Although I will likely avoid using the mesh->GetVertex by keeping the vertices’ coordinates in a separate array.)

Best wishes, Panu

However, you can reproduce the listing 5 mesh with this slightly modified version of your code from above:

#include "mfem.hpp"
#include <fstream>

using namespace std;
using namespace mfem;

int main(int argc, char ** argv)
{

   int quadCount = 7; // Number of quadrilaterals (excluding "fused")

   // Mesh params:      Dim, NVert,        NElem,       NBdrElem, spaceDim
   Mesh* mesh = new Mesh(2, 2*quadCount+4, quadCount+1, 2*quadCount+2, 2);

   // Create vertices
   for (int i = 0; i < 2; i++) // Loop over the 2 rows of vertices
   {
      for (int j = 0; j < quadCount+2; j++) // Loop over the elements
      {
         const double vertex[2] = {1.0*j, 1.0*i};
         mesh->AddVertex(vertex);
      }
   }

   // Create quads by looping over elements
   for (int q = 0; q < quadCount; q++)
   {
      const int quad[4] = {q, q+1, q+quadCount+3, q+quadCount+2};
      mesh->AddQuad(quad);
   }
   // Create the "fused" quad that connects to the initial vertices
   const int quadF[4] = {quadCount, 0, quadCount+2, 2*quadCount+2};
   mesh->AddQuad(quadF);

   // Create boundary at y = 0 by looping over elements
   for (int q = 0; q < quadCount; q++)
   {
      const int bd1[2] = {q, q+1};
      mesh->AddBdrSegment(bd1, 3);
   }
   // Create the "fused" boundaries
   const int bd1F[2] = {quadCount,0};
   mesh->AddBdrSegment(bd1F,3);

   // Create boundary at y = 1 by looping over elements
   for (int q = 0; q < quadCount; q++)
   {
      const int bd2[2] = {q+quadCount+3, q+quadCount+2};
      mesh->AddBdrSegment(bd2, 3);
   }
   // Create the "fused" boundaries
   const int bd2F[2] = {quadCount+2, 2*quadCount+2};
   mesh->AddBdrSegment(bd2F,3);

   // Create coordinate function
   L2_FECollection *fec = new L2_FECollection(1, 2);
   FiniteElementSpace *fes = new FiniteElementSpace(mesh, fec, 2,
                                                    Ordering::byVDIM);
   GridFunction *coords = new GridFunction(fes);
   cout << coords->Size() << endl;
   for (int q = 0; q < quadCount+1; q++)
   {
      // x and y coordinates
      (*coords)[8*q+0] = 1.0*q;     (*coords)[8*q+1] = 0.0;
      (*coords)[8*q+2] = 1.0*(q+1); (*coords)[8*q+3] = 0.0;
      (*coords)[8*q+4] = 1.0*q;     (*coords)[8*q+5] = 1.0;
      (*coords)[8*q+6] = 1.0*(q+1); (*coords)[8*q+7] = 1.0;
   }
   // Set the new coordinates and transfer ownership to the mesh
   mesh->NewNodes(*coords, true);

   // Write the mesh to a file
   {
      ofstream ofs("listing5.mesh");
      mesh->Print(ofs);
      ofs.close();
   }

   delete mesh;
}