[WIP] Stix miniapp in 2D

.gitignore

@@ -294,6 +294,9 @@ miniapps/performance/refined.mesh
 miniapps/performance/mesh.*
 miniapps/performance/sol.*
 
+miniapps/plasma/stix_r2d_dh
+miniapps/plasma/stix_r2d_eb
+
 miniapps/shifted/distance
 miniapps/shifted/ParaViewDistance
 miniapps/shifted/ParaViewLSF

fem/bilininteg.cpp

@@ -4361,7 +4361,7 @@ ScalarVectorProductInterpolator::AssembleElementMatrix2(
       }
    };
 
-   VShapeCoefficient dom_shape_coeff(*Q, dom_fe, Trans.GetSpaceDim());
+   VShapeCoefficient dom_shape_coeff(*Q, dom_fe, dom_fe.GetRangeDim());
 
    elmat.SetSize(ran_fe.GetDof(),dom_fe.GetDof());
 

fem/complex_fem.cpp

@@ -1415,6 +1415,192 @@ ParSesquilinearForm::Update(FiniteElementSpace *nfes)
    if ( pblfi ) { pblfi->Update(nfes); }
 }
 
+bool ParMixedSesquilinearForm::RealInteg()
+{
+   int nint = pblfr->GetTFBFI()->Size() + pblfr->GetDBFI()->Size() +
+              pblfr->GetBBFI()->Size() + pblfr->GetBTFBFI()->Size();
+   return (nint != 0);
+}
+
+bool ParMixedSesquilinearForm::ImagInteg()
+{
+   int nint = pblfi->GetTFBFI()->Size() + pblfi->GetDBFI()->Size() +
+              pblfi->GetBBFI()->Size() + pblfi->GetBTFBFI()->Size();
+   return (nint != 0);
+}
+
+ParMixedSesquilinearForm::ParMixedSesquilinearForm(ParFiniteElementSpace *tr_pf,
+                                                   ParFiniteElementSpace *te_pf,
+                                                   ComplexOperator::Convention
+                                                   convention)
+   : conv(convention),
+     pblfr(new ParMixedBilinearForm(tr_pf, te_pf)),
+     pblfi(new ParMixedBilinearForm(tr_pf, te_pf))
+{}
+
+ParMixedSesquilinearForm::ParMixedSesquilinearForm(ParFiniteElementSpace *tr_pf,
+                                                   ParFiniteElementSpace *te_pf,
+                                                   ParMixedBilinearForm *pbfr,
+                                                   ParMixedBilinearForm *pbfi,
+                                                   ComplexOperator::Convention convention)
+   : conv(convention),
+     pblfr(new ParMixedBilinearForm(tr_pf,te_pf,pbfr)),
+     pblfi(new ParMixedBilinearForm(tr_pf,te_pf,pbfi))
+{}
+
+ParMixedSesquilinearForm::~ParMixedSesquilinearForm()
+{
+   delete pblfr;
+   delete pblfi;
+}
+
+void ParMixedSesquilinearForm::Mult(const ParComplexGridFunction & x,
+                                    ParComplexLinearForm & y) const
+{
+   pblfr->Mult(x.real(), y.real());
+   pblfi->AddMult(x.imag(), y.real(), -1.0);
+
+   pblfr->Mult(x.imag(), y.imag());
+   pblfi->AddMult(x.real(), y.imag(), 1.0);
+
+   if (conv == ComplexOperator::Convention::BLOCK_SYMMETRIC)
+   {
+      y.imag() *= -1.0;
+   }
+
+}
+
+void ParMixedSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator
+                                                   *bfi_real,
+                                                   BilinearFormIntegrator
+                                                   *bfi_imag)
+{
+   if (bfi_real) { pblfr->AddDomainIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi->AddDomainIntegrator(bfi_imag); }
+}
+
+void
+ParMixedSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator
+                                                *bfi_real,
+                                                BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { pblfr->AddBoundaryIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi->AddBoundaryIntegrator(bfi_imag); }
+}
+
+void
+ParMixedSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator
+                                                *bfi_real,
+                                                BilinearFormIntegrator *bfi_imag,
+                                                Array<int> & bdr_marker)
+{
+   if (bfi_real) { pblfr->AddBoundaryIntegrator(bfi_real, bdr_marker); }
+   if (bfi_imag) { pblfi->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+ParMixedSesquilinearForm::AddTraceFaceIntegrator(BilinearFormIntegrator
+                                                 *bfi_real,
+                                                 BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { pblfr->AddTraceFaceIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi->AddTraceFaceIntegrator(bfi_imag); }
+}
+
+void
+ParMixedSesquilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator
+                                                    *bfi_real,
+                                                    BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { pblfr->AddBdrTraceFaceIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi->AddBdrTraceFaceIntegrator(bfi_imag); }
+}
+
+void
+ParMixedSesquilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator
+                                                    *bfi_real,
+                                                    BilinearFormIntegrator *bfi_imag,
+                                                    Array<int> &bdr_marker)
+{
+   if (bfi_real) { pblfr->AddBdrTraceFaceIntegrator(bfi_real, bdr_marker); }
+   if (bfi_imag) { pblfi->AddBdrTraceFaceIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+ParMixedSesquilinearForm::Assemble(int skip_zeros)
+{
+   pblfr->Assemble(skip_zeros);
+   pblfi->Assemble(skip_zeros);
+}
+
+void
+ParMixedSesquilinearForm::Finalize(int skip_zeros)
+{
+   pblfr->Finalize(skip_zeros);
+   pblfi->Finalize(skip_zeros);
+}
+
+ComplexHypreParMatrix *
+ParMixedSesquilinearForm::ParallelAssemble()
+{
+   return new ComplexHypreParMatrix(pblfr->ParallelAssemble(),
+                                    pblfi->ParallelAssemble(),
+                                    true, true, conv);
+}
+
+void ParMixedSesquilinearForm::FormRectangularSystemMatrix(
+   const Array<int> &trial_tdof_list,
+   const Array<int> &test_tdof_list,
+   OperatorHandle &A)
+{
+   OperatorHandle A_r, A_i;
+   if (RealInteg())
+   {
+      pblfr->FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list, A_r);
+   }
+   if (ImagInteg())
+   {
+      pblfi->FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list, A_i);
+   }
+   if (!RealInteg() && !ImagInteg())
+   {
+      MFEM_ABORT("Both Real and Imaginary part of the MixedSesquilinear form are empty");
+   }
+
+   // A = A_r + i A_i
+   A.Clear();
+   if ( A_r.Type() == Operator::Hypre_ParCSR ||
+        A_i.Type() == Operator::Hypre_ParCSR )
+   {
+      ComplexHypreParMatrix * A_hyp =
+         new ComplexHypreParMatrix(A_r.As<HypreParMatrix>(),
+                                   A_i.As<HypreParMatrix>(),
+                                   A_r.OwnsOperator(),
+                                   A_i.OwnsOperator(),
+                                   conv);
+      A.Reset<ComplexHypreParMatrix>(A_hyp, true);
+   }
+   else
+   {
+      ComplexOperator * A_op =
+         new ComplexOperator(A_r.As<Operator>(),
+                             A_i.As<Operator>(),
+                             A_r.OwnsOperator(),
+                             A_i.OwnsOperator(),
+                             conv);
+      A.Reset<ComplexOperator>(A_op, true);
+   }
+   A_r.SetOperatorOwner(false);
+   A_i.SetOperatorOwner(false);
+}
+
+void
+ParMixedSesquilinearForm::Update()
+{
+   if ( pblfr ) { pblfr->Update(); }
+   if ( pblfi ) { pblfi->Update(); }
+}
+
 #endif // MFEM_USE_MPI
 
 }

