Merge remote-tracking branch 'origin/master' into block-fem-dev

# Conflicts:
#	fem/linearform.hpp
#	fem/plinearform.hpp

fem/CMakeLists.txt

@@ -12,6 +12,7 @@
 set(SRCS
   bilinearform.cpp
   bilininteg.cpp
+  block_fem.cpp
   coefficient.cpp
   datacollection.cpp
   eltrans.cpp
@@ -34,6 +35,7 @@ set(SRCS
 set(HDRS
   bilinearform.hpp
   bilininteg.hpp
+  block_fem.hpp
   coefficient.hpp
   datacollection.hpp
   eltrans.hpp

fem/block_fem.cpp

@@ -0,0 +1,908 @@
+// 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.
+
+#include "block_fem.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+BlockGridFunction::BlockGridFunction(BlockFiniteElementSpace *f)
+   : BlockObject(f->GetNBlocks()), Vector(f->GetVSize()), fes(f)
+{
+   gf = new GridFunction*[nblocks];
+   for (int i=0; i<nblocks; i++)
+   {
+      gf[i] = new GridFunction(&f->GetBlock(i), &data[f->GetVSizeOffset(i)]);
+   }
+}
+
+BlockGridFunction::~BlockGridFunction()
+{
+   for (int i=0; i<nblocks; i++)
+   {
+      delete gf[i];
+   }
+   delete [] gf;
+}
+
+void
+BlockGridFunction::Update()
+{
+   /*
+    FiniteElementSpace * fes = gfr_->FESpace();
+
+    int vsize = fes->GetVSize();
+
+    const Operator *T = fes->GetUpdateOperator();
+    if (T)
+    {
+       // Update the individual GridFunction objects.  This will allocate
+       // new data arrays for each GridFunction.
+       gfr_->Update();
+       gfi_->Update();
+
+       // Our data array now contains old data as well as being the wrong size
+       // so reallocate it.
+       this->SetSize(2 * vsize);
+
+       // Create temporary vectors which point to the new data array
+       Vector gf_r(&data[0], vsize);
+       Vector gf_i(&data[vsize], vsize);
+
+       // Copy the updated GridFunctions into the new data array
+       gf_r = *gfr_;
+       gf_i = *gfi_;
+
+       // Replace the individual data arrays with pointers into the new data array
+       gfr_->NewDataAndSize(&data[0], vsize);
+       gfi_->NewDataAndSize(&data[vsize], vsize);
+    }
+    else
+    {
+       // The existing data will not be transferred to the new GridFunctions
+       // so delete it a allocate a new array
+       this->SetSize(2 * vsize);
+
+       // Point the individual GridFunctions to the new data array
+       gfr_->NewDataAndSize(&data[0], vsize);
+       gfi_->NewDataAndSize(&data[vsize], vsize);
+
+       // These updates will only set the proper 'sequence' value within
+       // the individual GridFunction objects because their sizes are
+       // already correct
+       gfr_->Update();
+       gfi_->Update();
+    }
+   */
+}
+/*
+void
+BlockGridFunction::ProjectCoefficient(Coefficient &real_coeff,
+             Coefficient &imag_coeff)
+{
+ gfr_->ProjectCoefficient(real_coeff);
+ gfi_->ProjectCoefficient(imag_coeff);
+}
+
+void
+BlockGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                      VectorCoefficient &imag_vcoeff)
+{
+ gfr_->ProjectCoefficient(real_vcoeff);
+ gfi_->ProjectCoefficient(imag_vcoeff);
+}
+
+void
+BlockGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
+                                         Coefficient &imag_coeff,
+                                         Array<int> &attr)
+{
+ gfr_->ProjectBdrCoefficient(real_coeff, attr);
+ gfi_->ProjectBdrCoefficient(imag_coeff, attr);
+}
+
+void
+BlockGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
+                                               VectorCoefficient &imag_vcoeff,
+                                               Array<int> &attr)
+{
+ gfr_->ProjectBdrCoefficientNormal(real_vcoeff, attr);
+ gfi_->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
+}
+
+void
+BlockGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
+                                                &real_vcoeff,
+                                                VectorCoefficient
+                                                &imag_vcoeff,
+                                                Array<int> &attr)
+{
+ gfr_->ProjectBdrCoefficientTangent(real_vcoeff, attr);
+ gfi_->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
+}
+*/
+
+BlockLinearForm::BlockLinearForm(BlockFiniteElementSpace *f)
+   : BlockObject(f->GetNBlocks()), Vector(f->GetVSize()), fes(f)
+{
+   lf = new LinearForm*[nblocks];
+   for (int i=0; i<nblocks; i++)
+   {
+      lf[i] = new LinearForm(&f->GetBlock(i), &data[f->GetVSizeOffset(i)]);
+   }
+}
+
+BlockLinearForm::~BlockLinearForm()
+{
+   for (int i=0; i<nblocks; i++)
+   {
+      delete lf[i];
+   }
+   delete [] lf;
+}
+
+void
+BlockLinearForm::AddDomainIntegrator(int index, LinearFormIntegrator *lfi)
+{
+   CheckIndex(index);
+   lf[index]->AddDomainIntegrator(lfi);
+}
+
+void
+BlockLinearForm::Update()
+{
+   /*
+    FiniteElementSpace *fes = lfr_->FESpace();
+
+    this->Update(fes);
+   */
+}
+
+void
+BlockLinearForm::Update(BlockFiniteElementSpace *f)
+{
+   /*
+    int vsize = fes->GetVSize();
+    SetSize(2 * vsize);
+
+    Vector lfr(&data[0], vsize);
+    Vector lfi(&data[vsize], vsize);
+
+    lfr_->Update(fes, lfr, 0);
+    lfi_->Update(fes, lfi, 0);
+   */
+}
+
+void
+BlockLinearForm::Assemble()
+{
+   for (int i=0; i<nblocks; i++)
+   {
+      lf[i]->Assemble();
+   }
+}
+
+double
+BlockLinearForm::operator()(const BlockGridFunction &gf) const
+{
+   double v = 0.0;
+   for (int i=0; i<nblocks; i++)
+   {
+      v += (*lf[i])(gf.GetBlock(i));
+   }
+   return v;
+}
+
+
+BlockBilinearForm::BlockBilinearForm(BlockFiniteElementSpace *f,
+                                     bool symmetric)
+   : TPBlockObject(f->GetNBlocks()),
+     trial_fes(f), test_fes(f), sym(symmetric)
+{
+   blf = new Matrix*[nblocks];
+   // blf = new BilinearForm*[nblocks];
+   // mblf = new MixedBilinearForm*[nblocks];
+   mixed = new bool[nblocks];
+   for (int i=0; i<nblocks; i++)
+   {
+      blf[i] = NULL;
+      mixed[i] = false;
+   }
+}
+
+BlockBilinearForm::~BlockBilinearForm()
+{
+   if ( !sym )
+   {
+      for (int i=0; i<nblocks; i++)
+      {
+         delete blf[i];
+      }
+   }
+   else
+   {
+      for (int r=0; r<nrows; r++)
+      {
+         for (int c=r; c<nrows; c++)
+         {
+            delete blf[r * nrows + c];
+         }
+      }
+   }
+   delete [] blf;
+}
+
+void BlockBilinearForm::initBilinearForm(int r, int c)
+{
+   int index = CheckIndex(r, c);
+   if ( blf[index] != NULL ) { return; }
+
+   if ( &trial_fes->GetBlock(r) == &test_fes->GetBlock(c) )
+   {
+      blf[index] = new BilinearForm(&trial_fes->GetBlock(r));
+   }
+   else
+   {
+      blf[index] = new MixedBilinearForm(&trial_fes->GetBlock(r),
+					 &test_fes->GetBlock(c));
+   }
+   if ( sym && r != c )
+   {
+      int indexT = CheckIndex(c, r);
+      blf[indexT] = blf[index];
+   }
+}
+
+void BlockBilinearForm::AddDomainIntegrator(int r, int c,
+                                            BilinearFormIntegrator *bfi)
+{
+   int index = CheckIndex(r, c);
+   if ( blf[index] == NULL )
+   {
+      initBilinearForm(r, c);
+   }
+   // if (bfi_real) { blfr_->AddDomainIntegrator(bfi_real); }
+   //  if (bfi_imag) { blfi_->AddDomainIntegrator(bfi_imag); }
+}
+
+void
+BlockBilinearForm::AddBoundaryIntegrator(int r, int c,
+                                         BilinearFormIntegrator *bfi)
+{
+   //   if (bfi_real) { blfr_->AddBoundaryIntegrator(bfi_real); }
+   // if (bfi_imag) { blfi_->AddBoundaryIntegrator(bfi_imag); }
+}
+
+void
+BlockBilinearForm::AddBoundaryIntegrator(int r, int c,
+                                         BilinearFormIntegrator *bfi,
+                                         Array<int> & bdr_marker)
+{
+   //  if (bfi_real) { blfr_->AddBoundaryIntegrator(bfi_real, bdr_marker); }
+   // if (bfi_imag) { blfi_->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+BlockBilinearForm::Assemble(int skip_zeros)
+{
+   // blfr_->Assemble(skip_zeros);
+   //  blfi_->Assemble(skip_zeros);
+}
+
+void
+BlockBilinearForm::Finalize(int skip_zeros)
+{
+   // blfr_->Finalize(skip_zeros);
+   // blfi_->Finalize(skip_zeros);
+}
+/*
+BlockSparseMatrix *
+BlockBilinearForm::AssembleCompSpMat()
+{
+// return new BlockSparseMatrix(&blfr_->SpMat(),
+//                              &blfi_->SpMat(),
+//                              false, false, conv_);
+}
+*/
+void
+BlockBilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
+                                    Vector &x, Vector &b,
+                                    OperatorHandle &A,
+                                    Vector &X, Vector &B,
+                                    int ci)
+{
+   /*
+    FiniteElementSpace * fes = blfr_->FESpace();
+
+    int vsize  = fes->GetVSize();
+    // int tvsize = pfes->GetTrueVSize();
+
+    double s = (conv_ == BlockOperator::HERMITIAN)?1.0:-1.0;
+
+    // Allocate temporary vectors
+    Vector b_0(vsize);  b_0 = 0.0;
+    // Vector B_0(tvsize); B_0 = 0.0;
+
+    // Extract the real and imaginary parts of the input vectors
+    MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
+    Vector x_r(x.GetData(), vsize);
+    Vector x_i(&(x.GetData())[vsize], vsize);
+
+    MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
+    Vector b_r(b.GetData(), vsize);
+    Vector b_i(&(b.GetData())[vsize], vsize);
+    b_i *= s;
+
+    SparseMatrix * A_r = new SparseMatrix;
+    SparseMatrix * A_i = new SparseMatrix;
+    Vector X_0, B_0;
+
+    b_0 = b_r;
+    blfr_->FormLinearSystem(ess_tdof_list, x_r, b_r, *A_r, X_0, B_0, ci);
+
+    int tvsize = B_0.Size();
+    X.SetSize(2 * tvsize);
+    B.SetSize(2 * tvsize);
+    Vector X_r(X.GetData(), tvsize);
+    Vector X_i(&(X.GetData())[tvsize], tvsize);
+    Vector B_r(B.GetData(), tvsize);
+    Vector B_i(&(B.GetData())[tvsize], tvsize);
+    X_r = X_0; B_r = B_0;
+
+    b_0 = 0.0;
+    blfi_->FormLinearSystem(ess_tdof_list, x_i, b_0, *A_i, X_0, B_0, false);
+    B_r -= B_0;
+
+    b_0 = b_i;
+    blfr_->FormLinearSystem(ess_tdof_list, x_i, b_0, *A_r, X_0, B_0, ci);
+    X_i = X_0; B_i = B_0;
+
+    b_0 = 0.0;
+    blfi_->FormLinearSystem(ess_tdof_list, x_r, b_0, *A_i, X_0, B_0, false);
+    B_i += B_0;
+
+    B_i *= s;
+    b_i *= s;
+
+    // A = A_r + i A_i
+    A.Clear();
+    BlockSparseMatrix * A_sp =
+       new BlockSparseMatrix(A_r, A_i, true, true, conv_);
+    A.Reset<BlockSparseMatrix>(A_sp, true);
+   */
+}
+
+void
+BlockBilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
+                                      Vector &x)
+{
+   /*
+    FiniteElementSpace * fes = blfr_->FESpace();
+
+    const SparseMatrix *P = fes->GetConformingProlongation();
+
+    int vsize  = fes->GetVSize();
+    int tvsize = X.Size() / 2;
+
+    Vector X_r(X.GetData(), tvsize);
+    Vector X_i(&(X.GetData())[tvsize], tvsize);
+
+    Vector x_r(x.GetData(), vsize);
+    Vector x_i(&(x.GetData())[vsize], vsize);
+
+    if (!P)
+    {
+       x = X;
+    }
+    else
+    {
+       // Apply conforming prolongation
+       P->Mult(X_r, x_r);
+       P->Mult(X_i, x_i);
+    }
+   */
+}
+
+void
+BlockBilinearForm::Update(BlockFiniteElementSpace *nfes)
+{
+   // if ( blfr_ ) { blfr_->Update(nfes); }
+   // if ( blfi_ ) { blfi_->Update(nfes); }
+}
+
+
+#ifdef MFEM_USE_MPI
+
+ParBlockGridFunction::ParBlockGridFunction(ParBlockFiniteElementSpace *pf)
+   : BlockObject(pf->GetNBlocks()), Vector(pf->GetVSize()), pfes(pf)
+{
+   pgf = new ParGridFunction*[nblocks];
+   for (int i=0; i<nblocks; i++)
+   {
+      pgf[i] = new ParGridFunction(&pf->GetBlock(i),
+                                   &data[pf->GetVSizeOffset(i)]);
+   }
+}
+
+void
+ParBlockGridFunction::Update()
+{
+   /*
+    ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
+
+    int vsize = pfes->GetVSize();
+
+    const Operator *T = pfes->GetUpdateOperator();
+    if (T)
+    {
+       // Update the individual GridFunction objects.  This will allocate
+       // new data arrays for each GridFunction.
+       pgfr_->Update();
+       pgfi_->Update();
+
+       // Our data array now contains old data as well as being the wrong size
+       // so reallocate it.
+       this->SetSize(2 * vsize);
+
+       // Create temporary vectors which point to the new data array
+       Vector gf_r(&data[0], vsize);
+       Vector gf_i(&data[vsize], vsize);
+
+       // Copy the updated GridFunctions into the new data array
+       gf_r = *pgfr_;
+       gf_i = *pgfi_;
+
+       // Replace the individual data arrays with pointers into the new data array
+       pgfr_->NewDataAndSize(&data[0], vsize);
+       pgfi_->NewDataAndSize(&data[vsize], vsize);
+    }
+    else
+    {
+       // The existing data will not be transferred to the new GridFunctions
+       // so delete it a allocate a new array
+       this->SetSize(2 * vsize);
+
+       // Point the individual GridFunctions to the new data array
+       pgfr_->NewDataAndSize(&data[0], vsize);
+       pgfi_->NewDataAndSize(&data[vsize], vsize);
+
+       // These updates will only set the proper 'sequence' value within
+       // the individual GridFunction objects because their sizes are
+       // already correct
+       pgfr_->Update();
+       pgfi_->Update();
+    }
+   */
+}
+/*
+void
+ParBlockGridFunction::ProjectCoefficient(Coefficient &real_coeff,
+                                           Coefficient &imag_coeff)
+{
+   pgfr_->ProjectCoefficient(real_coeff);
+   pgfi_->ProjectCoefficient(imag_coeff);
+}
+
+void
+ParBlockGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                           VectorCoefficient &imag_vcoeff)
+{
+   pgfr_->ProjectCoefficient(real_vcoeff);
+   pgfi_->ProjectCoefficient(imag_vcoeff);
+}
+
+void
+ParBlockGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
+                                              Coefficient &imag_coeff,
+                                              Array<int> &attr)
+{
+   pgfr_->ProjectBdrCoefficient(real_coeff, attr);
+   pgfi_->ProjectBdrCoefficient(imag_coeff, attr);
+}
+
+void
+ParBlockGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
+                                                    &real_vcoeff,
+                                                    VectorCoefficient
+                                                    &imag_vcoeff,
+                                                    Array<int> &attr)
+{
+   pgfr_->ProjectBdrCoefficientNormal(real_vcoeff, attr);
+   pgfi_->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
+}
+
+void
+ParBlockGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
+                                                     &real_vcoeff,
+                                                     VectorCoefficient
+                                                     &imag_vcoeff,
+                                                     Array<int> &attr)
+{
+   pgfr_->ProjectBdrCoefficientTangent(real_vcoeff, attr);
+   pgfi_->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
+}
+*/
+void
+ParBlockGridFunction::Distribute(const Vector *tv)
+{
+   /*
+    ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
+    HYPRE_Int size = pfes->GetTrueVSize();
+
+    double * tvd = tv->GetData();
+    Vector tvr(tvd, size);
+    Vector tvi(&tvd[size], size);
+
+    pgfr_->Distribute(tvr);
+    pgfi_->Distribute(tvi);
+   */
+}
+
+void
+ParBlockGridFunction::ParallelProject(Vector &tv) const
+{
+   /*
+    ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
+    HYPRE_Int size = pfes->GetTrueVSize();
+
+    double * tvd = tv.GetData();
+    Vector tvr(tvd, size);
+    Vector tvi(&tvd[size], size);
+
+    pgfr_->ParallelProject(tvr);
+    pgfi_->ParallelProject(tvi);
+   */
+}
+
+
+ParBlockLinearForm::ParBlockLinearForm(ParBlockFiniteElementSpace *pf)
+   : BlockObject(pf->GetNBlocks()), Vector(pf->GetVSize()), pfes(pf)
+{
+   plf = new ParLinearForm*[nblocks];
+   for (int i=0; i<nblocks; i++)
+   {
+      plf[i] = new ParLinearForm(&pf->GetBlock(i),
+                                 &data[pf->GetVSizeOffset(i)]);
+   }
+}
+
+ParBlockLinearForm::~ParBlockLinearForm()
+{
+   for (int i=0; i<nblocks; i++)
+   {
+      delete plf[i];
+   }
+   delete [] plf;
+}
+
+void
+ParBlockLinearForm::AddDomainIntegrator(int index,
+                                        LinearFormIntegrator *lfi)
+{
+   CheckIndex(index);
+   plf[index]->AddDomainIntegrator(lfi);
+}
+
+void
+ParBlockLinearForm::Update(ParBlockFiniteElementSpace *pf)
+{
+   /*
+    ParFiniteElementSpace *pfes = (pf!=NULL)?pf:plfr_->ParFESpace();
+    int vsize = pfes->GetVSize();
+    SetSize(2 * vsize);
+
+    Vector plfr(&data[0], vsize);
+    Vector plfi(&data[vsize], vsize);
+
+    plfr_->Update(pfes, plfr, 0);
+    plfi_->Update(pfes, plfi, 0);
+   */
+}
+
+void
+ParBlockLinearForm::Assemble()
+{
+   /*
+    plfr_->Assemble();
+    plfi_->Assemble();
+    if (conv_ == BlockOperator::BLOCK_SYMMETRIC)
+    {
+       *plfi_ *= -1.0;
+    }
+   */
+}
+
+void
+ParBlockLinearForm::ParallelAssemble(Vector &tv)
+{
+   /*
+    HYPRE_Int size = plfr_->ParFESpace()->GetTrueVSize();
+
+    double * tvd = tv.GetData();
+    Vector tvr(tvd, size);
+    Vector tvi(&tvd[size], size);
+
+    plfr_->ParallelAssemble(tvr);
+    plfi_->ParallelAssemble(tvi);
+   */
+}
+
+BlockVector *
+ParBlockLinearForm::ParallelAssemble()
+{
+   /*
+    const ParFiniteElementSpace * pfes = plfr_->ParFESpace();
+
+    HypreParVector * tv = new HypreParVector(pfes->GetComm(),
+                                             2*(pfes->GlobalTrueVSize()),
+                                             tdof_offsets_);
+
+    HYPRE_Int size = pfes->GetTrueVSize();
+
+    double * tvd = tv->GetData();
+    Vector tvr(tvd, size);
+    Vector tvi(&tvd[size], size);
+
+    plfr_->ParallelAssemble(tvr);
+    plfi_->ParallelAssemble(tvi);
+
+    return tv;
+   */
+}
+
+double
+ParBlockLinearForm::operator()(const ParBlockGridFunction &gf) const
+{
+   double v = 0.0;
+   for (int i=0; i<nblocks; i++)
+   {
+      v += (*plf[i])(gf.GetBlock(i));
+   }
+   return v;
+}
+
+
+ParBlockBilinearForm::ParBlockBilinearForm(ParBlockFiniteElementSpace *pf,
+                                           bool symmetric)
+   : TPBlockObject(pf->GetNBlocks()),
+     trial_pfes(pf), test_pfes(pf), sym(symmetric)
+{
+   pblf = new Matrix*[nblocks];
+   mixed = new bool[nblocks];
+   for (int i=0; i<nblocks; i++)
+   {
+      pblf[i] = NULL;
+      mixed[i] = false;
+   }
+}
+
+ParBlockBilinearForm::~ParBlockBilinearForm()
+{
+   if ( !sym )
+   {
+      for (int i=0; i<nblocks; i++)
+      {
+         delete pblf[i];
+      }
+   }
+   else
+   {
+      for (int r=0; r<nrows; r++)
+      {
+         for (int c=r; c<nrows; c++)
+         {
+            delete pblf[r * nrows + c];
+         }
+      }
+   }
+   delete [] pblf;
+}
+
+void ParBlockBilinearForm::initParBilinearForm(int r, int c)
+{
+   int index = CheckIndex(r, c);
+   if ( pblf[index] != NULL ) { return; }
+
+   if ( &trial_pfes->GetBlock(r) == &test_pfes->GetBlock(c) )
+   {
+      pblf[index] = new ParBilinearForm(&trial_pfes->GetBlock(r));
+   }
+   else
+   {
+      pblf[index] = new ParMixedBilinearForm(&trial_pfes->GetBlock(r),
+					     &test_pfes->GetBlock(c));
+   }
+   if ( sym && r != c )
+   {
+      int indexT = CheckIndex(c, r);
+      pblf[indexT] = pblf[index];
+   }
+}
+
+void ParBlockBilinearForm::AddDomainIntegrator(int r, int c,
+                                               BilinearFormIntegrator *bfi)
+{
+   // if (bfi_real) { pblfr_->AddDomainIntegrator(bfi_real); }
+   // if (bfi_imag) { pblfi_->AddDomainIntegrator(bfi_imag); }
+}
+
+void
+ParBlockBilinearForm::AddBoundaryIntegrator(int r, int c,
+                                            BilinearFormIntegrator *bfi)
+{
+   // if (bfi_real) { pblfr_->AddBoundaryIntegrator(bfi_real); }
+   // if (bfi_imag) { pblfi_->AddBoundaryIntegrator(bfi_imag); }
+}
+
+void
+ParBlockBilinearForm::AddBoundaryIntegrator(int r, int c,
+                                            BilinearFormIntegrator *bfi,
+                                            Array<int> & bdr_marker)
+{
+   // if (bfi_real) { pblfr_->AddBoundaryIntegrator(bfi_real, bdr_marker); }
+   // if (bfi_imag) { pblfi_->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+ParBlockBilinearForm::Assemble(int skip_zeros)
+{
+   // pblfr_->Assemble(skip_zeros);
+   //  pblfi_->Assemble(skip_zeros);
+}
+
+void
+ParBlockBilinearForm::Finalize(int skip_zeros)
+{
+   // pblfr_->Finalize(skip_zeros);
+   // pblfi_->Finalize(skip_zeros);
+}
+
+BlockOperator *
+ParBlockBilinearForm::ParallelAssemble()
+{
+   /*
+   return new BlockHypreParMatrix(pblfr_->ParallelAssemble(),
+                                     pblfi_->ParallelAssemble(),
+                                     true, true, conv_);
+   */
+}
+
+void
+ParBlockBilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
+                                       Vector &x, Vector &b,
+                                       OperatorHandle &A,
+                                       Vector &X, Vector &B,
+                                       int ci)
+{
+   /*
+    ParFiniteElementSpace * pfes = pblfr_->ParFESpace();
+
+    int tvs = pfes->TrueVSize();
+    cout << "TrueVSize returns " << tvs << endl;
+    cout << "GetVSize returns " << pfes->GetVSize() << endl;
+
+    int vsize = x.Size() / 2;
+    // int vsize  = pfes->GetVSize();
+    // int tvsize = pfes->GetTrueVSize();
+
+    cout << "x.Size/2 returns " << vsize << endl;
+
+    double s = (conv_ == BlockOperator::HERMITIAN)?1.0:-1.0;
+
+    // Allocate temporary vectors
+    Vector b_0(vsize);  b_0 = 0.0;
+    // Vector B_0(tvsize); B_0 = 0.0;
+
+    // Extract the real and imaginary parts of the input vectors
+    // MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
+    Vector x_r(x.GetData(), vsize);
+    Vector x_i(&(x.GetData())[vsize], vsize);
+
+    MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
+    Vector b_r(b.GetData(), vsize);
+    Vector b_i(&(b.GetData())[vsize], vsize);
+    b_i *= s;
+
+    OperatorHandle A_r, A_i;
+    Vector X_0, B_0;
+    cout << "pblfr fls 1" << endl << flush;
+    b_0 = b_r;
+    pblfr_->FormLinearSystem(ess_tdof_list, x_r, b_0, A_r, X_0, B_0, ci);
+
+    int tvsize = B_0.Size();
+    X.SetSize(2 * tvsize);
+    B.SetSize(2 * tvsize);
+    Vector X_r(X.GetData(), tvsize);
+    Vector X_i(&(X.GetData())[tvsize], tvsize);
+    Vector B_r(B.GetData(), tvsize);
+    Vector B_i(&(B.GetData())[tvsize], tvsize);
+    X_r = X_0; B_r = B_0;
+    cout << "pblfi fls 1" << endl << flush;
+    b_0 = 0.0;
+    pblfi_->FormLinearSystem(ess_tdof_list, x_i, b_0, A_i, X_0, B_0, false);
+    B_r -= B_0;
+    cout << "pblfr fls 2" << endl << flush;
+    b_0 = b_i;
+    pblfr_->FormLinearSystem(ess_tdof_list, x_i, b_0, A_r, X_0, B_0, ci);
+    X_i = X_0; B_i = B_0;
+    cout << "pblfi fls 2" << endl << flush;
+    b_0 = 0.0;
+    pblfi_->FormLinearSystem(ess_tdof_list, x_r, b_0, A_i, X_0, B_0, false);
+    B_i += B_0;
+
+    B_i *= s;
+    b_i *= s;
+
+    // A = A_r + i A_i
+    A.Clear();
+    if ( A_r.Type() == Operator::Hypre_ParCSR &&
+         A_i.Type() == Operator::Hypre_ParCSR )
+    {
+       BlockHypreParMatrix * A_hyp =
+          new BlockHypreParMatrix(A_r.As<HypreParMatrix>(),
+                                    A_i.As<HypreParMatrix>(),
+                                    A_r.OwnsOperator(),
+                                    A_i.OwnsOperator(),
+                                    conv_);
+       A.Reset<BlockHypreParMatrix>(A_hyp, true);
+    }
+    else
+    {
+       BlockOperator * A_op =
+          new BlockOperator(A_r.As<Operator>(),
+                              A_i.As<Operator>(),
+                              A_r.OwnsOperator(),
+                              A_i.OwnsOperator(),
+                              conv_);
+       A.Reset<BlockOperator>(A_op, true);
+    }
+   */
+}
+
+void
+ParBlockBilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
+                                         Vector &x)
+{
+   /*
+   ParFiniteElementSpace * pfes = pblfr_->ParFESpace();
+
+    const Operator &P = *pfes->GetProlongationMatrix();
+
+    int vsize  = pfes->GetVSize();
+    int tvsize = X.Size() / 2;
+
+    Vector X_r(X.GetData(), tvsize);
+    Vector X_i(&(X.GetData())[tvsize], tvsize);
+
+    Vector x_r(x.GetData(), vsize);
+    Vector x_i(&(x.GetData())[vsize], vsize);
+
+    // Apply conforming prolongation
+    P.Mult(X_r, x_r);
+    P.Mult(X_i, x_i);
+   */
+}
+
+void
+ParBlockBilinearForm::Update(ParBlockFiniteElementSpace *nfes)
+{
+   /*
+    if ( pblfr_ ) { pblfr_->Update(nfes); }
+    if ( pblfi_ ) { pblfi_->Update(nfes); }
+   */
+}
+
+
+#endif // MFEM_USE_MPI
+
+}

