minor

miniapps/dpg/pacoustics.cpp

@@ -180,6 +180,7 @@ int main(int argc, char *argv[])
 {
    Mpi::Init();
    int myid = Mpi::WorldRank();
+   int num_procs = Mpi::WorldSize();
    Hypre::Init();
 
    const char *mesh_file = "../../data/inline-quad.mesh";
@@ -261,6 +262,9 @@ int main(int argc, char *argv[])
    dim = mesh.Dimension();
    MFEM_VERIFY(dim > 1, "Dimension = 1 is not supported in this example");
 
+   // mesh.EnsureNCMesh(true);
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+
    CartesianPML * pml = nullptr;
    if (with_pml)
    {
@@ -269,8 +273,7 @@ int main(int argc, char *argv[])
       pml->SetOmega(omega);
    }
 
-   mesh.EnsureNCMesh(true);
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+
    mesh.Clear();
 
    Array<int> attr;
@@ -549,8 +552,13 @@ int main(int argc, char *argv[])
    int dof0 = 0;
 
    ParGridFunction p_r, p_i, u_r, u_i;
+   ParGridFunction p_sub_r, p_sub_i, p_exact_gf_r, p_exact_gf_i, p_th;
 
    ParaViewDataCollection * paraview_dc = nullptr;
+   ParaViewDataCollection * paraview_pml = nullptr;
+   ParaViewDataCollection * paraview_th = nullptr;
+   Array<int> dattr; dattr.Append(1);
+   ParSubMesh psubmesh = ParSubMesh::CreateFromDomain(pmesh,dattr);
 
    if (paraview)
    {
@@ -565,6 +573,20 @@ int main(int argc, char *argv[])
       paraview_dc->RegisterField("p_i",&p_i);
       paraview_dc->RegisterField("u_r",&u_r);
       paraview_dc->RegisterField("u_i",&u_i);
+      if (pml)
+      {
+         // paraview_pml = new ParaViewDataCollection(enum_str[prob], &psubmesh);
+         // paraview_pml->SetPrefixPath("ParaView/Acoustics/PML");
+         // paraview_pml->SetLevelsOfDetail(order);
+         // paraview_pml->SetCycle(0);
+         // paraview_pml->SetDataFormat(VTKFormat::BINARY);
+         // paraview_pml->SetHighOrderOutput(true);
+         // paraview_pml->SetTime(0.0); // set the time
+         // paraview_pml->RegisterField("p_sub_r",&p_sub_r);
+         // paraview_pml->RegisterField("p_sub_i",&p_sub_i);
+         // paraview_pml->RegisterField("p_exact_r",&p_exact_gf_r);
+         // paraview_pml->RegisterField("p_exact_i",&p_exact_gf_i);
+      }
    }
 
    if (static_cond) { a->EnableStaticCondensation(); }
@@ -718,6 +740,53 @@ int main(int argc, char *argv[])
 
       p_r.MakeRef(p_fes, x, 0);
       p_i.MakeRef(p_fes, x, offsets.Last());
+      if (it == pr)
+      {
+         ParSubMesh psubmesh = ParSubMesh::CreateFromDomain(pmesh,dattr);
+         ParFiniteElementSpace p_subfes(&psubmesh,p_fec);
+         p_sub_r.SetSpace(&p_subfes);
+         p_sub_i.SetSpace(&p_subfes);
+         psubmesh.Transfer(p_r,p_sub_r);
+         psubmesh.Transfer(p_i,p_sub_i);
+
+         // paraview_pml->SetMesh(&psubmesh);
+         // paraview_pml->SetCycle(it);
+         // paraview_pml->SetTime((double)it);
+         // paraview_pml->Save();
+
+         FunctionCoefficient p_ex_r2(p_exact_r);
+         FunctionCoefficient p_ex_i2(p_exact_i);
+         p_exact_gf_r.SetSpace(&p_subfes);
+         p_exact_gf_r.ProjectCoefficient(p_ex_r2);
+         p_exact_gf_i.SetSpace(&p_subfes);
+         p_exact_gf_i.ProjectCoefficient(p_ex_i2);
+
+
+         p_exact_gf_r-=p_sub_r;
+         p_exact_gf_i-=p_sub_i;
+         p_th.SetSpace(&p_subfes);
+         paraview_th = new ParaViewDataCollection(enum_str[prob], &psubmesh);
+         paraview_th->SetPrefixPath("ParaView/Acoustics/PML/TimeHarmonic");
+         paraview_th->SetLevelsOfDetail(order);
+         paraview_th->SetCycle(0);
+         paraview_th->SetDataFormat(VTKFormat::BINARY);
+         paraview_th->SetHighOrderOutput(true);
+         paraview_th->SetTime(0.0); // set the time
+         paraview_th->RegisterField("p_th",&p_th);
+         int num_frames = 64;
+         for (int j = 0; j<num_frames; j++)
+         {
+            p_th = 0.0;
+            double t = (double)(j % num_frames) / num_frames;
+            add(cos(2.0*M_PI*t), p_exact_gf_r, sin(2.0*M_PI*t), p_exact_gf_i, p_th);
+            // p_th = p_exact_gf_i;
+            paraview_th->SetMesh(&psubmesh);
+            paraview_th->SetCycle(j);
+            paraview_th->SetTime(t);
+            paraview_th->Save();
+         }
+      }
+
 
       u_r.MakeRef(u_fes,x, offsets[1]);
       u_i.MakeRef(u_fes,x, offsets.Last()+offsets[1]);
@@ -821,13 +890,88 @@ int main(int argc, char *argv[])
          pmesh.UniformRefinement();
       }
       if (pml) { pml->SetAttributes(&pmesh); }
