Improved the Jz computation.

cleanup.

fem/lininteg.cpp

@@ -197,6 +197,46 @@ void VectorDomainLFIntegrator::AssembleRHSElementVect(
    }
 }
 
+
+void DomainGradLFIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
+{
+   DenseMatrix Q_ir, dshape, adjJ;
+   Vector vec1, vec2, q_dot_ds;
+
+   int nd  = el.GetDof();
+   int dim = Q.GetVDim();
+   vec1.SetSize(dim);
+   vec2.SetSize(dim);
+   adjJ.SetSize(dim,dim);
+   dshape.SetSize(nd,dim);
+   elvect.SetSize(nd);
+   q_dot_ds.SetSize(nd);
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int intorder = el.GetOrder() + 1;
+      ir = &IntRules.Get(el.GetGeomType(), intorder);
+   }
+
+   Q.Eval(Q_ir, Tr, *ir);
+
+   elvect = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Tr.SetIntPoint(&ip);
+      el.CalcDShape(ip, dshape);
+      CalcAdjugate(Tr.Jacobian(), adjJ);
+      Q_ir.GetColumnReference(i, vec1);
+      vec1 *= ip.weight;
+      adjJ.Mult(vec1, vec2);
+      dshape.Mult(vec2, q_dot_ds);
+      elvect += q_dot_ds;
+   }
+}
+
 void VectorBoundaryLFIntegrator::AssembleRHSElementVect(
    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {

fem/lininteg.hpp

@@ -146,6 +146,28 @@ public:
    using LinearFormIntegrator::AssembleRHSElementVect;
 };
 
+
+/** Class for domain integration of L(v) := (f, grad v), where
+    f=(f1,...,fn) and v scalar. */
+class DomainGradLFIntegrator : public LinearFormIntegrator
+{
+private:
+   VectorCoefficient &Q;
+
+public:
+   /// Constructs a domain integrator with a given VectorCoefficient
+   DomainGradLFIntegrator(VectorCoefficient &QF) : Q(QF) { }
+
+   /** Given a particular Finite Element and a transformation (Tr)
+       computes the element right hand side element vector, elvect. */
+   virtual void AssembleRHSElementVect(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       Vector &elvect);
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+};
+
+
 /** Class for boundary integration of L(v) := (g, v), where
     f=(f1,...,fn) and v=(v1,...,vn). */
 class VectorBoundaryLFIntegrator : public LinearFormIntegrator

miniapps/electromagnetics/maglev/maglev.cpp

@@ -0,0 +1,308 @@
+//                       A 2D static magnetic levitation miniapp
+//
+// Compile with: make maglev
+//
+// Sample runs:  mpirun -np 4 ex1p -m ../data/square-disc.mesh
+
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+
+#include "maglev.hpp"
+#include "../../common/pfem_extras.hpp"
+#include <fstream>
+#include <iostream>
+#include <math.h>
+
+using namespace std;
+using namespace mfem;
+
+void solveForAz(MaglevProblemGeometry &problem, 
+                ParFiniteElementSpace &fes_h1,
+                ParGridFunction &x, 
+                ParGridFunction &conv, 
+                ParGridFunction &mag);
+
+void solveForJz(MaglevProblemGeometry &problem, 
+                ParFiniteElementSpace &fes_h1,
+                ParGridFunction &az,
+                ParGridFunction &jz);
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   int order = 1;
+   int num_refine = 4;
+   double domain_length = 3.0;
+   double domain_height = 2.0;
+   double ha_len = 2.0;
+   double ha_thick = 0.1;
+   int ha_num_magnets = 20;
+   double ha_yoff = 0.05;
+   double ha_vx = 0.0;
+   double ha_mu = 1.32e-6;
+   double ha_sigma = 6.67e5;
+   double ha_remanence = 1.0;
+   double conductor_thick = 0.025;
+   double conductor_yoff = 1.0;
+   double conductor_vx = -100.0;
+   double conductor_mu = 1.257e-6;
+   double conductor_sigma = 3.5e7;
+   double air_mu = 1.257e-6;
+   double air_sigma = 3e-15;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&num_refine, "-nr", "--num_refine",
+                  "Number of times to refine the mesh after initial setup.");
+   args.AddOption(&domain_length, "-dl", "--domain_length",
+                  "Length (x) of the simulation domain (m).");
+   args.AddOption(&domain_height, "-dh", "--domain_height",
+                  "Height (y) of the simulation domain (m).");
+   args.AddOption(&ha_len, "-hl", "--ha_len",
+                  "Physical length of the Halbach array (m).");
+   args.AddOption(&ha_thick, "-ht", "--ha_thick",
+                  "Physical thickness of the Halbach array (m).");
+   args.AddOption(&ha_num_magnets, "-hn", "--ha_num_magnets",
+                  "Number of magnets in the Halbach array.");
+   args.AddOption(&ha_yoff, "-hy", "--ha_yoff",
+                  "Distance between the Halbach array and the conductor (m).");
+   args.AddOption(&ha_vx, "-hv", "--ha_vx",
+                  "Velocity of the Halbach array in the x direction.");
+   args.AddOption(&ha_mu, "-hm", "--ha_mu",
+                  "Permeability of the Halbach array (H/m).");
+   args.AddOption(&ha_sigma, "-hs", "--ha_sigma",
+                  "Conductivity of the halbach array (S/m).");
+   args.AddOption(&ha_remanence, "-hr", "--ha_remanence",
+                  "Remanence of the halbach array (T).");
+   args.AddOption(&conductor_thick, "-ct", "--conductor_thick",
+                  "Physical thickness of the conductor (m).");
+   args.AddOption(&conductor_yoff, "-cy", "--conductor_yoff",
+                  "Distance from the bottom of the domain to the top of the conductor (m).");
+   args.AddOption(&conductor_vx, "-cv", "--conductor_vx",
+                  "Velocity of the conductor in the x direction (m/s).");
+   args.AddOption(&conductor_mu, "-cm", "--conductor_mu",
+                  "Permeability of the conductor (H/m).");
+   args.AddOption(&conductor_sigma, "-cs", "--conductor_sigma",
+                  "Conductivity of the conductor (S/m).");
+   args.AddOption(&air_mu, "-am", "--air_mu",
+                  "Permeability of the air (H/m).");
+   args.AddOption(&air_sigma, "-as", "--air_sigma",
+                  "Conductivity of the air (S/m).");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // Generate the mesh from the given inputs and set up the
+   // the object that holds on to all the problem definitions.
+   double dx = ha_len / double(ha_num_magnets);
+   double dy = min(min(ha_thick, conductor_thick), ha_yoff);
+   int nx = round(domain_length / dx);
+   int ny = round(domain_height / dy);
+   Mesh *mesh = new Mesh(nx, ny, Element::QUADRILATERAL, 1,
+                         domain_length, domain_height);
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int l = 0; l < num_refine; l++)
+   {
+      pmesh->UniformRefinement();
+   }
+   MaglevProblemGeometry problem(conductor_yoff - conductor_thick, conductor_yoff, conductor_vx,
+                                 conductor_mu, conductor_sigma,
+                                 domain_length/2.0 - ha_len/2.0, domain_length/2.0 + ha_len/2.0,
+                                 conductor_yoff+ha_yoff, conductor_yoff+ha_yoff+ha_thick,
+                                 ha_num_magnets, ha_vx, ha_mu, ha_sigma, ha_remanence/ha_mu,
+                                 air_mu, air_sigma);
+
+   // Define a parallel finite element space on the parallel mesh. Here we
+   // use continuous Lagrange finite elements of the specified order.
+   FiniteElementCollection *fec_h1 = new H1_FECollection(order, 2);
+   FiniteElementCollection *fec_rt = new RT_FECollection(order, 2);   
+   FiniteElementCollection *fec_l2 = new L2_FECollection(0, 2);
+   ParFiniteElementSpace *fes_h1 = new ParFiniteElementSpace(pmesh, fec_h1);
+   ParFiniteElementSpace *fes_rt = new ParFiniteElementSpace(pmesh, fec_rt);
+   ParFiniteElementSpace *fes_l2 = new ParFiniteElementSpace(pmesh, fec_l2);
+   ParFiniteElementSpace *fes_vl2 = new ParFiniteElementSpace(pmesh, fec_l2, 2);
+   HYPRE_Int ndof = fes_h1->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << ndof << endl;
+   }
+
+   ParGridFunction az(fes_h1);
+   ParGridFunction conv(fes_vl2);
+   ParGridFunction mag(fes_vl2);
+   solveForAz(problem, *fes_h1, az, conv, mag);
+
+   //Compute Bfield = curl A so we can visualize it
+   ParGridFunction bfield(fes_rt);
+   miniapps::ParDiscreteCurlOperator curl(fes_h1,fes_rt);
+   curl.Assemble();
+   curl.Finalize();
+   curl.Mult(az, bfield);
+
+   ParGridFunction jz(fes_h1);
+   solveForJz(problem, *fes_h1, az, jz);
+
+   //Make grid functions out of these coefficients for visualization
+   ParGridFunction vx(fes_l2);
+   ParGridFunction sigma(fes_l2);
+   ParGridFunction mu(fes_l2);
+   VxCoeff vx_coeff(problem);
+   SigmaCoeff sigma_coeff(problem);
+   MuCoeff mu_coeff(problem);
+   vx.ProjectCoefficient(vx_coeff);
+   sigma.ProjectCoefficient(sigma_coeff);
+   mu.ProjectCoefficient(mu_coeff);
+
+   VisItDataCollection visit_dc("maglev", pmesh);
+   visit_dc.RegisterField("Az", &az);
+   visit_dc.RegisterField("Jz", &jz);
+   visit_dc.RegisterField("B", &bfield);
+   visit_dc.RegisterField("convection_coeff", &conv);
+   visit_dc.RegisterField("mag_coeff", &mag);
+   visit_dc.RegisterField("vx", &vx);
+   visit_dc.RegisterField("sigma", &sigma);
+   visit_dc.RegisterField("mu", &mu);
+   visit_dc.Save();
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+
+//Solve Delta Az - mu sigma v dot grad Az = curl M
+//x:     The solution Az
+//conv:  The convection coefficient (for vizualization)
+//mag:   The magnitization coefficient (for vizualization)
+void solveForAz(MaglevProblemGeometry &problem, 
+                ParFiniteElementSpace &fes_h1,
+                ParGridFunction &az,
+                ParGridFunction &conv,
+                ParGridFunction &mag)
+{
+   // Set up the parallel bilinear form dc(.,.) on the finite element space
+   // corresponding to the Laplacian operator Delta Az - mu sigma v dot grad Az.
+   ConstantCoefficient one(-1.0);
+   ConvectionCoeff conv_coeff(problem);
+   ParBilinearForm dc(&fes_h1);
+   dc.AddDomainIntegrator(new DiffusionIntegrator(one));
+   dc.AddDomainIntegrator(new ConvectionIntegrator(conv_coeff));
+   dc.Assemble();
+
+   // Determine the list of true (i.e. parallel conforming) essential
+   // boundary dofs. In this example, the boundary conditions are defined
+   // by marking all the boundary attributes from the mesh as essential
+   // (Dirichlet) and converting them to a list of true dofs.
+   ParMesh *pmesh = fes_h1.GetParMesh();
+   Array<int> ess_tdof_list(fes_h1.GlobalTrueVSize());
+   ess_tdof_list = 0;
+   if (pmesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fes_h1.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // Set up the parallel linear form b(.) which corresponds to the
+   // right-hand side of the FEM linear system, which in this case is
+   // (1,phi_i) where phi_i are the basis functions in fes_h1.
+   ParLinearForm b(&fes_h1);
+   MagnetizationCoeff mag_coeff(problem);
+   b.AddDomainIntegrator(new DomainGradLFIntegrator(mag_coeff));
+   b.Assemble();
+
+   // Assemble the parallel bilinear form and the corresponding linear
+   // system, applying any necessary transformations such as
+   HypreParMatrix DC;
+   Vector B, AZ;
+   dc.FormLinearSystem(ess_tdof_list, az, b, DC, AZ, B);
+
+   // Define and apply a parallel GMRES solver for AX=B with the BoomerAMG
+   // preconditioner from hypre.
+   HypreParaSails sails(DC);
+   HypreGMRES gmres(DC);
+   gmres.SetTol(1e-7);
+   gmres.SetKDim(250);
+   gmres.SetMaxIter(10000);
+   gmres.SetPrintLevel(2);
+   gmres.SetPreconditioner(sails);
+   gmres.Mult(B, AZ);
+
+   // Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   dc.RecoverFEMSolution(AZ, b, az);
+
+   //Project these coefficients onto corresponding grid functions for vizualization
+   conv.ProjectCoefficient(conv_coeff);
+   mag.ProjectCoefficient(mag_coeff);
+}
+
+
+//Solve Jz = -Delta 1/mu Az
+void solveForJz(MaglevProblemGeometry &problem, 
+                ParFiniteElementSpace &fes_h1,
+                ParGridFunction &az,
+                ParGridFunction &jz)
+{
+   Array<int> ess_tdof_list(fes_h1.GlobalTrueVSize());
+   ess_tdof_list = 0;
+
+   //Set up M matrix for the left hand side
+   ConstantCoefficient one(-1.0);
+   ParBilinearForm m(&fes_h1);
+   m.AddDomainIntegrator(new MassIntegrator(one));
+   m.Assemble();
+
+   //Set up the biliinear for the K matrix for -Delta 1/mu on the right hand side
+   MuInvCoeff muinv_coeff(problem);
+   ParBilinearForm k(&fes_h1);
+   k.AddDomainIntegrator(new DiffusionIntegrator(muinv_coeff));
+   k.Assemble();
+
+   //Compute B = K AZ
+   ParGridFunction b(&fes_h1);
+   HypreParMatrix K;
+   HypreParVector *AZ = az.GetTrueDofs();
+   HypreParVector B(*AZ);
+   k.FormSystemMatrix(ess_tdof_list, K);
+   K.Mult(*AZ, B);
+
+   //Now Set up the linear system for M JZ = B and get jz
+   jz = 0.0;
+   HypreParVector *JZ = jz.GetTrueDofs();
+   HypreParMatrix M;
+   m.FormSystemMatrix(ess_tdof_list, M);
+   HypreBoomerAMG amg(M);
+   HypreGMRES gmres(M);
+   gmres.SetTol(1e-12);
+   gmres.SetKDim(250);
+   gmres.SetMaxIter(10000);
+   gmres.SetPrintLevel(2);
+   gmres.SetPreconditioner(amg);
+   gmres.Mult(B, *JZ);
+   jz = *JZ;
+
+   delete AZ;
+}

miniapps/electromagnetics/maglev/maglev.hpp

@@ -0,0 +1,318 @@
+#include "mfem.hpp"
+
+class MaglevProblemGeometry
+{
+   public:
+   MaglevProblemGeometry(double conductor_bottom, double conductor_top,
+                         double conductor_vx, double conductor_mu, double conductor_sigma,
+                         double ha_left, double ha_right, double ha_bottom, double ha_top,
+                         int ha_num_magnets, double ha_vx, double ha_mu, double ha_sigma, double ha_mag,
+                         double air_mu,  double air_sigma) :
+      m_conductor_bottom(conductor_bottom),
+      m_conductor_top(conductor_top),
+      m_conductor_vx(conductor_vx),
+      m_conductor_mu(conductor_mu),
+      m_conductor_sigma(conductor_sigma),
+      m_ha_left(ha_left),
+      m_ha_right(ha_right),
+      m_ha_bottom(ha_bottom),
+      m_ha_top(ha_top),
+      m_ha_num_magnets(ha_num_magnets),
+      m_ha_vx(ha_vx),
+      m_ha_mu(ha_mu),
+      m_ha_sigma(ha_sigma),
+      m_ha_mag(ha_mag),
+      m_air_vx(0.0),
+      m_air_mu(air_mu),
+      m_air_sigma(air_sigma) {m_magnet_size = (m_ha_right - m_ha_left) / double(m_ha_num_magnets);}
+
+   inline int getMagnetNumber(const mfem::Vector &x);
+   inline double getMu(const mfem::Vector &x);
+   inline double getSigma(const mfem::Vector &x);
+   inline double getVx(const mfem::Vector &x);
+   inline void getMuSigmaV(const mfem::Vector &x, mfem::Vector &out);
+   inline void getMPerp(const mfem::Vector &x, mfem::Vector &out);
+
+   private:
+   double m_conductor_bottom;
+   double m_conductor_top;
+   double m_conductor_vx;
+   double m_conductor_mu;   
+   double m_conductor_sigma;
+   double m_ha_left;
+   double m_ha_right;
+   double m_ha_bottom;
+   double m_ha_top;
+   int m_ha_num_magnets;
+   double m_ha_vx;
+   double m_ha_mu;
+   double m_ha_sigma;
+   double m_ha_mag;
+   double m_air_vx;
+   double m_air_mu;
+   double m_air_sigma;
+   double m_magnet_size;
+};
+
+
+class ConvectionCoeff : public mfem::VectorCoefficient
+{
+   public:
+   ConvectionCoeff(MaglevProblemGeometry &problem) :
+      VectorCoefficient(2),
+      m_problem(problem) {}
+
+   virtual void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
+                     const mfem::IntegrationPoint &ip);
+
+   private:
+   MaglevProblemGeometry &m_problem;
+};
+
+
+class MagnetizationCoeff : public mfem::VectorCoefficient
+{
+   public:
+   MagnetizationCoeff(MaglevProblemGeometry problem) :
+      VectorCoefficient(2),
+      m_problem(problem) {}
+
+   virtual void Eval(mfem::Vector &V, mfem::ElementTransformation &T,
+                     const mfem::IntegrationPoint &ip);
+
+   private:
+   MaglevProblemGeometry &m_problem;
+};
+
+
+class MuInvCoeff : public mfem::Coefficient
+{
+   public:
+   MuInvCoeff(MaglevProblemGeometry &problem) :
+      Coefficient(),
+      m_problem(problem) {}
+
+   virtual double Eval(mfem::ElementTransformation &T,
+                       const mfem::IntegrationPoint &ip);
+
+   private:
+   MaglevProblemGeometry &m_problem;
+};
+
+
+class VxCoeff : public mfem::Coefficient
+{
+   public:
+   VxCoeff(MaglevProblemGeometry &problem) :
+      Coefficient(),
+      m_problem(problem) {}
+
+   virtual double Eval(mfem::ElementTransformation &T,
+                       const mfem::IntegrationPoint &ip);
+
+   private:
+   MaglevProblemGeometry &m_problem;
+};
+
+
+class MuCoeff : public mfem::Coefficient
+{
+   public:
+   MuCoeff(MaglevProblemGeometry &problem) :
+      Coefficient(),
+      m_problem(problem) {}
+
+   virtual double Eval(mfem::ElementTransformation &T,
+                       const mfem::IntegrationPoint &ip);
+
+   private:
+   MaglevProblemGeometry &m_problem;
+};
+
+
+class SigmaCoeff : public mfem::Coefficient
+{
+   public:
+   SigmaCoeff(MaglevProblemGeometry &problem) :
+      Coefficient(),
+      m_problem(problem) {}
+
+   virtual double Eval(mfem::ElementTransformation &T,
+                       const mfem::IntegrationPoint &ip);
+
+   private:
+   MaglevProblemGeometry &m_problem;
+};
+
+
+inline int MaglevProblemGeometry::getMagnetNumber(const mfem::Vector &x)
+{
+   return int((x[0] - m_ha_left) / m_magnet_size);
+}
+
+
+inline double MaglevProblemGeometry::getVx(const mfem::Vector &x)
+{
+   double vx;
+   if (x[1] >= m_conductor_bottom && x[1] <= m_conductor_top)
+   {
+      vx = m_conductor_vx;
+   }
+   else if (x[0] >= m_ha_left && x[0] <= m_ha_right && x[1] >= m_ha_bottom && x[1] <= m_ha_top)
+   {
+      vx = m_ha_vx;
+   }
+   else
+   {
+      vx = m_air_vx;
+   }
+   return vx;
+}
+
+
+inline double MaglevProblemGeometry::getMu(const mfem::Vector &x)
+{
+   double mu;
+   if (x[1] >= m_conductor_bottom && x[1] <= m_conductor_top)
+   {
+      mu = m_conductor_mu;
+   }
+   else if (x[0] >= m_ha_left && x[0] <= m_ha_right && x[1] >= m_ha_bottom && x[1] <= m_ha_top)
+   {
+      mu = m_ha_mu;
+   }
+   else
+   {
+      mu = m_air_mu;
+   }
+   return mu;
+}
+
+
+inline double MaglevProblemGeometry::getSigma(const mfem::Vector &x)
+{
+   double sigma;
+   if (x[1] >= m_conductor_bottom && x[1] <= m_conductor_top)
+   {
+      sigma = m_conductor_sigma;
+   }
+   else if (x[0] >= m_ha_left && x[0] <= m_ha_right && x[1] >= m_ha_bottom && x[1] <= m_ha_top)
+   {
+      sigma = m_ha_sigma;
+   }
+   else
+   {
+      sigma = m_air_sigma;
+   }
+   return sigma;
+}
+
+
+inline void MaglevProblemGeometry::getMuSigmaV(const mfem::Vector &x, mfem::Vector &out)
+{
+   out = 0.0;
+   out[0] = getMu(x)*getSigma(x)*getVx(x);
+}
+
+
+inline void MaglevProblemGeometry::getMPerp(const mfem::Vector &x, mfem::Vector &out)
+{
+   out = 0.0;
+   if (x[0] >= m_ha_left && x[0] <= m_ha_right && x[1] >= m_ha_bottom && x[1] <= m_ha_top)
+   {
+      int mag = getMagnetNumber(x) % 4;
+      double a = 0.0;
+      double b = 0.0;
+      if (mag == 0)
+      {
+         a = -m_ha_mag;
+      }
+      else if (mag == 1)
+      {
+         b = m_ha_mag;
+      }
+      else if (mag == 2)
+      {
+         a = m_ha_mag;
+      }
+      else
+      {
+         b = -m_ha_mag;
+      }
+
+      //perpendicular to (a,b)
+      out[0] = -b;
+      out[1] = a;
+   }
+}
+
+
+void ConvectionCoeff::Eval(mfem::Vector &V, mfem::ElementTransformation &T,
+                           const mfem::IntegrationPoint &ip)
+{
+   double x[3];
+   mfem::Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+
+   V.SetSize(2);
+   m_problem.getMuSigmaV(transip, V);
+}
+
+
+void MagnetizationCoeff::Eval(mfem::Vector &V, mfem::ElementTransformation &T,
+                              const mfem::IntegrationPoint &ip)
+{
+   double x[3];
+   mfem::Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+
+   V.SetSize(2);
+   m_problem.getMPerp(transip, V);
+}
+
+
+//1/mu for the Jz = -laplace(1/mu Az) calculation
+double MuInvCoeff::Eval(mfem::ElementTransformation &T,
+                               const mfem::IntegrationPoint &ip)
+{
+   double x[3];
+   mfem::Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+   return 1.0 / m_problem.getMu(transip);
+}
+
+
+double VxCoeff::Eval(mfem::ElementTransformation &T,
+                     const mfem::IntegrationPoint &ip)
+{
+   double x[3];
+   mfem::Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+   return m_problem.getVx(transip);
+}
+
+
+double MuCoeff::Eval(mfem::ElementTransformation &T,
+                     const mfem::IntegrationPoint &ip)
+{
+   double x[3];
+   mfem::Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+   return m_problem.getMu(transip);
+}
+
+
+double SigmaCoeff::Eval(mfem::ElementTransformation &T,
+                        const mfem::IntegrationPoint &ip)
+{
+   double x[3];
+   mfem::Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+   return 1.0 / m_problem.getSigma(transip);
+}

miniapps/electromagnetics/maglev/makefile

@@ -0,0 +1,66 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../../..
+MFEM_BUILD_DIR ?= ../../..
+SRC = $(MFEM_DIR)/miniapps/electromagnetics/maglev/
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+TEST_MK = $(MFEM_DIR)/config/test.mk
+COMMON_O=../../common/pfem_extras.o
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_MINIAPPS = maglev
+PAR_MINIAPPS = 
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+.PRECIOUS: %.o
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(MINIAPPS)
+
+# Rules for building the miniapps
+$(MINIAPPS): \
+%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ ../../common/pfem_extras.o $(MFEM_LIBS)
+
+$(COMMON_O): \
+%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
+
+MFEM_TESTS = MINIAPPS
+include $(TEST_MK)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	rm -rf maglev
+	rm -rf *.000000 *.00000 maglev_000000 maglev_000000.mfem_root