fem/complex_fem.hpp

@@ -474,7 +474,12 @@ class ParComplexLinearForm : public Vector
    void AddDomainIntegrator(LinearFormIntegrator *lfi_real,
                             LinearFormIntegrator *lfi_imag);
 
-   /// Adds new Domain Integrator, restricted to specific attributes.
+   /** @brief Add new Domain Integrator, restricted to specific attributes.
+
+       Assumes ownership of @a lfi_real and @a lfi_imag.
+
+       The array @a attr_marker is stored internally as a pointer to the
+       given Array<int> object. */
    void AddDomainIntegrator(LinearFormIntegrator *lfi_real,
                             LinearFormIntegrator *lfi_imag,
                             Array<int> &elem_attr_marker);
@@ -682,6 +687,141 @@ class ParSesquilinearForm
    virtual ~ParSesquilinearForm();
 };
 
+/** Class for a parallel mixed sesquilinear form
+
+    A sesquilinear form is a generalization of a mixed bilinear form to
+    complex-valued fields. Sesquilinear forms are linear in the second argument
+    but the first argument involves a complex conjugate in the sense that:
+
+                a(alpha u, beta v) = conj(alpha) beta a(u, v)
+
+    The @a convention argument in the class's constructor is documented in the
+    mfem::ComplexOperator class found in linalg/complex_operator.hpp.
+
+    When supplying integrators to the ParSesquilinearForm either the real or
+    imaginary integrator can be NULL. This indicates that the corresponding
+    portion of the complex-valued material coefficient is equal to zero.
+*/
+class ParMixedSesquilinearForm
+{
+private:
+   ComplexOperator::Convention conv;
+
+   ParMixedBilinearForm *pblfr;
+   ParMixedBilinearForm *pblfi;
+
+   /* These methods check if the real/imag parts of the sesqulinear form are not
+      empty */
+   bool RealInteg();
+   bool ImagInteg();
+
+public:
+   ParMixedSesquilinearForm(ParFiniteElementSpace *tr_pf,
+                            ParFiniteElementSpace *te_pf,
+                            ComplexOperator::Convention
+                            convention = ComplexOperator::HERMITIAN);
+
+   /** @brief Create a ParMixedSesquilinearForm on the ParFiniteElementSpaces
+       @a tr_pf and @a te_pf, using the same integrators as the
+       ParMixedBilinearForms @a pbfr and @a pbfi .
+
+       The pointer @a pf is not owned by the newly constructed object.
+
+       The integrators are copied as pointers and they are not owned by the
+       newly constructed ParSesquilinearForm. */
+   ParMixedSesquilinearForm(ParFiniteElementSpace *tr_pf,
+                            ParFiniteElementSpace *te_pf,
+                            ParMixedBilinearForm *pbfr,
+                            ParMixedBilinearForm *pbfi,
+                            ComplexOperator::Convention
+                            convention = ComplexOperator::HERMITIAN);
+
+   ComplexOperator::Convention GetConvention() const { return conv; }
+   void SetConvention(const ComplexOperator::Convention &
+                      convention) { conv = convention; }
+
+   /// Set the desired assembly level.
+   /** Valid choices are:
+
+       - AssemblyLevel::LEGACYFULL (default)
+       - AssemblyLevel::FULL
+       - AssemblyLevel::PARTIAL
+       - AssemblyLevel::ELEMENT
+       - AssemblyLevel::NONE
+
+       This method must be called before assembly. */
+   void SetAssemblyLevel(AssemblyLevel assembly_level)
+   {
+      pblfr->SetAssemblyLevel(assembly_level);
+      pblfi->SetAssemblyLevel(assembly_level);
+   }
+
+   ParMixedBilinearForm & real() { return *pblfr; }
+   ParMixedBilinearForm & imag() { return *pblfi; }
+   const ParMixedBilinearForm & real() const { return *pblfr; }
+   const ParMixedBilinearForm & imag() const { return *pblfi; }
+
+   /// Matrix multiplication: $ y = M x $
+   void Mult(const ParComplexGridFunction & x,
+             ParComplexLinearForm & y) const;
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                            BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag);
+
+   /** @brief Adds new boundary Integrator, restricted to specific boundary
+       attributes.
+
+       Assumes ownership of @a bfi.
+
+       The array @a bdr_marker is stored internally as a pointer to the given
+       Array<int> object. */
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag,
+                              Array<int> &bdr_marker);
+
+   /// Adds new Face Integrator. Assumes ownership of @a bfi.
+   void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi_real,
+                               BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new boundary Face Integrator. Assumes ownership of @a bfi.
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi_real,
+                                  BilinearFormIntegrator *bfi_imag);
+
+   /** @brief Adds new boundary Face Integrator, restricted to specific boundary
+       attributes.
+
+       Assumes ownership of @a bfi.
+
+       The array @a bdr_marker is stored internally as a pointer to the given
+       Array<int> object. */
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi_real,
+                                  BilinearFormIntegrator *bfi_imag,
+                                  Array<int> &bdr_marker);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   /// Finalizes the matrix initialization.
+   void Finalize(int skip_zeros = 1);
+
+   /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
+   /** The returned matrix has to be deleted by the caller. */
+   ComplexHypreParMatrix *ParallelAssemble();
+
+   void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list,
+                                    const Array<int> &test_tdof_list,
+                                    OperatorHandle &A);
+
+   virtual void Update();
+
+   virtual ~ParMixedSesquilinearForm();
+};
+
 #endif // MFEM_USE_MPI
 
 }