-      for (int i =0; i<trial_fes.Size(); i++)
+      for (int i = 0; i<trial_fes.Size(); i++)
       {
          trial_fes[i]->Update(false);
       }
       a->Update();
    }
 
+   {
+      int num_frames = 32;
+
+      if (paraview)
+      {
+         ParGridFunction x_t(p_r);
+         ParGridFunction x_r(p_r);
+         ParGridFunction x_i(p_i);
+         ParaViewDataCollection * paraview_dct = new ParaViewDataCollection(
+            enum_str[prob], &pmesh);
+         ostringstream oss;
+         if (theta > 0.0)
+         {
+            oss << "ParaView/Acoustics/TimeHarmonicAMR";
+         }
+         else
+         {
+            oss << "ParaView/Acoustics/TimeHarmonicUniform";
+         }
+         paraview_dct->SetPrefixPath(oss.str());
+         paraview_dct->SetLevelsOfDetail(order);
+         paraview_dct->SetCycle(0);
+         paraview_dct->SetDataFormat(VTKFormat::BINARY);
+         paraview_dct->SetHighOrderOutput(true);
+         paraview_dct->SetTime(0.0); // set the time
+         paraview_dct->RegisterField("p_t",&x_t);
+
+         for (int j = 0; j<num_frames; j++)
+         {
+            x_t = 0.0;
+            double t = (double)(j % num_frames) / num_frames;
+            add(cos(2.0*M_PI*t), x_r, sin(2.0*M_PI*t), x_i, x_t);
+            paraview_dct->SetCycle(j);
+            paraview_dct->SetTime(t);
+            paraview_dct->Save();
+         }
+         delete paraview_dct;
+      }
+      MPI_Barrier(MPI_COMM_WORLD); // try to prevent streams from mixing
+
+
+      ParGridFunction p_t(p_r);
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      string keys;
+      keys = (dim == 3) ? "keys macF\n" : keys = "keys amrRljcUUuu\n";
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n"
+               << "solution\n" << pmesh << p_t << keys << "autoscale off\n"
+               << "window_title 'Harmonic Solution (t = 0.0 T)'"
+               << "pause\n" << flush;
+
+      if (myid == 0)
+      {
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+
+      int i = 0;
+      while (sol_sock)
+      {
+         double t = (double)(i % num_frames) / num_frames;
+         ostringstream oss;
+         oss << "Harmonic Solution (t = " << t << " T)";
+
+         add(cos(2.0*M_PI*t), p_r, sin(2.0*M_PI*t), p_i, p_t);
+         sol_sock << "parallel " << num_procs << " " << myid << "\n";
+         sol_sock << "solution\n" << pmesh << p_t
+                  << "window_title '" << oss.str() << "'" << flush;
+         i++;
+      }
+   }
+
+
    if (paraview)
    {
       delete paraview_dc;