fem/block_fem.hpp

@@ -0,0 +1,537 @@
+// 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.
+
+#ifndef MFEM_BLOCK_FEM
+#define MFEM_BLOCK_FEM
+
+#include "gridfunc.hpp"
+#include "linearform.hpp"
+#include "bilinearform.hpp"
+#ifdef MFEM_USE_MPI
+#include "pgridfunc.hpp"
+#include "plinearform.hpp"
+#include "pbilinearform.hpp"
+#endif
+
+namespace mfem
+{
+
+class BlockObject
+{
+protected:
+   int nblocks;
+
+   BlockObject(int num_blocks) : nblocks(num_blocks) {}
+
+   inline void CheckIndex(int index) const
+   {
+      MFEM_ASSERT(index >= 0 && index < nblocks,
+                  "Out of bounds access: " << index << ", size = " << nblocks);
+   }
+
+public:
+
+   inline int GetNBlocks() const { return nblocks; }
+
+};
+
+/// Tensor Product Block Object
+class TPBlockObject
+{
+protected:
+   int nrows;
+   int ncols;
+   int nblocks;
+
+   // Square tensor product
+   TPBlockObject(int size)
+      : nrows(size), ncols(size), nblocks(size * size) {}
+   // Rectangular tensor product
+   TPBlockObject(int num_rows, int num_cols)
+      : nrows(num_rows), ncols(num_cols), nblocks(nrows * ncols) {}
+
+   inline int CheckIndex(int r, int c) const
+   {
+      MFEM_ASSERT(r >= 0 && r < nrows && c >= 0 && c < ncols,
+                  "Out of bounds access: (" << r << "," << c
+                  << "), dimensions = " << nrows << " x " << ncols);
+      return r * ncols + c;
+   }
+
+public:
+   inline int GetNBlocks() const { return nblocks; }
+   inline int GetNRows() const { return nrows; }
+   inline int GetNColumns() const { return ncols; }
+
+};
+
+class BlockFiniteElementSpace : public BlockObject
+{
+private:
+   Array<int> vsize_offsets;
+   FiniteElementSpace **fes;
+   bool owns_fes;
+
+public:
+   /** Construct FiniteElementSpaces with different finite element
+       collections but the same vdim and ordering.  @note The internally
+       constructed FiniteElementSpace objects will be owned, and free'd, by
+       the BlockFiniteElementSpace object. */
+   BlockFiniteElementSpace(int num_blocks, Mesh *mesh,
+                           FiniteElementCollection ** fec,
+                           int vdim = 1, int ordering = Ordering::byNODES);
+
+   /** Construct from an arbitrary set of finite element spaces.
+       @note The FiniteElementSpace objects must be maintained by the caller. */
+   BlockFiniteElementSpace(int num_blocks, FiniteElementSpace ** fespace);
+
+   virtual ~BlockFiniteElementSpace();
+
+   /** The following methods might be useful but let's hold off and see
+       which are actually needed. */
+
+   /// Returns the mesh
+   // inline Mesh *GetMesh() const { return mesh; }
+
+   /// Returns number of degrees of freedom.
+   // inline int GetNDofs() const { return ndofs; }
+
+   /// Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
+   inline int GetVSize() const { return vsize_offsets[nblocks]; }
+
+   /// Return the number of vector true (conforming) dofs.
+   // virtual int GetTrueVSize() const { return GetConformingVSize(); }
+
+   /// Returns the number of conforming ("true") degrees of freedom
+   /// (if the space is on a nonconforming mesh with hanging nodes).
+   // int GetNConformingDofs() const;
+
+   // int GetConformingVSize() const { return vdim * GetNConformingDofs(); }
+
+   void Update();
+
+   inline int GetVSizeOffset(int index) const
+   { return vsize_offsets[index]; }
+
+   FiniteElementSpace & GetBlock(int index)
+   { CheckIndex(index); return *fes[index]; }
+   const FiniteElementSpace & GetBlock(int index) const
+   { CheckIndex(index); return *fes[index]; }
+};
+
+/// Class for grid functions defined in blocks
+class BlockGridFunction : public BlockObject, public Vector
+{
+private:
+   BlockFiniteElementSpace *fes;
+   GridFunction **gf;
+
+public:
+
+   /* @brief Construct a BlockGridFunction associated with the
+      FiniteElementSpace @a *f. */
+   BlockGridFunction(BlockFiniteElementSpace *f);
+
+   void Update();
+
+   /// Assign constant values to the BlockGridFunction data.
+   BlockGridFunction &operator=(double value);
+
+   /*
+    virtual void ProjectCoefficient(Coefficient &real_coeff,
+                                    Coefficient &imag_coeff);
+    virtual void ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                    VectorCoefficient &imag_vcoeff);
+
+    virtual void ProjectBdrCoefficient(Coefficient &real_coeff,
+                                       Coefficient &imag_coeff,
+                                       Array<int> &attr);
+    virtual void ProjectBdrCoefficientNormal(VectorCoefficient &real_coeff,
+                                             VectorCoefficient &imag_coeff,
+                                             Array<int> &attr);
+    virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff,
+                                              VectorCoefficient &imag_coeff,
+                                              Array<int> &attr);
+   */
+   BlockFiniteElementSpace *FESpace() { return fes; }
+   const BlockFiniteElementSpace *FESpace() const { return fes; }
+   /*
+    GridFunction & real() { return *gfr_; }
+    GridFunction & imag() { return *gfi_; }
+    const GridFunction & real() const { return *gfr_; }
+    const GridFunction & imag() const { return *gfi_; }
+   */
+   GridFunction & GetBlock(int index)
+   { CheckIndex(index); return *gf[index]; }
+   const GridFunction & GetBlock(int index) const
+   { CheckIndex(index); return *gf[index]; }
+
+   /// Destroys grid function.
+   virtual ~BlockGridFunction();
+};
+
+class BlockLinearForm : public BlockObject, public Vector
+{
+private:
+   BlockFiniteElementSpace *fes;
+   LinearForm **lf;
+
+public:
+
+   BlockLinearForm(BlockFiniteElementSpace *f);
+
+   virtual ~BlockLinearForm();
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(int index, LinearFormIntegrator *lfi);
+
+   BlockFiniteElementSpace *FESpace() const { return fes; }
+   /*
+    LinearForm & real() { return *lfr_; }
+    LinearForm & imag() { return *lfi_; }
+    const LinearForm & real() const { return *lfr_; }
+    const LinearForm & imag() const { return *lfi_; }
+   */
+   LinearForm & GetBlock(int index)
+   { CheckIndex(index); return *lf[index]; }
+   const LinearForm & GetBlock(int index) const
+   { CheckIndex(index); return *lf[index]; }
+
+   void Update();
+   void Update(BlockFiniteElementSpace *f);
+
+   /// Assembles the linear form i.e. sums over all domain/bdr integrators.
+   void Assemble();
+
+   double operator()(const BlockGridFunction &gf) const;
+
+};
+
+// Class for block-structured bilinear forms
+class BlockBilinearForm : public TPBlockObject
+{
+private:
+   BlockFiniteElementSpace *trial_fes; // Row space
+   BlockFiniteElementSpace *test_fes; // Col space
+   Matrix **blf;
+   bool *mixed;
+   bool sym;
+
+   void initBilinearForm(int r, int c);
+
+public:
+   BlockBilinearForm(BlockFiniteElementSpace *fes, bool symmetric = false);
+   BlockBilinearForm(BlockFiniteElementSpace *trial_fes,
+                     BlockFiniteElementSpace *test_fes);
+   /*
+    BilinearForm & real() { return *blfr_; }
+    BilinearForm & imag() { return *blfi_; }
+    const BilinearForm & real() const { return *blfr_; }
+    const BilinearForm & imag() const { return *blfi_; }
+   */
+   bool isBlockMixed(int r, int c) const
+   { int index = CheckIndex(r, c); return mixed[index]; }
+
+   /** The following methods can return NULL if the requested block is empty */
+   BilinearForm * GetSquareBlock(int r, int c)
+   { int index = CheckIndex(r, c); return (BilinearForm*)blf[index]; }
+   const BilinearForm * GetSquareBlock(int r, int c) const
+   { int index = CheckIndex(r, c); return (BilinearForm*)blf[index]; }
+
+   MixedBilinearForm * GetMixedBlock(int r, int c)
+   { int index = CheckIndex(r, c); return (MixedBilinearForm*)blf[index]; }
+   const MixedBilinearForm * GetMixedBlock(int r, int c) const
+   { int index = CheckIndex(r, c); return (MixedBilinearForm*)blf[index]; }
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(int r, int c, BilinearFormIntegrator *bfi);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(int r, int c, BilinearFormIntegrator *bfi);
+
+   /// Adds new Boundary Integrator, restricted to specific boundary attributes.
+   void AddBoundaryIntegrator(int r, int c, BilinearFormIntegrator *bfi,
+                              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^t A P.
+   /** The returned matrix has to be deleted by the caller. */
+   // BlockSparseMatrix *AssembleCompSpMat();
+
+   /// Return the parallel FE space associated with the ParBilinearForm.
+   BlockFiniteElementSpace *TrialFESpace() const { return trial_fes; }
+   BlockFiniteElementSpace *RowFESpace() const { return trial_fes; }
+   BlockFiniteElementSpace *TestFESpace() const { return test_fes; }
+   BlockFiniteElementSpace *ColumnFESpace() const { return test_fes; }
+
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                         OperatorHandle &A, Vector &X, Vector &B,
+                         int copy_interior = 0);
+
+   /** Call this method after solving a linear system constructed using the
+       FormLinearSystem method to recover the solution as a ParGridFunction-size
+       vector in x. Use the same arguments as in the FormLinearSystem call. */
+   void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
+
+   void Update(BlockFiniteElementSpace *nfes = NULL);
+
+   virtual ~BlockBilinearForm();
+};
+
+#ifdef MFEM_USE_MPI
+
+class ParBlockFiniteElementSpace : public BlockObject
+{
+private:
+   Array<int> vsize_offsets;
+   Array<int> truev_offsets;
+   ParFiniteElementSpace **pfes;
+   bool owns_pfes;
+
+public:
+   /** Construct FiniteElementSpaces with different finite element
+       collections but the same vdim and ordering.  @note The internally
+       constructed FiniteElementSpace objects will be owned, and free'd, by
+       the BlockFiniteElementSpace object. */
+   ParBlockFiniteElementSpace(int num_blocks, ParMesh *mesh,
+                              FiniteElementCollection ** fec,
+                              int vdim = 1, int ordering = Ordering::byNODES);
+
+   /** Construct from an arbitrary set of finite element spaces.
+       @note The FiniteElementSpace objects must be maintained by the caller. */
+   ParBlockFiniteElementSpace(int num_blocks, ParFiniteElementSpace ** fespace);
+
+   virtual ~ParBlockFiniteElementSpace();
+
+   /** The following methods might be useful but let's hold off and see
+       which are actually needed. */
+
+   /// Returns the mesh
+   // inline Mesh *GetMesh() const { return mesh; }
+
+   /// Returns number of degrees of freedom.
+   // inline int GetNDofs() const { return ndofs; }
+
+   /// Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
+   inline int GetVSize() const { return vsize_offsets[nblocks]; }
+
+   /// Return the number of vector true (conforming) dofs.
+   // virtual int GetTrueVSize() const { return GetConformingVSize(); }
+
+   /// Returns the number of conforming ("true") degrees of freedom
+   /// (if the space is on a nonconforming mesh with hanging nodes).
+   // int GetNConformingDofs() const;
+
+   // int GetConformingVSize() const { return vdim * GetNConformingDofs(); }
+
+   void Update();
+
+   inline int GetVSizeOffset(int index) const
+   { return vsize_offsets[index]; }
+
+   inline int GetTrueVSizeOffset(int index) const
+   { return truev_offsets[index]; }
+
+   ParFiniteElementSpace & GetBlock(int index)
+   { CheckIndex(index); return *pfes[index]; }
+   const ParFiniteElementSpace & GetBlock(int index) const
+   { CheckIndex(index); return *pfes[index]; }
+};
+
+/// Class for complex-valued grid function - Vector with associated FE space.
+class ParBlockGridFunction : public BlockObject, public Vector
+{
+private:
+   ParBlockFiniteElementSpace *pfes;
+   ParGridFunction **pgf;
+
+protected:
+   void Destroy();
+
+public:
+
+   /* @brief Construct a ParBlockGridFunction associated with the
+      ParFiniteElementSpace @a *f. */
+   ParBlockGridFunction(ParBlockFiniteElementSpace *pf);
+
+   void Update();
+
+   /// Assign constant values to the ParBlockGridFunction data.
+   ParBlockGridFunction &operator=(double value);
+
+   /*
+    virtual void ProjectCoefficient(Coefficient &real_coeff,
+                                    Coefficient &imag_coeff);
+    virtual void ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                    VectorCoefficient &imag_vcoeff);
+
+    virtual void ProjectBdrCoefficient(Coefficient &real_coeff,
+                                       Coefficient &imag_coeff,
+                                       Array<int> &attr);
+    virtual void ProjectBdrCoefficientNormal(VectorCoefficient &real_coeff,
+                                             VectorCoefficient &imag_coeff,
+                                             Array<int> &attr);
+    virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff,
+                                              VectorCoefficient &imag_coeff,
+                                              Array<int> &attr);
+   */
+   void Distribute(const Vector *tv);
+   void Distribute(const Vector &tv) { Distribute(&tv); }
+
+   /// Returns the vector restricted to the true dofs.
+   void ParallelProject(Vector &tv) const;
+
+   ParBlockFiniteElementSpace *ParFESpace() { return pfes; }
+   const ParBlockFiniteElementSpace *ParFESpace() const { return pfes; }
+   /*
+    ParGridFunction & real() { return *pgfr_; }
+    ParGridFunction & imag() { return *pgfi_; }
+    const ParGridFunction & real() const { return *pgfr_; }
+    const ParGridFunction & imag() const { return *pgfi_; }
+   */
+   ParGridFunction & GetBlock(int index)
+   { CheckIndex(index); return *pgf[index]; }
+   const ParGridFunction & GetBlock(int index) const
+   { CheckIndex(index); return *pgf[index]; }
+
+   /// Destroys grid function.
+   virtual ~ParBlockGridFunction() { Destroy(); }
+
+};
+
+class ParBlockLinearForm : public BlockObject, public Vector
+{
+private:
+   ParBlockFiniteElementSpace *pfes;
+   ParLinearForm ** plf;
+
+public:
+
+   ParBlockLinearForm(ParBlockFiniteElementSpace *pf);
+
+   virtual ~ParBlockLinearForm();
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(int index, LinearFormIntegrator *lfi);
+
+   ParBlockFiniteElementSpace *ParFESpace() const { return pfes; }
+   /*
+    ParLinearForm & real() { return *plfr_; }
+    ParLinearForm & imag() { return *plfi_; }
+    const ParLinearForm & real() const { return *plfr_; }
+    const ParLinearForm & imag() const { return *plfi_; }
+   */
+   ParLinearForm & GetBlock(int index)
+   { CheckIndex(index); return *plf[index]; }
+   const ParLinearForm & GetBlock(int index) const
+   { CheckIndex(index); return *plf[index]; }
+
+   void Update(ParBlockFiniteElementSpace *pf = NULL);
+
+   /// Assembles the linear form i.e. sums over all domain/bdr integrators.
+   void Assemble();
+
+   /// Assemble the vector on the true dofs, i.e. P^t v.
+   void ParallelAssemble(Vector &tv);
+
+   /// Returns the vector assembled on the true dofs, i.e. P^t v.
+   BlockVector *ParallelAssemble();
+
+   double operator()(const ParBlockGridFunction &gf) const;
+
+};
+
+// Class for parallel block-structured bilinear forms
+class ParBlockBilinearForm : public TPBlockObject
+{
+private:
+   ParBlockFiniteElementSpace *trial_pfes; // Row space
+   ParBlockFiniteElementSpace *test_pfes; // Col space
+   Matrix **pblf;
+   bool *mixed;
+   bool sym;
+
+   void initParBilinearForm(int r, int c);
+
+public:
+   ParBlockBilinearForm(ParBlockFiniteElementSpace *pf, bool symmetric = false);
+   ParBlockBilinearForm(ParBlockFiniteElementSpace *trial_pf,
+                        ParBlockFiniteElementSpace *test_pf);
+   /*
+    ParBilinearForm & real() { return *pblf[index]; }
+    ParBilinearForm & imag() { return *pblfi_; }
+    const ParBilinearForm & real() const { return *pblfr_; }
+    const ParBilinearForm & imag() const { return *pblfi_; }
+   */
+   bool isBlockMixed(int r, int c) const
+   { int index = CheckIndex(r, c); return mixed[index]; }
+  
+   /** The following methods can return NULL if the requested block is empty */
+   ParBilinearForm * GetSquareBlock(int r, int c)
+   { int index = CheckIndex(r, c); return (ParBilinearForm*)pblf[index]; }
+   const ParBilinearForm * GetSquareBlock(int r, int c) const
+   { int index = CheckIndex(r, c); return (ParBilinearForm*)pblf[index]; }
+
+   ParMixedBilinearForm * GetMixedBlock(int r, int c)
+   { int index = CheckIndex(r, c); return (ParMixedBilinearForm*)pblf[index]; }
+   const ParMixedBilinearForm * GetMixedBlock(int r, int c) const
+   { int index = CheckIndex(r, c); return (ParMixedBilinearForm*)pblf[index]; }
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(int r, int c, BilinearFormIntegrator *bfi);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(int r, int c, BilinearFormIntegrator *bfi);
+
+   /// Adds new Boundary Integrator, restricted to specific boundary attributes.
+   void AddBoundaryIntegrator(int r, int c, BilinearFormIntegrator *bfi,
+                              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^t A P.
+   /** The returned matrix has to be deleted by the caller. */
+   BlockOperator *ParallelAssemble();
+
+   /// Return the parallel FE space associated with the ParBilinearForm.
+   ParBlockFiniteElementSpace *TrialParFESpace() const { return trial_pfes; }
+   ParBlockFiniteElementSpace *RowParFESpace() const { return trial_pfes; }
+   ParBlockFiniteElementSpace *TestParFESpace() const { return test_pfes; }
+   ParBlockFiniteElementSpace *ColumnParFESpace() const { return test_pfes; }
+
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                         OperatorHandle &A, Vector &X, Vector &B,
+                         int copy_interior = 0);
+
+   /** Call this method after solving a linear system constructed using the
+       FormLinearSystem method to recover the solution as a ParGridFunction-size
+       vector in x. Use the same arguments as in the FormLinearSystem call. */
+   void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
+
+   void Update(ParBlockFiniteElementSpace *nfes = NULL);
+
+   virtual ~ParBlockBilinearForm();
+};
+
+#endif // MFEM_USE_MPI
+
+}
+
+#endif // MFEM_BLOCK_FEM

fem/linearform.hpp

@@ -66,6 +66,13 @@ public:
    LinearForm(FiniteElementSpace *f) : Vector(f->GetVSize())
    { fes = f; extern_lfs = 0; }
 
+   /// Creates linear form associated with FE space @a *f and data @a data.
+   /** The pointers @f and @a data are not owned by the newly constructed
+       object. */ 
+   LinearForm (FiniteElementSpace * f, double * data)
+      : Vector (data, f -> GetVSize())
+   { fes = f; extern_lfs = 0; }
+
    /** @brief Create a LinearForm on the FiniteElementSpace @a f, using the
        same integrators as the LinearForm @a lf.
 

fem/plinearform.hpp

@@ -45,6 +45,12 @@ public:
    /** The pointer @a pf is not owned by the newly constructed object. */
    ParLinearForm(ParFiniteElementSpace *pf) : LinearForm(pf) { pfes = pf; }
 
+   /// Creates ParLinearForm associated with FE space @a *f and data @a data.
+   /** The pointers @f and @a data are not owned by the newly constructed
+       object. */ 
+   ParLinearForm(ParFiniteElementSpace *pf, double * data)
+      : LinearForm(pf, data) { pfes = pf; }
+
    /** @brief Create a ParLinearForm on the ParFiniteElementSpace @a *pf, using
        the same integrators as the ParLinearForm @a *plf.