Merge branch 'master' into DST-maxwell-solver-dev

examples/maxwell-solver-dev/ParDST/DofMapsDST.cpp

@@ -0,0 +1,960 @@
+#include "DofMapsDST.hpp"
+
+double testcoeff(const Vector & x)
+{
+   return sin(3*M_PI*(x.Sum()));
+}
+
+int get_rank(int tdof, std::vector<int> & tdof_offsets)
+{
+   int size = tdof_offsets.size();
+   if (size == 1) { return 0; }
+   std::vector<int>::iterator up;
+   up=std::upper_bound(tdof_offsets.begin(), tdof_offsets.end(),tdof); // 
+   return std::distance(tdof_offsets.begin(),up)-1;
+}
+
+void ComputeTdofOffsets(const MPI_Comm & comm, const ParFiniteElementSpace * pfes, 
+                        std::vector<int> & tdof_offsets)
+{
+   int num_procs;
+   MPI_Comm_size(comm, &num_procs);
+   tdof_offsets.resize(num_procs);
+   int mytoffset = pfes->GetMyTDofOffset();
+   MPI_Allgather(&mytoffset,1,MPI_INT,&tdof_offsets[0],1,MPI_INT,comm);
+}
+
+void GetSubdomainijk(int ip, const Array<int> nxyz, Array<int> & ijk)
+{
+   ijk.SetSize(3);
+   ijk[2] = ip/(nxyz[0]*nxyz[1]);
+   ijk[1] = (ip-ijk[2]*nxyz[0]*nxyz[1])/nxyz[0];
+   ijk[0] = (ip-ijk[2]*nxyz[0]*nxyz[1])%nxyz[0];
+}
+void GetDirectionijk(int id, Array<int> & ijk)
+{
+   ijk.SetSize(3);
+   int n = 3;
+   ijk[2] = id/(n*n) - 1;
+   ijk[1] = (id-(ijk[2]+1)*n*n)/n - 1;
+   ijk[0] = (id-(ijk[2]+1)*n*n)%n - 1;
+}
+
+int GetSubdomainId(const Array<int> nxyz, Array<int> & ijk)
+{
+   int dim=ijk.Size();
+   int k = (dim==2)? 0 : ijk[2];
+   return k*nxyz[1]*nxyz[0] + ijk[1]*nxyz[0] + ijk[0];
+}
+
+int GetDirectionId(const Array<int> & ijk)
+{
+   int n = 3;
+   int dim = ijk.Size();
+   int k = (dim == 2) ? -1 : ijk[2];
+   return (k+1)*n*n + (ijk[1]+1)*n + ijk[0]+1; 
+}
+
+void DofMaps::Init()
+{
+   comm = pfes->GetComm();
+   MPI_Comm_size(comm, &num_procs);
+   MPI_Comm_rank(comm, &myid);
+
+   dim = pfes->GetParMesh()->Dimension();
+   ComputeTdofOffsets(comm, pfes, tdof_offsets);
+   myelemoffset = part->myelem_offset;
+   mytoffset = pfes->GetMyTDofOffset();
+   subdomain_rank = part->subdomain_rank;
+   nrsubdomains = part->nrsubdomains;
+   nxyz.SetSize(3);
+   for (int i = 0; i<3; i++) { nxyz[i] = part->nxyz[i]; }
+
+   //compute sign factors for tdofs
+   int lsize = pfes->GetVSize();
+   int tsize = pfes->GetTrueVSize();
+   tdof_sign.SetSize(tsize); 
+   for (int i = 0; i<lsize; i++)
+   {
+      int j = pfes->GetGlobalTDofNumber(i);
+      if (j<mytoffset || j>=mytoffset+tsize) continue;
+      tdof_sign[j-mytoffset] = pfes->GetDofSign(i);
+   }
+}
+
+DofMaps::DofMaps(ParFiniteElementSpace *pfes_, ParMeshPartition * part_, bool CompFlag_)
+: pfes(pfes_), part(part_), CompFlag(CompFlag_)
+{
+   Init();
+   Setup();
+}
+
+void DofMaps::Setup()
+{
+   // Setup the local FiniteElementSpaces
+   const FiniteElementCollection * fec = pfes->FEColl();
+   fes.SetSize(nrsubdomains);
+   for (int i = 0; i<nrsubdomains; i++)
+   {
+      fes[i] = nullptr; // initialize with null on all procs
+      if (myid == subdomain_rank[i])
+      {
+         fes[i] = new FiniteElementSpace(part->subdomain_mesh[i],fec);
+      }
+   }
+   // cout << "Computing Overlap Tdofs" << endl;
+   SubdomainToSubdomainMapsSetup();
+   // TestSubdomainToSubdomainMaps();
+
+   SubdomainToGlobalMapsSetup();
+   // TestSubdomainToGlobalMaps();
+}
+
+void DofMaps::SubdomainToSubdomainMapsSetup()
+{
+   ComputeOvlpElems();
+   ComputeOvlpTdofs();
+}
+
+void DofMaps::AddElementToOvlpLists(int l, int iel, 
+                            const Array<bool> & neg, const Array<bool> & pos)
+{
+   int kbeg = (dim == 2) ? 0 : -1;
+   int kend = (dim == 2) ? 0 :  1;
+   Array<int> dijk(3);
+   for (int k = kbeg; k<=kend; k++)
+   {
+      if (dim == 3)
+      {
+         if (k == -1 && !neg[2]) continue;   
+         if (k ==  1 && !pos[2]) continue;   
+      }   
+
+      for (int j = -1; j<=1; j++)
+      {
+         if (j== -1 && !neg[1]) continue;   
+         if (j==  1 && !pos[1]) continue;   
+         for (int i = -1; i<=1; i++)
+         {
+            // cases to skip 
+            if (i==-1 && !neg[0]) continue;   
+            if (i== 1 && !pos[0]) continue;   
+
+            if (i==0 && j==0 && k == 0) continue;
+            dijk[0] = i; dijk[1] = j; dijk[2] = (dim==2)?-1 : k;  
+            int DirId = GetDirectionId(dijk);
+            OvlpElems[l][DirId].Append(iel);
+         }
+      }
+   }
+}
+
+void DofMaps::ComputeOvlpElems()
+{
+   // first compute the element in the overlaps
+   OvlpElems.resize(nrsubdomains);
+   int nlayers = 2*part->OvlpNlayers; 
+   // loop through subdomains
+   for (int l = 0; l<nrsubdomains; l++)
+   {
+      if (myid == subdomain_rank[l])
+      {
+         Array<int> ijk;
+         GetSubdomainijk(l,nxyz,ijk);
+         Mesh * mesh = part->subdomain_mesh[l];
+         OvlpElems[l].resize(pow(3,dim)); 
+         Vector pmin, pmax;
+         mesh->GetBoundingBox(pmin,pmax);
+         double h = part->MeshSize;
+         // loop through the elements in the mesh and assign them to the 
+         // appropriate lists of overlaps
+         for (int iel=0; iel< mesh->GetNE(); iel++)
+         {
+            // Get element center
+            Vector center(dim);
+            int geom = mesh->GetElementBaseGeometry(iel);
+            ElementTransformation * tr = mesh->GetElementTransformation(iel);
+            tr->Transform(Geometries.GetCenter(geom),center);
+
+            Array<bool> pos(dim); pos = false;
+            Array<bool> neg(dim); neg = false;
+            // loop through dimensions   
+            for (int d=0;d<dim; d++)
+            {
+               if (ijk[d]>0 && center[d] < pmin[d]+h*nlayers)
+               {
+                  neg[d] = true;
+               }
+
+               if (ijk[d]<nxyz[d]-1 && center[d] > pmax[d]-h*nlayers) 
+               {
+                  pos[d] = true;
+               }
+            }
+            // Add the element to the appropriate lists
+            AddElementToOvlpLists(l,iel,neg,pos);        
+         }   
+      }
+   }
+}
+
+void DofMaps::ComputeOvlpTdofs()
+{
+   OvlpTDofs.resize(nrsubdomains);
+   int nrneighbors = pow(3,dim); // including its self
+
+   // loop through subdomains
+   for (int l = 0; l<nrsubdomains; l++)
+   {
+      if (myid != subdomain_rank[l]) continue;
+      int ntdofs = fes[l]->GetTrueVSize();
+      Array<int> tdof_marker(ntdofs); 
+      OvlpTDofs[l].resize(nrneighbors);
+      // loop through neighboring directions/neighbors
+      for (int d=0; d<nrneighbors; d++)
+      {
+         tdof_marker = 0;
+         Array<int> tdoflist;
+         // Get the direction
+         Array<int> dijk;
+         GetDirectionijk(l,dijk);
+         int nel = OvlpElems[l][d].Size();
+         Array<int>Elems = OvlpElems[l][d];
+         for (int iel = 0; iel<nel; ++iel)
+         {
+            int jel = Elems[iel];   
+            Array<int> ElemDofs;
+
+            fes[l]->GetElementDofs(jel,ElemDofs);
+            int ndof = ElemDofs.Size();
+            for (int i = 0; i<ndof; ++i)
+            {
+               int dof_ = ElemDofs[i];
+               int dof = (dof_ >= 0) ? dof_ : abs(dof_) - 1;
+               if (!tdof_marker[dof])
+               {
+                  tdoflist.Append(dof); // dofs of ip0 in ovlp
+                  tdof_marker[dof] = 1;
+               }
+            }
+         }
+         OvlpTDofs[l][d] = tdoflist; 
+         if (CompFlag)
+         {
+            for (int i=0; i<tdoflist.Size(); i++)
+            {
+               tdoflist[i] += fes[l]->GetTrueVSize();
+            }
+            OvlpTDofs[l][d].Append(tdoflist);
+         }
+      }
+   }
+}
+
+void DofMaps::PrintOvlpTdofs()
+{
+   int nrneighbors = pow(3,dim); // including its self
+   if (myid == 0)
+   {
+      for (int i = 0; i<nrsubdomains; i++)
+      {
+         if (myid != subdomain_rank[i]) continue;
+         Array<int> ijk;
+         GetSubdomainijk(i,nxyz,ijk);
+         cout << "subdomain = " ; ijk.Print();
+         cout << "myid = " << myid << endl;
+         cout << "ip   = " << i << endl;
+         for (int d = 0; d<nrneighbors; d++)
+         {
+            Array<int> dijk;
+            GetDirectionijk(d,dijk);
+            cout << "direction = " ; dijk.Print();
+
+            if (OvlpTDofs[i][d].Size())
+            {
+               cout << "OvlpTdofs = " ; 
+               OvlpTDofs[i][d].Print(cout,OvlpTDofs[i][d].Size() );
+            }
+         }
+      }
+   }
+}
+
+void DofMaps::TransferToNeighbors(const Array<int> & SubdomainIds, const Array<Vector *> & x,
+std::vector<std::vector<Vector * >> & OvlpSol)
+{
+   // 2D for now....
+   MFEM_VERIFY(SubdomainIds.Size() == x.Size(), "TransferToNeighbors: Size inconsistency");
+   int nrsendIds = SubdomainIds.Size();
+   int nrneighbors = pow(3,dim);
+   MPI_Request *recv_requests = new MPI_Request[nrsendIds*nrneighbors];
+   MPI_Request *send_requests = new MPI_Request[nrsendIds*nrneighbors];
+   MPI_Status  *recv_statuses = new MPI_Status[nrsendIds*nrneighbors];
+   MPI_Status  *send_statuses = new MPI_Status[nrsendIds*nrneighbors];
+   Array<Vector * > send_buffer(nrsendIds*nrneighbors);
+   Array<Vector * > recv_buffer(nrsendIds*nrneighbors);
+   int send_counter = 0;
+   int recv_counter = 0;
+   for (int is = 0; is<nrsendIds; is++)
+   {
+      int i0 = SubdomainIds[is];
+      Array<int> ijk;
+      GetSubdomainijk(i0,nxyz,ijk);
+      for (int d=0;d<nrneighbors; d++)
+      {
+         Array<int>directions;
+         GetDirectionijk(d,directions);
+
+         if (dim == 2 && directions[0] == 0 && directions[1] == 0) continue;
+         if (dim == 3 && directions[0] == 0 
+                      && directions[1] == 0 
+                      && directions[2] == 0) continue;
+         int i = ijk[0] + directions[0];
+         if (i<0 || i>=nxyz[0]) continue;
+         int j = ijk[1] + directions[1];
+         if (j<0 || j>=nxyz[1]) continue;
+         int k = (dim ==3 ) ? ijk[2] + directions[2] : 0;
+         if (k<0 || k>=nxyz[2]) continue;
+         Array<int>ijk1(3);
+         ijk1[0] = i;
+         ijk1[1] = j;
+         ijk1[2] = k;
+         int i1 = GetSubdomainId(nxyz,ijk1);
+         if (myid == subdomain_rank[i0])
+         {
+            Array<int> tdofs0 = OvlpTDofs[i0][d]; // map of dofs in the overlap
+            send_buffer[send_counter] = new Vector(tdofs0.Size());
+            x[is]->GetSubVector(tdofs0,*send_buffer[send_counter]);
+            // Destination rank
+            int dest = subdomain_rank[i1];
+            int tag = i0 * nrneighbors + d;
+
+            int count = tdofs0.Size();
+            MPI_Isend(send_buffer[send_counter]->GetData(),count,MPI_DOUBLE,dest,
+                      tag,comm,&send_requests[send_counter]);
+            send_counter++;
+
+         }
+         if (myid == subdomain_rank[i1])
+         {
+            Array<int> direction1(3); direction1 = -1;
+            for (int dd=0;dd<dim;dd++) 
+            { 
+               direction1[dd] = -directions[dd];
+            }
+            int d1 = GetDirectionId(direction1);
+
+            int count = OvlpTDofs[i1][d1].Size();
+            recv_buffer[recv_counter] = new Vector(count); 
+            int src = subdomain_rank[i0];
+            int tag = i0 * nrneighbors + d;
+            MPI_Irecv(recv_buffer[recv_counter]->GetData(), count,MPI_DOUBLE,src,
+                      tag,comm, &recv_requests[recv_counter]);
+            recv_counter++;
+         }
+      }   
+   }
+   MPI_Waitall(send_counter, send_requests, send_statuses);
+   MPI_Waitall(recv_counter, recv_requests, recv_statuses);
+
+   delete [] send_statuses;
+   delete [] send_requests;
+   delete [] recv_statuses;
+   delete [] recv_requests;
+
+   for (int i = 0; i<send_counter; i++)
+   {
+      delete send_buffer[i];
+   }
+   send_buffer.DeleteAll();
+
+
+   // Extract the transfered solutions
+   recv_counter = 0;
+   for (int is = 0; is<nrsendIds; is++)
+   {
+      int i0 = SubdomainIds[is];
+      Array<int> ijk;
+      GetSubdomainijk(i0,nxyz,ijk);
+      for (int d=0;d<nrneighbors; d++)
+      {
+         Array<int>directions;
+         GetDirectionijk(d,directions);
+         if (dim == 2 && directions[0] == 0 && directions[1] == 0) continue;
+         if (dim == 3 && directions[0] == 0 
+                      && directions[1] == 0 
+                      && directions[2] == 0) continue;
+         int i = ijk[0] + directions[0];
+         if (i<0 || i>=nxyz[0]) continue;
+         int j = ijk[1] + directions[1];
+         if (j<0 || j>=nxyz[1]) continue;
+         int k = (dim ==3 ) ? ijk[2] + directions[2] : 0;
+         if (k<0 || k>=nxyz[2]) continue;
+
+         Array<int>ijk1(3);
+         ijk1[0] = i;
+         ijk1[1] = j;
+         ijk1[2] = k;
+         int i1 = GetSubdomainId(nxyz,ijk1);
+         if (myid == subdomain_rank[i1])
+         {
+            Array<int> direction1(3); direction1 = -1;
+            for (int d=0;d<dim;d++) 
+            { 
+               direction1[d] = -directions[d];
+            }
+            int d1 = GetDirectionId(direction1);
+            Array<int> tdofs1 = OvlpTDofs[i1][d1];
+            if (!OvlpSol[i1][d1])
+            {
+               OvlpSol[i1][d1] = new Vector(2*fes[i1]->GetTrueVSize());
+            }
+            *OvlpSol[i1][d1] = 0.0;
+            OvlpSol[i1][d1]->SetSubVector(tdofs1,*recv_buffer[recv_counter]);
+            recv_counter++;
+         }
+      }   
+   }
+   for (int i = 0; i<recv_counter; i++)
+   {
+      delete recv_buffer[i];
+   }
+   recv_buffer.DeleteAll();
+}
+
+void DofMaps::TestSubdomainToSubdomainMaps()
+{
+      // testing inter-subdomain communication
+   FunctionCoefficient c1(testcoeff);
+   int nrsub = nrsubdomains;
+   Array<int> subdomain_ids(nrsub);
+   Array<Vector*> x(nrsub);
+   for (int i = 0; i<nrsub; i++)
+   {
+      x[i] = nullptr;
+      subdomain_ids[i] = i;
+      if (fes[i])
+      {
+         ComplexGridFunction gf(fes[i]);
+         gf = 0.0;
+         gf.ProjectCoefficient(c1,c1);
+         x[i] = new Vector(2*fes[i]->GetTrueVSize());
+         *x[i] = gf;
+      }
+   }
+
+   std::vector<std::vector<Vector * >> OvlpSol;
+
+   OvlpSol.resize(nrsubdomains);
+   int nrneighbors = pow(3,dim);
+   for (int ip = 0; ip<nrsubdomains; ip++)
+   {
+      if (myid == subdomain_rank[ip])
+      {
+         OvlpSol[ip].resize(nrneighbors);
+      }
+   }
+
+
+   TransferToNeighbors(subdomain_ids,x,OvlpSol);
+
+   string keys = "keys amrRljc\n";
+   for (int i0 = 0 ; i0< nrsubdomains; i0++)
+   {
+      if (fes[i0])
+      {
+         ComplexGridFunction gf0(fes[i0]);
+         for (int d = 0; d<nrneighbors; d++)
+         {
+            if(OvlpSol[i0][d])
+            {
+               Array<int>dijk;
+               GetDirectionijk(d,dijk);
+               Array<int>ijk;
+               GetSubdomainijk(i0,nxyz,ijk);
+               ostringstream oss;
+               oss << "myid: " << myid 
+                     << ", subdomain: (" << ijk[0] << "," << ijk[1] <<")"
+                     << ", direction: (" << dijk[0] << "," << dijk[1] <<")";
+
+               gf0 = 0.0;
+               gf0.real().SetVector(*OvlpSol[i0][d],0);
+               gf0.imag().SetVector(*OvlpSol[i0][d],fes[i0]->GetTrueVSize());
+               char vishost[] = "localhost";
+               int  visport   = 19916;
+               socketstream sol_sock(vishost, visport);
+               sol_sock.precision(8);
+               sol_sock << "solution\n" << *(part->subdomain_mesh[i0]) << gf0.real() 
+               << keys
+               << "window_title '" << oss.str() << "'" << flush;
+            }
+         }
+      }
+   }
+   for (int i = 0; i<nrsub; i++)
+   {
+      delete x[i];
+   }
+}
+
+void DofMaps::SubdomainToGlobalMapsSetup()
+{
+   // workspace for MPI_AlltoAll
+   send_count.SetSize(num_procs);  send_count = 0;
+   send_displ.SetSize(num_procs);  send_displ = 0; 
+   recv_count.SetSize(num_procs);  recv_count = 0;
+   recv_displ.SetSize(num_procs);  recv_displ = 0;
+
+   // 1. Communicate to the subdomain rank the list of tdofs 
+   // a. Compute send count
+   for (int ip = 0; ip<nrsubdomains; ++ip)
+   {
+      // avoid any communication if on subdomain rank
+      int nel = part->local_element_map[ip].Size();
+
+      for (int iel = 0; iel<nel; iel++)
+      {
+         int elem_idx = part->local_element_map[ip][iel] - myelemoffset;
+         // int ndofs = local_tdofs[ip].Size();
+         int ndofs = pfes->GetFE(elem_idx)->GetDof();
+
+         send_count[subdomain_rank[ip]] += 2 + ndofs; 
+      }
+   }
+   // b. Compute receive count
+   MPI_Alltoall(send_count,1,MPI_INT,recv_count,1,MPI_INT,comm);
+   for (int k=0; k<num_procs-1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+   }
+   sbuff_size = send_count.Sum();
+   rbuff_size = recv_count.Sum();
+   // c. Allocate and fill the send buffer
+   Array<int> sendbuf(sbuff_size);  sendbuf = 0;
+   Array<int> soffs(num_procs); soffs = 0;
+   for (int ip = 0; ip<nrsubdomains; ++ip)
+   {
+      int nel = part->local_element_map[ip].Size();
+      for (int iel = 0; iel<nel; iel++)
+      {
+         int elem_idx = part->local_element_map[ip][iel] - myelemoffset;
+         Array<int>ElemDofs;
+         pfes->GetElementDofs(elem_idx,ElemDofs);
+         int ndofs = ElemDofs.Size();
+
+         int j = send_displ[subdomain_rank[ip]] + soffs[subdomain_rank[ip]];
+         sendbuf[j] = ip;
+         sendbuf[j+1] = ndofs;
+
+         for (int k = 0; k < ndofs ; ++k)
+         {
+            int edof_ = ElemDofs[k];
+            int edof = (edof_ >= 0) ? edof_ : abs(edof_) - 1;
+            sendbuf[j+2+k] = pfes->GetGlobalTDofNumber(edof);
+         }
+         soffs[subdomain_rank[ip]] +=  2 + ndofs;
+      }
+   }
+
+   // d. Communication
+   Array<int> recvbuf(rbuff_size);
+   MPI_Alltoallv(sendbuf, send_count, send_displ, MPI_INT, recvbuf,
+                 recv_count, recv_displ, MPI_INT, comm);
+
+   // 3. Extract from recv_buffer
+   std::vector<Array<int>> global_tdofs(nrsubdomains);
+   int k=0;
+   while (k<rbuff_size)
+   {
+      int ip = recvbuf[k++]; 
+      int ndofs = recvbuf[k++];
+      for (int i = 0; i < ndofs; ++i)
+      {
+         global_tdofs[ip].Append(recvbuf[i+k]); 
+      }
+      k += ndofs;
+   }
+
+   SubdomainGTrueDofs.resize(nrsubdomains);
+   // 4. Construct SubdomainTdof to Global mesh tdof maps
+   for (int ip=0; ip<nrsubdomains; ++ip)
+   {
+      if (myid != subdomain_rank[ip]) continue;
+      int nrdof = fes[ip]->GetTrueVSize();
+
+      SubdomainGTrueDofs[ip].SetSize(nrdof);
+      int nel = part->element_map[ip].Size();
+      int k = 0;
+      for (int iel = 0; iel<nel; ++iel)
+      {
+         Array<int> elem_dofs;
+         fes[ip]->GetElementDofs(iel,elem_dofs);
+         int ndof = elem_dofs.Size();
+         for (int i = 0; i<ndof; ++i)
+         {
+            int edof_ = elem_dofs[i];
+            int edof = (edof_ >= 0) ? edof_ : abs(edof_) - 1;
+            // rearranging dofs from serial fespace to pfes ordering
+            SubdomainGTrueDofs[ip][edof] = global_tdofs[ip][k++];
+         }
+      }
+   }
+
+   // 5. Communicate SubdomainGTrueDofs to participating ranks 
+   send_count = 0;  send_displ = 0;
+   recv_count = 0;  recv_displ = 0;
+
+   for (int ip = 0; ip < nrsubdomains; ++ip)
+   {
+      if (myid != subdomain_rank[ip]) continue;
+      int ndofs = SubdomainGTrueDofs[ip].Size();
+      for (int i = 0; i<ndofs; ++i)
+      {
+         int tdof = SubdomainGTrueDofs[ip][i];
+         int rank = get_rank(tdof,tdof_offsets);
+         if (rank == subdomain_rank[ip]) continue; // <--------------
+         send_count[rank] += 2; // 1 for the dof and 1 for the ip that goes to 
+      }
+   }
+
+   // communicate so that recv_count is constructed
+   MPI_Alltoall(send_count,1,MPI_INT,recv_count,1,MPI_INT,comm);
+   //
+   for (int k=0; k<num_procs-1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+   }
+   sbuff_size = send_count.Sum();
+   rbuff_size = recv_count.Sum();
+
+   sendbuf.SetSize(sbuff_size);  
+   sendbuf = 0; soffs = 0;
+
+   for (int ip = 0; ip < nrsubdomains; ip++)
+   {
+      if (myid != subdomain_rank[ip]) continue;
+      int ndofs = SubdomainGTrueDofs[ip].Size();
+      // loop through dofs
+      for (int i = 0; i<ndofs; ++i)
+      {
+         int tdof = SubdomainGTrueDofs[ip][i];
+         int irank = get_rank(tdof,tdof_offsets);
+         if (irank == subdomain_rank[ip]) continue; // <--------------
+         int j = send_displ[irank] + soffs[irank];
+         sendbuf[j] = ip;
+         sendbuf[j+1] = SubdomainGTrueDofs[ip][i];
+         soffs[irank] +=  2 ;
+      }
+   }
+
+   recvbuf.SetSize(rbuff_size);
+   MPI_Alltoallv(sendbuf, send_count, send_displ, MPI_INT, recvbuf,
+                 recv_count, recv_displ, MPI_INT, comm);
+
+   // List of tdofs owned by the processor for subdomains not owned
+   SubdomainLTrueDofs.resize(nrsubdomains); 
+   for (int k=0; k<rbuff_size/2; k++)
+   {
+      int ip = recvbuf[2*k];
+      int tdof = recvbuf[2*k+1];
+      SubdomainLTrueDofs[ip].Append(tdof);
+   }    
+}
+
+// Restriction of global residual to subdomain residuals
+void DofMaps::GlobalToSubdomains(const Vector & y, Array<Vector*> & x)
+{
+   send_count = 0; send_displ = 0;
+   recv_count = 0; recv_displ = 0;
+
+   // Compute send_counts
+   int m = (CompFlag) ? 2 : 1 ;
+   for (int ip = 0; ip < nrsubdomains; ip++)
+   {
+      if (myid == subdomain_rank[ip]) continue;  // <---------------
+      int ndofs = SubdomainLTrueDofs[ip].Size();
+      send_count[subdomain_rank[ip]] += m * ndofs;
+   }
+
+    // communicate so that recv_count is constructed
+   MPI_Alltoall(send_count,1,MPI_INT,recv_count,1,MPI_INT,comm);
+
+   for (int k=0; k<num_procs-1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+   }
+   sbuff_size = send_count.Sum();
+   rbuff_size = recv_count.Sum();
+
+   Array<double> sendbuf(sbuff_size); sendbuf = 0;
+   Array<int> soffs(num_procs); soffs = 0;
+
+   for (int ip = 0; ip < nrsubdomains; ip++)
+   {
+      if (myid == subdomain_rank[ip]) continue;  // <---------------
+      int ndofs = SubdomainLTrueDofs[ip].Size();
+      for (int i = 0; i<ndofs; i++)
+      {
+         int tdof = SubdomainLTrueDofs[ip][i];
+         int j = send_displ[subdomain_rank[ip]] + soffs[subdomain_rank[ip]];
+         soffs[subdomain_rank[ip]] +=m;
+         int k = tdof - mytoffset;
+         // sendbuf[j] = y[k];
+         sendbuf[j] = tdof_sign[k]*y[k];
+         if (CompFlag)
+         {  // if complex valued 
+            int tsize = pfes->GetTrueVSize();
+            // sendbuf[j+1] = y[k+tsize];
+            sendbuf[j+1] = tdof_sign[k]*y[k+tsize];
+         }
+      }
+   }
+
+   // communication
+   Array<double> recvbuf(rbuff_size);
+   MPI_Alltoallv(sendbuf, send_count, send_displ, MPI_DOUBLE, recvbuf,
+                 recv_count, recv_displ, MPI_DOUBLE, comm);
+   Array<int> roffs(num_procs);
+   roffs = 0;
+   // Now each process will construct the res vector
+   x.SetSize(nrsubdomains);
+   for (int ip = 0; ip < nrsubdomains; ip++)
+   {
+      if (myid != subdomain_rank[ip]) continue;
+      int ndof = SubdomainGTrueDofs[ip].Size();
+      if (!x[ip]) x[ip] = new Vector(m*ndof); 
+      *x[ip] = 0.0;
+      // extract the data from receiv buffer
+      for (int i=0; i<ndof; i++)
+      {
+         // pick up the tdof and find its rank
+         int tdof = SubdomainGTrueDofs[ip][i];
+         int tdof_rank = get_rank(tdof,tdof_offsets);
+         if (tdof_rank != subdomain_rank[ip])  // <---------------
+         {
+            int k = recv_displ[tdof_rank] + roffs[tdof_rank];
+            roffs[tdof_rank] += m;
+            (*x[ip])[i] = recvbuf[k];
+            if (CompFlag)
+            {  
+               (*x[ip])[i+ndof] = recvbuf[k+1];
+            }
+         }
+         else
+         {
+            int k = tdof - mytoffset;
+            // (*x[ip])[i] = y[k];
+            (*x[ip])[i] = tdof_sign[k]*y[k];
+            if (CompFlag)
+            {  
+               int gtsize = pfes->GetTrueVSize();
+               (*x[ip])[i+ndof] = tdof_sign[k]*y[k+gtsize];
+            }
+         }
+      }
+   }
+}
+
+// Prolongation of subdomain solutions to the global solution
+void DofMaps::SubdomainsToGlobal(const Array<Vector*> & x, Vector & y)
+{
+   send_count = 0; send_displ = 0;
+   recv_count = 0; recv_displ = 0;
+
+   // Compute send_counts
+   int m = (CompFlag) ? 2 : 1 ;
+   for (int ip = 0; ip < nrsubdomains; ip++)
+   {
+      if (myid != subdomain_rank[ip]) continue;
+      int ndofs = SubdomainGTrueDofs[ip].Size();
+      for (int i=0; i<ndofs; i++)
+      {
+         // pick up the tdof and find its rank
+         int tdof = SubdomainGTrueDofs[ip][i];
+         int tdof_rank = get_rank(tdof,tdof_offsets);
+         if (tdof_rank == subdomain_rank[ip]) continue;
+         send_count[tdof_rank] +=m;
+      }
+   }
+
+   MPI_Alltoall(send_count,1,MPI_INT,recv_count,1,MPI_INT,comm);
+
+   for (int k=0; k<num_procs-1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+   }
+   sbuff_size = send_count.Sum();
+   rbuff_size = recv_count.Sum();
+
+   Array<double> sendbuf(sbuff_size); sendbuf = 0;
+   Array<int> soffs(num_procs); soffs = 0;
+
+   for (int ip = 0; ip < nrsubdomains; ip++)
+   {
+      if (myid != subdomain_rank[ip]) continue;
+      int ndofs = SubdomainGTrueDofs[ip].Size();
+      // loop through dofs
+      for (int i=0; i<ndofs; i++)
+      {
+         //  pick up the dof and find its tdof_rank
+         int tdof = SubdomainGTrueDofs[ip][i];
+         int tdof_rank = get_rank(tdof,tdof_offsets);
+         // offset
+         if (tdof_rank == subdomain_rank[ip]) continue;
+         int k = send_displ[tdof_rank] + soffs[tdof_rank];
+         soffs[tdof_rank] +=m;
+         sendbuf[k] = (*x[ip])[i];
+         if (CompFlag)
+         {
+            sendbuf[k+1] = (*x[ip])[i+ndofs];
+         }
+      }
+   }
+
+   Array<double> recvbuf(rbuff_size);
+   Array<int> roffs(num_procs); roffs = 0;
+   MPI_Alltoallv(sendbuf, send_count, send_displ, MPI_DOUBLE, recvbuf,
+                 recv_count, recv_displ, MPI_DOUBLE, comm);
+
+   for (int ip = 0; ip < nrsubdomains; ip++)
+   {
+      if (myid == subdomain_rank[ip]) 
+      {
+         int ndofs = SubdomainGTrueDofs[ip].Size();
+         for (int i = 0; i<ndofs; i++)
+         {
+            int tdof = SubdomainGTrueDofs[ip][i];
+            int k = tdof - mytoffset;
+            if (k<0 || k>=pfes->GetTrueVSize()) continue;
+            y[k] += tdof_sign[k] * (*x[ip])[i];
+            if (CompFlag)
+            {  
+               int gtsize = pfes->GetTrueVSize();
+               y[k+gtsize] += tdof_sign[k]*(*x[ip])[i+ndofs];
+            }
+         }
+      }
+      else
+      {
+         int ndofs = SubdomainLTrueDofs[ip].Size();
+         for (int i = 0; i<ndofs; i++)
+         {
+            int tdof = SubdomainLTrueDofs[ip][i];
+            int k = tdof - mytoffset;
+            int j = recv_displ[subdomain_rank[ip]] + roffs[subdomain_rank[ip]];
+            roffs[subdomain_rank[ip]] +=m;
+            y[k] += tdof_sign[k] * recvbuf[j];
+            if (CompFlag)
+            {
+               int tsize = pfes->GetTrueVSize();
+               y[k+tsize] += tdof_sign[k]*recvbuf[j+1];
+            }
+         }
+      }
+   }              
+}
+
+void DofMaps::TestSubdomainToGlobalMaps()
+{
+   cout << "Testing Subdomain To Global Maps" << endl;
+   FunctionCoefficient c1(testcoeff);
+   Array<Vector*> x(nrsubdomains);
+   Vector y(pfes->GetTrueVSize()); y = 0.0;
+   for (int i = 0 ; i<nrsubdomains; i++)
+   {
+      if (myid != subdomain_rank[i]) continue;
+      x[i] = new Vector(fes[i]->GetTrueVSize());
+      GridFunction gf(fes[i]);
+      gf = 0.0;
+
+      if (i==3) gf.ProjectCoefficient(c1);
+      *x[i] = gf;
+   }
+
+   SubdomainsToGlobal(x,y);
+
+   // cout << "1: myid = " <<  myid << ", y = "; y.Print();
+
+   string keys = (dim==2) ? "keys amrRljc\n": "keys m\n";
+   ParGridFunction pgf(pfes);
+
+   const Operator &P = *pfes->GetProlongationMatrix();
+   P.Mult(y, pgf);
+   
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock(vishost, visport);
+   sol_sock.precision(8);
+   sol_sock << "parallel " << num_procs << " " << myid << "\n"
+               << "solution\n" << *pfes->GetParMesh() << pgf 
+               << keys << flush;
+
+   ParGridFunction pgf1(pfes);
+   pgf1.ProjectCoefficient(c1);
+   Vector y1(pfes->GetTrueVSize());
+   const SparseMatrix * R = pfes->GetRestrictionMatrix();
+
+   R->Mult(pgf1,y1);
+   // P.MultTranspose(pgf1,y1);
+   Array<Vector*> x1;
+   GlobalToSubdomains(y1,x1);
+
+
+   // for (int i = 0 ; i<nrsubdomains; i++)
+   // {
+   //    if (myid != subdomain_rank[i]) continue;
+   //    ostringstream mesh_name;
+   //    mesh_name << "output/mesh." << setfill('0') << setw(6) << i;
+   //    ofstream mesh_ofs(mesh_name.str().c_str());
+   //    mesh_ofs.precision(8);
+   //    fes[i]->GetMesh()->Print(mesh_ofs);
+   //    GridFunction gf(fes[i]);
+   //    gf = x1[i];
+   //    ostringstream gf_name;
+   //    gf_name << "output/gf." << setfill('0') << setw(6) << i;
+   //    ofstream gf_ofs(gf_name.str().c_str());
+   //    gf_ofs.precision(8);
+   //    gf.Save(gf_ofs);
+   // }
+
+
+
+   int nrsub = nrsubdomains;
+   for (int i = 0 ; i<nrsub; i++)
+   {
+      if (myid == subdomain_rank[i])
+      {
+         socketstream sol_sock1(vishost, visport);
+         sol_sock1.precision(8);
+         sol_sock1 << "parallel " << nrsub << " " << i << "\n";
+         GridFunction gf(fes[i]);
+         GridFunction gf1(fes[i]);
+         gf1.ProjectCoefficient(c1);
+         gf = *x1[i];
+         gf1-=gf;
+         cout << "ip, Diff norm = " <<i<<", " << gf1.Norml2() << endl;
+         sol_sock1 << "solution\n" << *fes[i]->GetMesh() << gf
+                   << keys << flush;
+      }
+      MPI_Barrier(MPI_COMM_WORLD);
+   }
+
+   socketstream gf_sock(vishost, visport);
+   gf_sock.precision(8);
+   gf_sock << "parallel " << num_procs << " " << myid << "\n"
+           << "solution\n" << *pfes->GetParMesh() << pgf1 
+           << keys << flush;
+}
+
+
+DofMaps::~DofMaps()
+{
+   for (int i = 0; i<nrsubdomains; i++)
+   {
+      delete fes[i];
+   }
+}

examples/maxwell-solver-dev/ParDST/DofMapsDST.hpp

@@ -0,0 +1,109 @@
+#pragma once
+#include "../common/Utilities.hpp"
+#include "../common/PML.hpp"
+using namespace std;
+using namespace mfem;
+
+double testcoeff(const Vector & x);
+int get_rank(int tdof, std::vector<int> & tdof_offsets);
+
+
+void ComputeTdofOffsets(const MPI_Comm & comm, const ParFiniteElementSpace * pfes, 
+                        std::vector<int> & tdof_offsets);
+
+void GetSubdomainijk(int ip, const Array<int> nxyz, Array<int> & ijk);
+void GetDirectionijk(int id, Array<int> & ijk);
+int GetSubdomainId(const Array<int> nxyz, Array<int> & ijk);
+int GetDirectionId(const Array<int> & ijk);
+
+
+// class handling two types of dof maps
+// 1. Subdomain truedofs ---> Global truedofs
+// 2. Subdomain truedofs ---> Neighbor truedofs 
+class DofMaps
+{
+private:
+   // The FE space of the problem (H1/Hcurl)
+   ParFiniteElementSpace *pfes = nullptr;
+
+   // The given partition of the parmesh
+   ParMeshPartition *part = nullptr;
+   // partition in x-y-z
+   Array<int> nxyz;
+
+   // MPI parameters
+   MPI_Comm comm = MPI_COMM_WORLD;
+   int num_procs, myid;
+
+   // true dof offset and element offset of the processor
+   vector<int> tdof_offsets;
+   int mytoffset;
+   int myelemoffset;
+   
+   int dim;
+   // Total number of subdomains
+   int nrsubdomains;
+   
+   // Array specifying the subdomain rank
+   Array<int> subdomain_rank;
+
+   // Complex flag
+   bool CompFlag;
+
+   // sign factors
+   Array<int> tdof_sign;
+   // Initializing mpi and helper parameters
+   void Init();
+
+   // 1. Setting up the subdomains FE spaces
+   // 2. Setting up the subdomains-to-subdomains maps
+   // 3. Setting up the subdomain-to-global maps
+   void Setup();
+
+   // -----------------------------------------------
+   //           Subdomain to Subdomain maps
+   // -----------------------------------------------
+   std::vector<std::vector<Array<int>>> OvlpElems;
+   void AddElementToOvlpLists(int l, int iel, 
+                            const Array<bool> & neg, 
+                            const Array<bool> & pos);
+   std::vector<std::vector<Array<int>>> OvlpTDofs;
+   void SubdomainToSubdomainMapsSetup();
+   void ComputeOvlpElems();
+   void ComputeOvlpTdofs();
+   void PrintOvlpTdofs();
+
+   // -----------------------------------------------
+   //           Subdomain to Global maps
+   // -----------------------------------------------
+   std::vector<Array<int>> SubdomainGTrueDofs; // Subdomain Tdofs to Global Tdofs
+   std::vector<Array<int>> SubdomainLTrueDofs; // Subdomain Tdofs to Local (on rank) Tdofs
+
+   Array<int> send_count, send_displ;  
+   Array<int> recv_count, recv_displ;  
+   int sbuff_size = 0;
+   int rbuff_size = 0;
+   void SubdomainToGlobalMapsSetup();
+
+   // Testing
+   void TestSubdomainToGlobalMaps();
+   void TestSubdomainToSubdomainMaps();
+
+public:
+   // constructor
+
+   // FiniteElementSpaces of the subdomains
+   Array<FiniteElementSpace *> fes;
+
+   DofMaps(ParFiniteElementSpace *fespace_, ParMeshPartition * part_, bool CompFlag_ = false);
+   ~DofMaps();
+   // Transfering from subdomains SubdomainIds to all their neighbors
+   void TransferToNeighbors(const Array<int> & SubdomainIds, const Array<Vector *> & x, 
+   std::vector<std::vector<Vector * >> & OvlpSol);
+
+   // Prolongation of subdomain solutions to the global solution
+   void SubdomainsToGlobal(const Array<Vector*> & x, Vector & y);
+   // Restriction of global residual to subdomain residuals
+   // bool comp: true for complex valued problems
+   void GlobalToSubdomains(const Vector & y, Array<Vector*> & x);
+};

examples/maxwell-solver-dev/ParDST/ParDST.cpp

@@ -0,0 +1,849 @@
+//Parallel Diagonal Source Transfer Preconditioner
+
+#include "ParDST.hpp"
+
+ParDST::ParDST(ParSesquilinearForm * bf_, Array2D<double> & Pmllength_, 
+         double omega_, Coefficient * Q_,  int nrlayers_ , int nx_, int ny_, int nz_)
+   : Solver(2*bf_->ParFESpace()->GetTrueVSize(), 2*bf_->ParFESpace()->GetTrueVSize()), 
+     bf(bf_), Pmllength(Pmllength_), omega(omega_), 
+     Q(Q_), nrlayers(nrlayers_)
+{
+   nx = nx_; ny = ny_; nz = nz_;
+   Init();
+}
+ParDST::ParDST(ParSesquilinearForm * bf_, Array2D<double> & Pmllength_, 
+         double omega_, VectorCoefficient * VQ_,  int nrlayers_ , int nx_, int ny_, int nz_)
+   : Solver(2*bf_->ParFESpace()->GetTrueVSize(), 2*bf_->ParFESpace()->GetTrueVSize()), 
+     bf(bf_), Pmllength(Pmllength_), omega(omega_), 
+     VQ(VQ_), nrlayers(nrlayers_)
+{
+   nx = nx_; ny = ny_; nz = nz_;
+   Init();
+}
+ParDST::ParDST(ParSesquilinearForm * bf_, Array2D<double> & Pmllength_, 
+         double omega_, MatrixCoefficient * MQ_,  int nrlayers_ , int nx_, int ny_, int nz_)
+   : Solver(2*bf_->ParFESpace()->GetTrueVSize(), 2*bf_->ParFESpace()->GetTrueVSize()), 
+     bf(bf_), Pmllength(Pmllength_), omega(omega_), 
+     MQ(MQ_), nrlayers(nrlayers_)
+{
+   nx = nx_; ny = ny_; nz = nz_;
+   Init();
+}
+
+void ParDST::Init()
+{
+   pfes = bf->ParFESpace();
+   fec = pfes->FEColl();
+ 
+   comm = pfes->GetComm();
+   MPI_Comm_size(comm, &num_procs);
+   MPI_Comm_rank(comm, &myid);
+
+   //1. Indentify problem ... Helmholtz or Maxwell
+   prob_kind = fec->GetContType();
+   if (myid == 0)
+   {
+      cout << " 1. Indentify problem to be solved ... " << endl;
+      if (prob_kind == 0) cout << "    Helmholtz" << endl;
+      if (prob_kind == 1) cout << "    Maxwell" << endl; 
+   }
+
+   //2. Create the parallel mesh partition
+   pmesh = pfes->GetParMesh();
+   dim = pmesh->Dimension();
+   if (myid == 0)
+   {
+      cout << "\n 2. Generating ParMesh partitioning ... " << endl;
+   }
+   ovlpnrlayers = nrlayers+1;
+   part = new ParMeshPartition(pmesh,nx,ny,nz,ovlpnrlayers);
+   nxyz.SetSize(3);
+   nxyz[0] = nx = part->nxyz[0]; 
+   nxyz[1] = ny = part->nxyz[1];  
+   nxyz[2] = nz = part->nxyz[2];
+
+   nrsubdomains = part->nrsubdomains;
+   SubdomainRank = part->subdomain_rank;
+
+   for (int ip = 0; ip<nrsubdomains; ip++)
+   {
+      if (myid == SubdomainRank[ip])
+      {
+         RankSubdomains.Append(ip);
+      }
+   }
+
+   cout << " myid: " << myid 
+        << ", nrsubdomains: " << RankSubdomains.Size() << endl;
+
+   MPI_Barrier(comm);
+
+   if (myid == 0)
+   {
+      cout << "    Done ! " << endl;
+   }
+   //3. Setup info for sweeps
+   if (myid == 0)
+   {
+      cout << "\n 3. Computing sweeps info ..." << endl; 
+   }
+   sweeps = new Sweep(dim);
+   if (myid == 0)
+   {
+      cout << "    Done ! " << endl;
+   }
+   //4. Create LocalToGlobal maps 
+   //   (local GridFunctions/Vector to Global ParGridFunction/Vector) 
+   if (myid == 0)
+   {
+      cout << "\n 4. Computing true dofs maps ..." << endl; 
+   }
+
+   // if (myid == SubdomainRank[0])
+   // {
+   //    cout << "myid = " << myid << endl;
+   //    char vishost[] = "localhost";
+   //    int visport = 19916;
+   //    socketstream mesh_sock1(vishost, visport);
+   //    mesh_sock1.precision(8);
+   //    mesh_sock1 << "mesh\n"
+   //               << *part->subdomain_mesh[0] << "window_title 'Subdomain'" << flush;
+   //    part->subdomain_mesh[0]->Print();
+   
+   // }
+   bool comp = true;
+
+   dmaps = new DofMaps(pfes,part, comp);
+   
+   if (myid == 0)
+   {
+      cout << "    Done ! " << endl;
+   }
+   // 4. Setting up the local problems 
+   if (myid == 0)
+   {
+      cout << "\n 5. Setting up the subdomain problems ..." << endl; 
+   }
+
+   SetupSubdomainProblems();
+   
+   if (myid == 0)
+   {
+      cout << "    Done ! " << endl;
+   }
+   if (myid == 0)
+   {
+      cout << "\n 6. Mark subdomain overlap truedofs ..." << endl; 
+   }
+   MarkSubdomainOverlapDofs(comp);
+   if (myid == 0)
+   {
+      cout << "    Done ! " << endl;
+   }
+}
+
+void ParDST::Mult(const Vector &r, Vector &z) const
+{
+   // Initialize transfered residuals to 0.0;
+   for (int ip=0; ip<nrsubdomains; ip++)
+   {
+      if (myid != SubdomainRank[ip]) continue;
+      for (int i=0;i<sweeps->nsweeps; i++)
+      {
+         *f_transf[ip][i] = 0.0;
+      }
+   }
+
+   // restrict given residual to subdomains
+   dmaps->GlobalToSubdomains(r,f_orig);
+
+   for (int ip=0; ip<nrsubdomains; ip++)
+   {
+      if (myid != SubdomainRank[ip]) continue;
+      Array<int> ijk(3);
+      GetSubdomainijk(ip,nxyz,ijk);
+      Array2D<int> direct(dim,2); direct = 0;
+      for (int d=0;d<dim; d++)
+      {
+         if (ijk[d] > 0) direct[d][0] = 1; 
+         if (ijk[d] < part->nxyz[d]-1) direct[d][1] = 1; 
+      }
+      GetChiRes(*f_orig[ip],ip,direct);
+   }
+
+   z = 0.0; 
+   int nsteps;
+   switch(dim)
+   {
+      case 1: nsteps = nx; break;
+      case 2: nsteps = nx+ny-1; break;
+      default: nsteps = nx+ny+nz-2; break;
+   }
+   int nsweeps = sweeps->nsweeps;
+   // 1. Loop through sweeps
+   if (dim == 3 && nz == 1) { nsweeps = 4; } // x-y partition only;
+   for (int l=0; l<nsweeps; l++)
+   {  
+      // 2. loop through diagonals/steps of each sweep   
+      for (int s = 0; s<nsteps; s++)
+      {
+         Array2D<int> subdomains;
+         GetStepSubdomains(l,s,subdomains);
+         int nsubdomains = subdomains.NumRows();
+
+         // 3. Loop through the subdomains on the diagonal
+         Array<int> subdomain_ids;
+         for (int sb=0; sb < nsubdomains; sb++)
+         {
+            Array<int> ijk(dim); ijk = 0;
+            for (int d=0; d<dim; d++) ijk[d] = subdomains[sb][d]; 
+            int ip = GetSubdomainId(nxyz,ijk);
+            subdomain_ids.Append(ip);
+            if (myid != SubdomainRank[ip]) continue;
+
+            int n = dmaps->fes[ip]->GetTrueVSize();
+            Vector res_local(2*n); res_local = 0.0;
+
+            if (l==0)  { res_local += *f_orig[ip]; }
+            res_local += *f_transf[ip][l];
+            if (res_local.Norml2() < 1e-12) 
+            {
+               *subdomain_sol[ip] = 0.0;
+               continue;
+            }
+            PmlMatInv[ip]->Mult(res_local, *subdomain_sol[ip]);
+         }
+         // 4. Transfer solutions to neighbors so that the subdomain
+         // residuals are updated
+         TransferSources(l,subdomain_ids);
+      }
+      // 5. Update the global solution 
+      dmaps->SubdomainsToGlobal(subdomain_sol,z);
+   }
+}
+
+void ParDST::SetupSubdomainProblems()
+{
+   sqf.SetSize(nrsubdomains);
+   Optr.SetSize(nrsubdomains);
+   PmlMat.SetSize(nrsubdomains);
+   PmlMatInv.SetSize(nrsubdomains);
+   f_orig.SetSize(nrsubdomains);
+   f_transf.SetSize(nrsubdomains);
+   subdomain_sol.SetSize(nrsubdomains);
+   for (int ip=0; ip<nrsubdomains; ip++)
+   {
+      sqf[ip] = nullptr;
+      f_orig[ip] = nullptr;
+      subdomain_sol[ip] = nullptr;
+      PmlMat[ip] = nullptr;
+      PmlMatInv[ip] = nullptr;
+      Optr[ip] = nullptr;
+
+      if (myid != SubdomainRank[ip]) continue;
+      subdomain_sol[ip] = new Vector(2*dmaps->fes[ip]->GetTrueVSize());
+      if (prob_kind == 0)
+      {
+         SetHelmholtzPmlSystemMatrix(ip);
+      }
+      else if (prob_kind == 1)
+      {
+         SetMaxwellPmlSystemMatrix(ip);
+      }
+      PmlMat[ip] = Optr[ip]->As<ComplexSparseMatrix>();
+
+      PmlMatInv[ip] = new ComplexUMFPackSolver;
+      PmlMatInv[ip]->Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+      PmlMatInv[ip]->SetOperator(*PmlMat[ip]);
+
+      // HYPRE_Int rowstarts[2]; rowstarts[0] = 0;
+      // rowstarts[1] = dmaps->fes[ip]->GetTrueVSize();
+      // HypreParMatrix * HypreMat_r =
+      //    new HypreParMatrix(MPI_COMM_SELF,rowstarts[1],rowstarts,
+      //                       &(PmlMat[ip]->real()));
+      // HypreParMatrix * HypreMat_i =
+      //    new HypreParMatrix(MPI_COMM_SELF,rowstarts[1],rowstarts,
+      //                       &(PmlMat[ip]->imag()));
+      // ComplexHypreParMatrix * HypreMat =
+      //    new ComplexHypreParMatrix(HypreMat_r,HypreMat_i,true,true);
+      // PmlMatInv[ip] = new ComplexMUMPSSolver;
+      // PmlMatInv[ip]->SetOperator(*HypreMat);
+      // delete HypreMat;
+      int ndofs = dmaps->fes[ip]->GetTrueVSize();
+      f_transf[ip].SetSize(sweeps->nsweeps);
+      for (int i=0;i<sweeps->nsweeps; i++)
+      {
+         f_transf[ip][i] = new Vector(2*ndofs);
+      }
+   }
+}
+
+
+
+void ParDST::SetHelmholtzPmlSystemMatrix(int ip)
+{
+   MFEM_VERIFY(part->subdomain_mesh[ip], "Null mesh pointer");
+   Mesh * mesh = part->subdomain_mesh[ip];
+   double h = part->MeshSize;
+   Array2D<double> length(dim,2);
+   length = h*(nrlayers);
+
+   Array<int> ijk;
+   GetSubdomainijk(ip,nxyz,ijk);
+   int i = ijk[0];
+   int j = ijk[1];
+   int k = ijk[2];
+
+   if (i == 0 )    length[0][0] = Pmllength[0][0];
+   if (i == nx-1 ) length[0][1] = Pmllength[0][1];
+   if (dim  > 1)
+   {
+      if (j == 0 )    length[1][0] = Pmllength[1][0];
+      if (j == ny-1 ) length[1][1] = Pmllength[1][1];
+   }
+   if (dim  == 3)
+   {
+      if (k == 0 )    length[2][0] = Pmllength[2][0];
+      if (k == nz-1 ) length[2][1] = Pmllength[2][1];
+   }
+   
+   CartesianPML pml(mesh, length);
+   pml.SetOmega(omega);
+
+   Array <int> ess_tdof_list;
+   if (mesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      dmaps->fes[ip]->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient sigma(-pow(omega, 2));
+   PmlMatrixCoefficient c1_re(dim,pml_detJ_JT_J_inv_Re,&pml);
+   PmlMatrixCoefficient c1_im(dim,pml_detJ_JT_J_inv_Im,&pml);
+   PmlCoefficient detJ_re(pml_detJ_Re,&pml);
+   PmlCoefficient detJ_im(pml_detJ_Im,&pml);
+   ProductCoefficient c2_re0(sigma, detJ_re);
+   ProductCoefficient c2_im0(sigma, detJ_im);
+   ProductCoefficient c2_re(c2_re0, *Q);
+   ProductCoefficient c2_im(c2_im0, *Q);
+   sqf[ip] = new SesquilinearForm (dmaps->fes[ip],bf->GetConvention());
+
+   sqf[ip]->AddDomainIntegrator(new DiffusionIntegrator(c1_re),
+                         new DiffusionIntegrator(c1_im));
+   sqf[ip]->AddDomainIntegrator(new MassIntegrator(c2_re),
+                         new MassIntegrator(c2_im));
+   sqf[ip]->Assemble(0);
+
+   Optr[ip] = new OperatorPtr;
+   sqf[ip]->FormSystemMatrix(ess_tdof_list,*Optr[ip]);
+}
+
+void ParDST::SetMaxwellPmlSystemMatrix(int ip)
+{
+   MFEM_VERIFY(part->subdomain_mesh[ip], "Null mesh pointer");
+   Mesh * mesh = part->subdomain_mesh[ip];
+   double h = part->MeshSize;
+   Array2D<double> length(dim,2);
+   length = h*(nrlayers);
+
+   Array<int> ijk;
+   GetSubdomainijk(ip,nxyz,ijk);
+   int i = ijk[0];
+   int j = ijk[1];
+   int k = ijk[2];
+
+   if (i == 0 )    length[0][0] = Pmllength[0][0];
+   if (i == nx-1 ) length[0][1] = Pmllength[0][1];
+   if (dim  > 1)
+   {
+      if (j == 0 )    length[1][0] = Pmllength[1][0];
+      if (j == ny-1 ) length[1][1] = Pmllength[1][1];
+   }
+   if (dim  == 3)
+   {
+      if (k == 0 )    length[2][0] = Pmllength[2][0];
+      if (k == nz-1 ) length[2][1] = Pmllength[2][1];
+   }
+   
+   CartesianPML pml(mesh, length);
+   pml.SetOmega(omega);
+   Array <int> ess_tdof_list;
+   if (mesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      dmaps->fes[ip]->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   ConstantCoefficient omeg(-pow(omega, 2));
+   int cdim = (dim == 2) ? 1 : dim;
+
+   PmlMatrixCoefficient pml_c1_Re(cdim,detJ_inv_JT_J_Re, &pml);
+   PmlMatrixCoefficient pml_c1_Im(cdim,detJ_inv_JT_J_Im, &pml);
+
+   PmlMatrixCoefficient pml_c2_Re(dim, detJ_JT_J_inv_Re,&pml);
+   PmlMatrixCoefficient pml_c2_Im(dim, detJ_JT_J_inv_Im,&pml);
+   ScalarMatrixProductCoefficient c2_Re0(omeg,pml_c2_Re);
+   ScalarMatrixProductCoefficient c2_Im0(omeg,pml_c2_Im);
+
+   MatrixCoefficient * c2_Re=nullptr;
+   MatrixCoefficient * c2_Im=nullptr;
+
+   if (Q)
+   {
+      c2_Re = new ScalarMatrixProductCoefficient(*Q,c2_Re0);
+      c2_Im = new ScalarMatrixProductCoefficient(*Q,c2_Im0);
+   }
+   else if (VQ)
+   {
+      MFEM_ABORT("Vector Coeffiecient not supported ");
+   }
+   else if (MQ)
+   {
+      c2_Re = new MatrixMatrixProductCoefficient(c2_Re0,*MQ);
+      c2_Im = new MatrixMatrixProductCoefficient(c2_Im0,*MQ);
+   }
+   
+   sqf[ip] = new SesquilinearForm(dmaps->fes[ip],bf->GetConvention());
+
+   sqf[ip]->AddDomainIntegrator(new CurlCurlIntegrator(pml_c1_Re),
+                         new CurlCurlIntegrator(pml_c1_Im));
+   sqf[ip]->AddDomainIntegrator(new VectorFEMassIntegrator(*c2_Re),
+                         new VectorFEMassIntegrator(*c2_Im));
+   sqf[ip]->Assemble(0);
+
+   Optr[ip] = new OperatorPtr;
+   sqf[ip]->FormSystemMatrix(ess_tdof_list,*Optr[ip]);
+   delete c2_Re;
+   delete c2_Im;
+}
+
+
+void ParDST::MarkSubdomainOverlapDofs(const bool comp)
+{
+   // First mark the elements
+         // cout<< "Compute Overlap Elements (in each possible direction) " << endl;
+      // Lists of elements
+      // x,y,z = +/- 1 ovlp 
+   NovlpElems.resize(nrsubdomains);
+
+   for (int ip = 0; ip<nrsubdomains; ip++)
+   {
+      if (myid != SubdomainRank[ip]) continue;
+      Array<int> ijk;
+      GetSubdomainijk(ip,nxyz,ijk);
+
+      Mesh * mesh = dmaps->fes[ip]->GetMesh();
+      NovlpElems[ip].resize(2*dim);
+
+      Vector pmin, pmax;
+      mesh->GetBoundingBox(pmin,pmax);
+      double h = part->MeshSize;
+      // Loop through elements
+      for (int iel=0; iel<mesh->GetNE(); iel++)
+      {
+         // Get element center
+         Vector center(dim);
+         int geom = mesh->GetElementBaseGeometry(iel);
+         ElementTransformation * tr = mesh->GetElementTransformation(iel);
+         tr->Transform(Geometries.GetCenter(geom),center);
+
+         // Assign elements to the appropriate lists
+         for (int d=0;d<dim; d++)
+         {
+            if (ijk[d]>0)
+            {
+               if (center[d] >= pmin[d]+h*ovlpnrlayers) 
+               {
+                  NovlpElems[ip][d].Append(iel);
+               }
+            }
+            else
+            {
+               NovlpElems[ip][d].Append(iel);
+            }
+               
+            if (ijk[d]<nxyz[d]-1)
+            {
+               if (center[d] <= pmax[d]-h*ovlpnrlayers) 
+               {
+                  NovlpElems[ip][dim+d].Append(iel);
+               }
+            }
+            else
+            {
+               NovlpElems[ip][dim+d].Append(iel);
+            }
+         }
+      }
+   }
+
+      // mark dofs
+   NovlpDofs.resize(nrsubdomains);
+   int mm = (comp) ? 2 : 1; // complex or real valued
+   for (int ip = 0; ip<nrsubdomains; ip++)
+   {
+      if (myid != SubdomainRank[ip]) continue;
+      FiniteElementSpace * fes = dmaps->fes[ip];
+      // Loop through the marked elements
+      NovlpDofs[ip].resize(2*dim);
+      int n = fes->GetTrueVSize();
+      Array<int> marker(n);
+      for (int d=0;d<2*dim; d++)
+      {
+         marker = 0;
+         int m = 0;
+         int melems = NovlpElems[ip][d].Size(); 
+         for (int iel=0; iel<melems; iel++)
+         {
+            Array<int> ElemDofs;
+            int el =  NovlpElems[ip][d][iel];
+            fes->GetElementDofs(el,ElemDofs);
+            int ndof = ElemDofs.Size();
+            for (int i = 0; i<ndof; ++i)
+            {
+               int eldof = ElemDofs[i];
+               int tdof = (eldof >= 0) ? eldof : abs(eldof) - 1; 
+               if (marker[tdof] == 1) continue;
+               marker[tdof] = 1;
+               m++;
+            }
+         }
+         int k = mm*(n-m);
+         NovlpDofs[ip][d].SetSize(k);
+         int l = 0;
+         for (int i = 0; i<n; i++)
+         {
+            if (marker[i]==0) 
+            {
+               NovlpDofs[ip][d][l] = i;  // real dofs
+               if (comp) 
+               {
+                  NovlpDofs[ip][d][l+k/2] = i+fes->GetTrueVSize();
+               }
+               l++;
+            }
+         }
+      }
+   }
+}
+
+void ParDST::GetChiRes(Vector & res, int ip, Array2D<int> direct) const
+{
+   for (int d=0; d<dim; d++)
+   {
+      // negative direction
+      if (direct[d][0]==1) res.SetSubVector(NovlpDofs[ip][d],0.0);
+      // possitive direction
+      if (direct[d][1]==1) res.SetSubVector(NovlpDofs[ip][d+dim],0.0);
+   }
+}
+
+
+
+void ParDST::PlotLocal(Vector & sol, socketstream & sol_sock, int ip) const
+{
+   FiniteElementSpace * fes = dmaps->fes[ip];
+   Mesh * mesh = fes->GetMesh();
+   GridFunction gf(fes);
+   double * data = sol.GetData();
+   gf.SetData(data);
+   
+   string keys;
+   keys = "keys mrRljc\n";
+   sol_sock << "solution\n" << *mesh << gf << keys << flush;
+}
+
+
+void ParDST::GetStepSubdomains(const int sweep, const int step, Array2D<int> & subdomains) const
+{
+   Array<int> aux;
+   switch(dim)
+   {
+      case 2: 
+      for (int i=nx-1;i>=0; i--)
+      {  
+         int j;
+         switch (sweep)
+         {
+            case 0:  j = step-i;         break;
+            case 1:  j = step-nx+i+1;    break;
+            case 2:  j = nx+i-step-1;    break;
+            default: j = nx+ny-i-step-2; break;
+         }
+         if (j<0 || j>=ny) continue;
+         aux.Append(i); aux.Append(j);
+      }
+      break; 
+      default:
+      for (int i=nx-1;i>=0; i--)
+      {  
+         for (int j=ny-1;j>=0; j--)
+         {
+            int k;
+            switch (sweep)
+            {
+               case 0:  k = step-i-j;            break;
+               case 1:  k = step-nx+i+1-j;       break;
+               case 2:  k = step-ny+j+1-i;       break;
+               case 3:  k = step-nx-ny+i+j+2;    break;
+               case 4:  k = i+j+nz-1-step;       break;
+               case 5:  k = nx+nz-i+j-step-2;    break;
+               case 6:  k = ny+nz+i-j-step-2;    break;
+               default: k = nx+ny+nz-i-j-step-3; break;
+            }
+            if (k<0 || k>=nz) continue;   
+            aux.Append(i); aux.Append(j); aux.Append(k); 
+         }
+      }
+      break;
+   }
+
+   int nrows = aux.Size()/dim;
+   int ncols = dim;
+
+   subdomains.SetSize(nrows,ncols);
+   for (int r=0;r<nrows; r++)
+   {
+      for (int c=0; c<ncols; c++)
+      {
+         int k = r*ncols + c;
+         subdomains[r][c] = aux[k];
+      }
+   }
+}
+
+
+void ParDST::TransferSources(int sweep, const Array<int> & subdomain_ids) const
+{
+   OvlpSol.resize(nrsubdomains);
+   int nrneighbors = pow(3,dim);
+   for (int ip = 0; ip<nrsubdomains; ip++)
+   {
+      if (myid == SubdomainRank[ip])
+      {
+         OvlpSol[ip].resize(nrneighbors);
+      }
+   }
+   int m = subdomain_ids.Size();
+   Array<Vector *> x(m);
+   for (int i = 0; i<m; i++)
+   {
+      x[i] = nullptr;
+      int ip = subdomain_ids[i];
+      if (myid != SubdomainRank[ip]) continue;
+      x[i] = new Vector(subdomain_sol[ip]->GetData(),subdomain_sol[ip]->Size());
+   }
+   dmaps->TransferToNeighbors(subdomain_ids,x,OvlpSol);
+   for (int i = 0; i<m; i++)
+   {
+      delete x[i]; x[i] = nullptr;
+   }
+   // Update residuals
+//   Find all neighbors of patch ip0
+   for (int is = 0; is<m; is++)
+   {
+      int ip0 = subdomain_ids[is];
+      Array<int> ijk;
+      Array<int> ijk1(3);
+      GetSubdomainijk(ip0,nxyz,ijk);
+      Array<int> directions(3);   
+      for (int i=-1; i<2; i++)
+      {
+         int i1 = ijk[0] + i;
+         if (i1 <0 || i1>=nx) continue;
+         directions[0] = i;
+         ijk1[0] = i1;
+         for (int j=-1; j<2; j++)
+         {
+            int j1 = ijk[1] + j;
+            if (j1 <0 || j1>=ny) continue;
+            directions[1] = j;
+            ijk1[1] = j1;
+            int kbeg = (dim == 2) ? 0 : -1;
+            int kend = (dim == 2) ? 1 :  2;
+            for (int k=kbeg; k<kend; k++)
+            {
+               int k1 = ijk[2] + k;
+               if (k1 <0 || k1>=nz) continue;
+               directions[2] = (dim == 3) ? k : -1 ;
+               if (i==0 && j==0 && k==0) continue;
+
+               int l = GetSweepToTransfer(sweep,directions);
+               if (l == -1) continue;
+               ijk1[2] = k1;
+               int ip1 = GetSubdomainId(nxyz,ijk1);
+
+               if (myid != SubdomainRank[ip1]) continue;
+               Array<int>directions1(3); directions1 = -1;
+               for (int i = 0; i<dim; i++) directions1[i] = -directions[i];
+               int dir = GetDirectionId(directions1);
+               int n = dmaps->fes[ip1]->GetTrueVSize();
+               Vector res(2*n);
+               PmlMat[ip1]->Mult(*OvlpSol[ip1][dir],res);
+
+               Array2D<int> direct(dim,2); direct = 0;
+               for (int d = 0; d<dim; d++)
+               {
+                  if (directions[d]==1)  direct[d][0] = 1;
+                  if (directions[d]==-1) direct[d][1] = 1;
+               }   
+               GetChiRes(res,ip1,direct);
+               *f_transf[ip1][l] -= res;
+            }
+         }  
+      }
+   }
+
+   for (int ip = 0; ip<nrsubdomains; ip++)
+   {
+      if (myid == SubdomainRank[ip])
+      {
+         for (int i = 0; i<nrneighbors; i++)
+         {
+            if (OvlpSol[ip][i])
+            {
+               delete OvlpSol[ip][i];
+            }
+         }
+         OvlpSol[ip].clear();
+      }
+   }
+}
+
+int ParDST::GetSweepToTransfer(const int s, Array<int> directions) const
+{
+   int l1=-1;
+   int nsweeps = sweeps->nsweeps;
+   Array<int> sweep0;   
+   sweeps->GetSweep(s,sweep0);
+   switch (dim)
+   {
+      case 2:
+      for (int l=s; l<nsweeps; l++)
+      {
+         // Rule 1: the transfer source direction has to be similar with 
+         // the sweep direction
+         Array<int> sweep1;   
+         sweeps->GetSweep(l,sweep1);
+         int ddot = 0;
+         for (int d=0; d<dim; d++) ddot+= sweep1[d] * directions[d];
+         if (ddot <= 0) continue;
+
+         // Rule 2: The horizontal or vertical transfer source cannot be used
+         // Case of horizontal or vertical transfer source 
+         // (it can't be both 0 cause it's skipped)
+         if (directions[0]==0 || directions[1] == 0) 
+         {
+            if (sweep0[0] == -sweep1[0] && sweep0[1] == -sweep1[1]) continue;
+         }
+         l1 = l;
+         break;
+      }
+      break;
+      default:
+      for (int l=s; l<nsweeps; l++)
+      {
+         // Rule 1: (similar directions) the transfer source direction has to be similar with 
+         // the sweep direction
+         Array<int> sweep1;   
+         sweeps->GetSweep(l,sweep1);
+         int ddot = 0;
+         bool similar = true;
+         for (int d=0; d<dim; d++) 
+         {
+            if (sweep1[d] * directions[d] < 0) similar = false;
+            ddot+= sweep1[d] * directions[d];
+         }
+         if (!similar || ddot<=0) continue; // not similar
+
+         // Rule 2: (oposite directions) the transfer source direction has to be similar with 
+         // the sweep direction
+         // 
+         // check any of the projections onto the planes
+         // (xy, xz, yz)
+
+         if ( (directions[0]==0 && directions[1] != 0) || 
+              (directions[0]!=0 && directions[1] == 0) ||
+              (directions[0]==0 && directions[2] != 0) || 
+              (directions[0]!=0 && directions[2] == 0) ||
+              (directions[2]==0 && directions[1] != 0) || 
+              (directions[2]!=0 && directions[1] == 0) ) 
+         {
+            if (sweep0[0] == -sweep1[0] && 
+                sweep0[1] == -sweep1[1] && 
+                sweep0[2] == -sweep1[2]) continue;
+         }
+         l1 = l;
+         break;
+      }
+      break;
+   }
+
+   return l1;   
+}
+
+void ParDST::CorrectOrientation(int ip,Vector &x) const
+{
+   FiniteElementSpace * fespace = dmaps->fes[ip];
+   Mesh * mesh = fespace->GetMesh();
+   int nrelems = mesh->GetNE();
+   // GridFunction test; 
+   // test.SetFromTrueDofs(x)
+   Array<int> signs(fespace->GetTrueVSize()); signs = 0;
+   for (int iel=0; iel<nrelems; iel++)
+   {
+      Array<int> ElemDofs;
+      fespace->GetElementDofs(iel,ElemDofs);
+      int ndofs = ElemDofs.Size();
+      ElemDofs.Print();
+      for (int i = 0; i< ndofs; i++)
+      {
+         int pdof_ = ElemDofs[i];
+         if (pdof_ < 0) 
+         {
+            signs[abs(pdof_)-1] += 1.0 ; 
+         }
+         else
+         {
+            signs[pdof_] -= 1.0 ; 
+         }
+      }
+   }
+
+   cout << "signs = " ; signs.Print();
+   for (int i = 0; i<fespace->GetTrueVSize(); i++)
+   {
+      if (signs[i]<0) 
+      {
+         x(i) *= -1.0;
+         x(i+fespace->GetTrueVSize()) *= -1.0;
+      }
+   }
+}
+
+
+ParDST::~ParDST()
+{
+
+   for (int ip=0; ip<nrsubdomains; ip++)
+   {
+      delete Optr[ip];
+      delete subdomain_sol[ip];
+      delete PmlMatInv[ip];
+      delete sqf[ip];
+      if (myid != SubdomainRank[ip]) continue;
+      for (int i=0;i<sweeps->nsweeps; i++)
+      {
+         delete f_transf[ip][i];
+      }
+      delete f_orig[ip];
+   }
+   f_orig.DeleteAll();
+   delete dmaps;
+   delete sweeps;
+   delete part;
+
+}

examples/maxwell-solver-dev/ParDST/ParDST.hpp

@@ -0,0 +1,70 @@
+#pragma once
+#include "../common/Utilities.hpp"
+#include "../common/PML.hpp"
+#include "DofMapsDST.hpp"
+using namespace std;
+using namespace mfem;
+
+
+class ParDST : public Solver//
+{
+private:
+   MPI_Comm comm = MPI_COMM_WORLD;
+   int num_procs, myid;
+   // Constructor inputs
+   int prob_kind;
+   ParSesquilinearForm *bf=nullptr;
+   ParFiniteElementSpace * pfes = nullptr;
+   ParMesh * pmesh = nullptr;
+   ParMeshPartition * part = nullptr;
+   Array<int> SubdomainRank;
+   Array<int> RankSubdomains;
+   const FiniteElementCollection * fec = nullptr;
+   Array2D<double> Pmllength;
+   int dim = 2;
+   double omega = 0.5;
+   Coefficient * Q=nullptr;
+   VectorCoefficient * VQ=nullptr;
+   MatrixCoefficient * MQ=nullptr;
+   int nrlayers;
+   int ovlpnrlayers;
+   int nrsubdomains = 0;
+   int nx,ny,nz;
+   Array<int> nxyz;
+   Sweep * sweeps = nullptr;
+   DofMaps * dmaps = nullptr;
+   Array< SesquilinearForm * > sqf;
+   Array< OperatorPtr * > Optr;
+   Array<ComplexSparseMatrix *> PmlMat;
+   Array<ComplexUMFPackSolver *> PmlMatInv;
+//    Array<ComplexMUMPSSolver *> PmlMatInv;
+   mutable Array<Vector *> f_orig;
+   mutable Array<Array<Vector * >> f_transf;
+   mutable Array<Vector * > subdomain_sol;
+   mutable std::vector<std::vector<Vector * >> OvlpSol;
+   void SetupSubdomainProblems();
+   std::vector<std::vector<Array<int>>> NovlpElems;
+   std::vector<std::vector<Array<int>>> NovlpDofs;
+   void MarkSubdomainOverlapDofs(const bool comp = false);
+   void SetHelmholtzPmlSystemMatrix(int ip);   
+   void SetMaxwellPmlSystemMatrix(int ip);
+   void GetChiRes(Vector & res, int ip, Array2D<int> direct) const;
+   void PlotLocal(Vector & sol, socketstream & sol_sock, int ip) const;
+   void GetStepSubdomains(const int sweep, const int step, Array2D<int> & subdomains) const;
+   void TransferSources(int sweep, const Array<int> & subdomain_ids) const;
+   int GetSweepToTransfer(const int s, Array<int> directions) const;
+   void CorrectOrientation(int ip, Vector & x) const;
+   void Init();
+public:
+   ParDST(ParSesquilinearForm * bf_, Array2D<double> & Pmllength_, 
+       double omega_, Coefficient * Q_, int nrlayers_, int nx_=2, int ny_=2, int nz_=2);
+   ParDST(ParSesquilinearForm * bf_, Array2D<double> & Pmllength_, 
+       double omega_, VectorCoefficient * VQ_, int nrlayers_, int nx_=2, int ny_=2, int nz_=2);
+   ParDST(ParSesquilinearForm * bf_, Array2D<double> & Pmllength_, 
+       double omega_, MatrixCoefficient * MQ_, int nrlayers_, int nx_=2, int ny_=2, int nz_=2);              
+   virtual void SetOperator(const Operator &op) {}
+   virtual void Mult(const Vector &r, Vector &z) const;
+   virtual ~ParDST();
+};
+
+

examples/maxwell-solver-dev/common/MeshPartition.hpp

@@ -0,0 +1,151 @@
+#pragma once
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+using namespace std;
+using namespace mfem;
+
+
+struct UniqueIndexGenerator
+{
+   int counter = 0;
+   std::unordered_map<int,int> idx;
+   int Get(int i)
+   {
+      std::unordered_map<int,int>::iterator f = idx.find(i);
+      if (f == idx.end())
+      {
+         idx[i] = counter;
+         return counter++;
+      }
+      else
+      {
+         return (*f).second;
+      }
+   }
+   void Reset()
+   {
+      counter = 0;
+      idx.clear();
+   }
+};
+
+
+double GetUniformMeshElementSize(Mesh * mesh);
+Mesh * ExtendMesh(Mesh * mesh, const Array<int> & directions);
+
+class CartesianMeshPartition 
+{
+private:
+   Mesh *mesh=nullptr;
+public:
+   int nrpatch;
+   int nxyz[3];
+   double MeshSize;
+   std::vector<Array<int>> element_map;
+   Array3D<int>subdomains;
+   // constructor
+   CartesianMeshPartition(Mesh * mesh_,int & nx, int & ny, int & nz);
+   ~CartesianMeshPartition() {};
+};
+
+class OverlappingCartesianMeshPartition 
+{
+private:
+   Mesh *mesh=nullptr;
+public:
+   int nrpatch;
+   double MeshSize;
+   int nxyz[3];
+   std::vector<Array<int>> element_map;
+   Array3D<int> subdomains;
+   // constructor
+   OverlappingCartesianMeshPartition(Mesh * mesh_,int & nx, int & ny, int & nz);
+   OverlappingCartesianMeshPartition(Mesh * mesh_,int & nx, int & ny, int & nz, int ovlp_nlayers);
+   ~OverlappingCartesianMeshPartition() {};
+};
+
+class STPOverlappingCartesianMeshPartition // Special layered partition for STP
+{
+private:
+   Mesh *mesh=nullptr;
+public:
+   int nrpatch;
+   int nx, ny, nz;
+   std::vector<Array<int>> element_map;
+   // constructor
+   STPOverlappingCartesianMeshPartition(Mesh * mesh_);
+   ~STPOverlappingCartesianMeshPartition() {};
+};
+
+class MeshPartition
+{
+private:
+   Mesh *mesh=nullptr;
+   void AddElementToMesh(Mesh * mesh,mfem::Element::Type elem_type,int * ind);
+   void GetNumVertices(int type, mfem::Element::Type & elem_type, int & nrvert);
+   void PrintElementMap();
+public:
+   int nrpatch;
+   double MeshSize;
+   std::vector<Array<int>> element_map;
+   Array3D<int> subdomains;
+   Array<Mesh *> patch_mesh;
+   int partition_kind;
+   int nxyz[3];
+   // constructor
+   MeshPartition(Mesh * mesh_, int part, int mx=1, int my=1, int mz=1, int ovl_nlayers=0);
+   ~MeshPartition();
+};
+
+void SaveMeshPartition(Array<Mesh * > meshes, 
+                       string mfilename="output/mesh.",
+                       string sfilename="output/sol.");
+
+
+#ifdef MFEM_USE_MPI
+
+class CartesianParMeshPartition 
+{
+private:
+   ParMesh *pmesh=nullptr;
+public:
+   int nrsubdomains;
+   int nxyz[3];
+   double MeshSize;
+   std::vector<Array<int>> local_element_map;
+   Array<int> subdomain_rank;
+   Array3D<int>subdomains;
+   // constructor
+   CartesianParMeshPartition(ParMesh * pmesh_,int & nx, int & ny, int & nz, 
+                             int ovlp_nlayers);
+   ~CartesianParMeshPartition() {};
+};
+
+class ParMeshPartition
+{
+private:
+   MPI_Comm comm;
+   ParMesh *pmesh=nullptr;
+   void AddElementToMesh(Mesh * mesh,mfem::Element::Type elem_type,int * ind);
+   void GetNumVertices(int type, mfem::Element::Type & elem_type, int & nrvert);
+   void PrintElementMap();
+public:
+   int nrsubdomains;
+   int OvlpNlayers;
+   int myelem_offset = 0;
+   double MeshSize;
+   std::vector<Array<int>> element_map;
+   std::vector<Array<int>> local_element_map;
+   Array3D<int> subdomains;
+   Array<Mesh *> subdomain_mesh;
+   Array<int> subdomain_rank;
+   int partition_kind;
+   int nxyz[3];
+   // constructor
+   ParMeshPartition(ParMesh * pmesh_, int mx=1, int my=1, int mz=1, int ovl_nlayers=0);
+   void SaveMeshPartition();
+   ~ParMeshPartition();
+};
+
+#endif

examples/maxwell-solver-dev/common/PML.cpp

@@ -0,0 +1,570 @@
+#include "PML.hpp"
+
+CartesianPML::CartesianPML(Mesh *mesh_, Array2D<double> length_)
+   : mesh(mesh_), length(length_)
+{
+   dim = mesh->Dimension();
+   SetBoundaries();
+}
+
+void CartesianPML::SetBoundaries()
+{
+   comp_dom_bdr.SetSize(dim, 2);
+   dom_bdr.SetSize(dim, 2);
+   // initialize
+   for (int i = 0; i < dim; i++)
+   {
+      dom_bdr(i, 0) = infinity();
+      dom_bdr(i, 1) = -infinity();
+   }
+
+   for (int i = 0; i < mesh->GetNBE(); i++)
+   {
+      Array<int> bdr_vertices;
+      mesh->GetBdrElementVertices(i, bdr_vertices);
+      for (int j = 0; j < bdr_vertices.Size(); j++)
+      {
+         for (int k = 0; k < dim; k++)
+         {
+            dom_bdr(k, 0) = min(dom_bdr(k, 0), mesh->GetVertex(bdr_vertices[j])[k]);
+            dom_bdr(k, 1) = max(dom_bdr(k, 1), mesh->GetVertex(bdr_vertices[j])[k]);
+         }
+      }
+   }
+
+#ifdef MFEM_USE_MPI
+   ParMesh * pmesh = dynamic_cast<ParMesh *>(mesh);
+   if (pmesh)
+   {
+      for (int d=0; d<dim; d++)
+      {
+         MPI_Allreduce(MPI_IN_PLACE,&dom_bdr(d,0),1,MPI_DOUBLE,MPI_MIN,pmesh->GetComm());
+         MPI_Allreduce(MPI_IN_PLACE,&dom_bdr(d,1),1,MPI_DOUBLE,MPI_MAX,pmesh->GetComm());
+      }
+   }
+#endif
+
+   for (int i = 0; i < dim; i++)
+   {
+      comp_dom_bdr(i, 0) = dom_bdr(i, 0) + length(i, 0);
+      comp_dom_bdr(i, 1) = dom_bdr(i, 1) - length(i, 1);
+   }
+}
+
+void CartesianPML::SetAttributes(Mesh *mesh_)
+{
+   int nrelem = mesh_->GetNE();
+   elems.SetSize(nrelem);
+
+   for (int i = 0; i < nrelem; ++i)
+   {
+      elems[i] = 1;
+      bool in_pml = false;
+      Element *el = mesh_->GetElement(i);
+      Array<int> vertices;
+      // Initialize Attribute
+      el->SetAttribute(1);
+      el->GetVertices(vertices);
+      int nrvert = vertices.Size();
+      // Check if any vertex is in the pml
+      for (int iv = 0; iv < nrvert; ++iv)
+      {
+         int vert_idx = vertices[iv];
+         double *coords = mesh_->GetVertex(vert_idx);
+         for (int comp = 0; comp < dim; ++comp)
+         {
+            if (coords[comp] > comp_dom_bdr(comp, 1) ||
+                coords[comp] < comp_dom_bdr(comp, 0))
+            {
+               in_pml = true;
+               break;
+            }
+         }
+      }
+      if (in_pml)
+      {
+         elems[i] = 0;
+         el->SetAttribute(2);
+      }
+   }
+   mesh_->SetAttributes();
+}
+
+void CartesianPML::StretchFunction(const Vector &x,
+                                   vector<complex<double>> &dxs, double omega)
+{
+   complex<double> zi = complex<double>(0., 1.);
+
+   double n = 2.0;
+   double c = 10.0;
+   // double c = log(omega);
+   double coeff;
+   // Stretch in each direction independently
+   for (int i = 0; i < dim; ++i)
+   {
+      dxs[i] = 1.0;
+      if (x(i) >= comp_dom_bdr(i, 1))
+      {
+         coeff = n * c / omega / pow(length(i, 1), n);
+         dxs[i] = 1.0 + zi * coeff * abs(pow(x(i) - comp_dom_bdr(i, 1), n - 1.0));
+      }
+      if (x(i) <= comp_dom_bdr(i, 0))
+      {
+         coeff = n * c / omega / pow(length(i, 0), n);
+         dxs[i] = 1.0 + zi * coeff * abs(pow(x(i) - comp_dom_bdr(i, 0), n - 1.0));
+      }
+   }
+}
+
+ToroidPML::ToroidPML(Mesh *mesh_)
+   : mesh(mesh_)
+{
+   dim = mesh->Dimension();
+   zlim.SetSize(2);
+   rlim.SetSize(2);
+   alim.SetSize(2);
+   zpml_thickness.SetSize(2);
+   rpml_thickness.SetSize(2);
+   apml_thickness.SetSize(2);
+   SetBoundaries();
+}
+
+void ToroidPML::SetBoundaries()
+{
+   mesh->EnsureNodes();
+   int nrnodes = mesh->GetNodalFESpace()->GetTrueVSize()/dim;
+   double zmin = infinity();
+   double zmax = -infinity();
+   double rmin = infinity();
+   double rmax = -infinity();
+   double amin = infinity(); // in degrees
+   double amax = -infinity(); // in degrees
+   for (int i = 0; i<nrnodes; i++)
+   {
+      Vector coord(dim);
+      mesh->GetNode(i,coord);
+      for (int d = 0; d<dim; d++)
+      {
+         if (abs(coord[d])<1e-13) coord[d] = 0.0;
+      }
+      // Find r and a for this point
+      double x = coord[0];
+      double y = coord[1];
+      double z = 0.0;
+      if (dim == 3) z = coord[2];
+      double a = GetAngle(x,y);
+      double r = sqrt(x*x + y*y);
+      
+      zmin = min(zmin,z);
+      zmax = max(zmax,z);
+      rmin = min(rmin,r);
+      rmax = max(rmax,r);
+      amin = min(amin,a);
+      amax = max(amax,a);
+   }
+
+   zlim[0] = zmin;
+   zlim[1] = zmax;
+   rlim[0] = rmin;
+   rlim[1] = rmax;
+   alim[0] = amin;
+   alim[1] = amax;
+
+
+#ifdef MFEM_USE_MPI
+   ParMesh * pmesh = dynamic_cast<ParMesh *>(mesh);
+   if (pmesh)
+   {
+      MPI_Allreduce(MPI_IN_PLACE,&zlim[0],1,MPI_DOUBLE,MPI_MIN,pmesh->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE,&zlim[1],1,MPI_DOUBLE,MPI_MAX,pmesh->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE,&rlim[0],1,MPI_DOUBLE,MPI_MIN,pmesh->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE,&rlim[1],1,MPI_DOUBLE,MPI_MAX,pmesh->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE,&alim[0],1,MPI_DOUBLE,MPI_MIN,pmesh->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE,&alim[1],1,MPI_DOUBLE,MPI_MAX,pmesh->GetComm());
+   }
+#endif
+}
+
+void ToroidPML::SetAttributes(Mesh *mesh_)
+{
+   int nrelem = mesh_->GetNE();
+   elems.SetSize(nrelem);
+
+   // Loop through the elements and identify which of them are in the PML
+   for (int i = 0; i < nrelem; ++i)
+   {
+      // initialize with 1
+      elems[i] = 1;
+      Element *el = mesh_->GetElement(i);
+      // Initialize attribute
+      el->SetAttribute(1);
+
+      Array<int> vertices;
+      el->GetVertices(vertices);
+      int nrvert = vertices.Size();
+      // Check if any vertex is in the pml
+      bool in_pml = false;
+      for (int iv = 0; iv < nrvert; ++iv)
+      {
+         int vert_idx = vertices[iv];
+         double *coords = mesh_->GetVertex(vert_idx);
+         double x = coords[0];
+         double y = coords[1];
+         double a = GetAngle(x,y);
+         double r = sqrt(x*x + y*y);
+
+         if (astretch)
+         {
+            if ( (a <= alim[0]+apml_thickness[0]) ||
+                 (a >= alim[1]-apml_thickness[1]) )
+            {
+               in_pml = true;
+               break;  
+            }
+         }
+
+         if (rstretch)
+         {
+            if ( (r <= rlim[0]+rpml_thickness[0]) ||
+                 (r >= rlim[1]-rpml_thickness[1]) )
+            {
+               in_pml = true;
+               break;  
+            }
+         }
+      }
+      if (in_pml)
+      {
+         elems[i] = 0;
+         el->SetAttribute(2);
+      }
+
+      // Vector center;
+      // mesh_->GetElementCenter(i,center);
+      // double x = center[0];
+      // double y = center[1];
+      // double a = GetAngle(x,y);
+      // double r = sqrt(x*x + y*y);
+      // // check upper and lower bound
+      // if (astretch)
+      // {
+      //    if ( (a <= alim[0]+apml_thickness[0]) ||
+      //      (a >= alim[1]-apml_thickness[1]) )
+      //    {
+      //       elems[i] = 0;
+      //       el->SetAttribute(2);
+      //    }
+      // }
+      // if (rstretch)
+      // {
+      //    if ( (r <= rlim[0]+rpml_thickness[0]) ||
+      //         (r >= rlim[1]-rpml_thickness[1]) )
+      //    {
+      //       elems[i] = 0;
+      //       el->SetAttribute(2);
+      //    }
+      // }
+   }
+   mesh_->SetAttributes();
+}
+
+
+double ToroidPML::GetAngle(const double x, const double y)
+{
+   // Find r and a for this point
+   double arad;
+   if (x == 0.0)
+   {
+      arad = (y > 0.0)? M_PI/2.0 : 3.0 * M_PI/2.0;
+   }
+   else
+   {
+      arad = atan(y/x);
+      int k = 0;
+      if (x<0)
+      {
+         k = 1;
+      }
+      else if (y<0)
+      {
+         k = 2;
+      }
+      arad += k*M_PI;
+   }
+   return arad * 180.0/M_PI;
+}
+
+// void ToroidPML::StretchFunction(const Vector &X,
+//                                    vector<complex<double>> &dxs, double omega)
+void ToroidPML::StretchFunction(const Vector &X, ComplexDenseMatrix & J, double omega)
+
+{
+   complex<double> zi = complex<double>(0., 1.);
+
+   double n = 2.0;
+   double c = 10.0;
+   // double c = log(omega);
+   // Stretch in the azimuthal direction
+   double x = X[0];
+   double y = X[1];
+   if (abs(x) < 1e-12) x = 0.0;
+   if (abs(y) < 1e-12) y = 0.0;
+   double a = GetAngle(x,y);
+   double r = sqrt(x*x + y*y);
+   // dxs[0] = 1.0;
+   // dxs[1] = 1.0;
+   J = 0.0;
+   J(0,0) = 1.0;
+   J(1,1) = 1.0;
+   if (dim == 3) J(2,2) = 1.0;
+
+   if (astretch)
+   {
+      double th = a * M_PI/180.0;
+      double thl, thL, thH;
+      bool in_pml = false;
+      // negative direction 
+      if (a <= alim[0]+apml_thickness[0])
+      {
+         in_pml = true;
+         thL = alim[1] * M_PI/180.0;
+         thH = apml_thickness[1] * M_PI/180.0;
+         thl = thL + thH;
+      }
+      // positive direction 
+      if (a >= alim[1]-apml_thickness[1])
+      {
+         in_pml = true;
+         thL = alim[1] * M_PI/180.0;
+         thH = apml_thickness[1] * M_PI/180.0;
+         thl = thL - thH;
+      }
+      // double c1 = min(20.0*M_PI/180.0,thH);
+      if (in_pml)
+      {
+         double c1 = thH;
+         double coeff = n * c / omega / pow(c1,n);
+         double f_th = pow(th - thl,n-1);
+         double th_x = - y / (r * r);
+         double th_y =   x / (r * r);
+
+         J(0,0) = 1.0 + zi * coeff * abs(f_th * th_x);
+         J(0,1) = zi * f_th * th_y;
+         J(1,0) = zi * f_th * th_x;
+         J(1,1) = 1.0 + zi * coeff * abs(f_th * th_y);
+      }
+   }
+   // Stretch in the radial direction
+   if (rstretch)
+   {  // negative
+      double rl, rL, rH;
+      bool in_pml = false;
+      if (r <= rlim[0]+rpml_thickness[0])
+      {
+         in_pml = true;
+         rL = rlim[0];     
+         rH = rpml_thickness[0];
+         rl = rL + rH;
+      }
+      // positive direction 
+      if (r >= rlim[1]-rpml_thickness[1])
+      {
+         in_pml = true;
+         rL = rlim[1];     
+         rH = rpml_thickness[1];
+         rl = rL - rH;
+      }
+
+      if (in_pml)
+      {
+         double coeff = n * c / omega / pow (rH,n);
+         double f_r =  pow(r-rl,n-1.0);
+         double r_x = x / r;
+         double r_y = y / r;
+
+         J(0,0) = 1.0 + zi * coeff * abs(f_r*r_x);
+         J(0,1) = zi * f_r * r_y;
+         J(1,0) = zi * f_r * r_x;
+         J(1,1) = 1.0 + zi * coeff * abs(f_r*r_y);
+      }
+   }
+}
+
+double pml_detJ_Re(const Vector & x, CartesianPML * pml)
+{
+   int dim = pml->dim;
+   double omega = pml->omega;
+   std::vector<std::complex<double>> dxs(dim);
+   complex<double> det(1.0,0.0);
+   pml->StretchFunction(x, dxs, omega);
+   for (int i=0; i<dim; ++i) det *= dxs[i];
+   return det.real();
+}
+
+double pml_detJ_Im(const Vector & x, CartesianPML * pml)
+{
+   int dim = pml->dim;
+   double omega = pml->omega;
+   std::vector<std::complex<double>> dxs(dim);
+   complex<double> det(1.0,0.0);
+   pml->StretchFunction(x, dxs, omega);
+   for (int i=0; i<dim; ++i) det *= dxs[i];
+   return det.imag();
+}
+
+void pml_detJ_JT_J_inv_Re(const Vector & x, CartesianPML * pml , DenseMatrix & M)
+{
+   int dim = pml->dim;
+   double omega = pml->omega;
+   std::vector<std::complex<double>> dxs(dim);
+   complex<double> det(1.0,0.0);
+   pml->StretchFunction(x, dxs, omega);
+
+   for (int i = 0; i<dim; ++i)
+   {
+      det *= dxs[i];
+   }
+
+   M=0.0;
+   for (int i = 0; i<dim; ++i)
+   {
+      M(i,i) = (det / pow(dxs[i],2)).real();
+   }
+}
+
+void pml_detJ_JT_J_inv_Im(const Vector & x, CartesianPML * pml , DenseMatrix & M)
+{
+   int dim = pml->dim;
+   double omega = pml->omega;
+
+   std::vector<std::complex<double>> dxs(dim);
+   complex<double> det = 1.0;
+   pml->StretchFunction(x, dxs, omega);
+
+   for (int i = 0; i<dim; ++i)
+   {
+      det *= dxs[i];
+   }
+
+   M=0.0;
+   for (int i = 0; i<dim; ++i)
+   {
+      M(i,i) = (det / pow(dxs[i],2)).imag();
+   }
+}
+
+
+void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+{
+   int dim = pml->dim;
+   double omega = pml->omega;
+   vector<complex<double>> dxs(dim);
+   complex<double> det(1.0, 0.0);
+   pml->StretchFunction(x, dxs, omega);
+
+   for (int i = 0; i < dim; ++i)
+   {
+      det *= dxs[i];
+   }
+
+   M = 0.0;
+   for (int i = 0; i < dim; ++i)
+   {
+      M(i, i) = (det / pow(dxs[i], 2)).real();
+   }
+}
+
+void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+{
+   int dim = pml->dim;
+   double omega = pml->omega;
+   vector<complex<double>> dxs(dim);
+   complex<double> det = 1.0;
+   pml->StretchFunction(x, dxs, omega);
+
+   for (int i = 0; i < dim; ++i)
+   {
+      det *= dxs[i];
+   }
+
+   M = 0.0;
+   for (int i = 0; i < dim; ++i)
+   {
+      M(i, i) = (det / pow(dxs[i], 2)).imag();
+   }
+}
+
+void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+{
+   int dim = pml->dim;
+   double omega = pml->omega;
+   vector<complex<double>> dxs(dim);
+   complex<double> det = 1.0;
+   pml->StretchFunction(x, dxs, omega);
+
+   for (int i = 0; i < dim; ++i)
+   {
+      det *= dxs[i];
+   }
+
+   M = 0.0;
+   for (int i = 0; i < dim; ++i)
+   {
+      M(i, i) = abs(det / pow(dxs[i], 2));
+   }
+}
+
+void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+{
+   int dim = pml->dim;
+   double omega = pml->omega;
+   vector<complex<double>> dxs(dim);
+   complex<double> det(1.0, 0.0);
+   pml->StretchFunction(x, dxs, omega);
+
+   for (int i = 0; i < dim; ++i)
+   {
+      det *= dxs[i];
+   }
+
+   // in the 2D case the coefficient is scalar 1/det(J)
+   if (dim == 2)
+   {
+      M = (1.0 / det).real();
+   }
+   else
+   {
+      M = 0.0;
+      for (int i = 0; i < dim; ++i)
+      {
+         M(i, i) = (pow(dxs[i], 2) / det).real();
+      }
+   }
+}
+
+void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+{
+   int dim = pml->dim;
+   double omega = pml->omega;
+   vector<complex<double>> dxs(dim);
+   complex<double> det = 1.0;
+   pml->StretchFunction(x, dxs, omega);
+
+   for (int i = 0; i < dim; ++i)
+   {
+      det *= dxs[i];
+   }
+
+   if (dim == 2)
+   {
+      M = (1.0 / det).imag();
+   }
+   else
+   {
+      M = 0.0;
+      for (int i = 0; i < dim; ++i)
+      {
+         M(i, i) = (pow(dxs[i], 2) / det).imag();
+      }
+   }
+}

examples/maxwell-solver-dev/common/PML.hpp

@@ -0,0 +1,179 @@
+#pragma once
+#include "mfem.hpp"
+#include "complex_linalg.hpp"
+#include <fstream>
+#include <iostream>
+using namespace std;
+using namespace mfem;
+
+
+// Class for setting up a simple Cartesian PML region
+class CartesianPML
+{
+private:
+   Mesh *mesh;
+
+   // Length of the PML Region in each direction
+   Array2D<double> length;
+
+   // Computational Domain Boundary
+   Array2D<double> comp_dom_bdr;
+
+   // Domain Boundary
+   Array2D<double> dom_bdr;
+
+   // Integer Array identifying elements in the pml
+   // 0: in the pml, 1: not in the pml
+   Array<int> elems;
+
+   // Compute Domain and Computational Domain Boundaries
+   void SetBoundaries();
+
+public:
+   // Constructor
+   CartesianPML(Mesh *mesh_,Array2D<double> length_);
+
+   int dim;
+   double omega;
+   // Return Computational Domain Boundary
+   Array2D<double> GetCompDomainBdr() {return comp_dom_bdr;}
+
+   // Return Domain Boundary
+   Array2D<double> GetDomainBdr() {return dom_bdr;}
+
+   // Return Marker list for elements
+   Array<int> * GetMarkedPMLElements() {return &elems;}
+
+   // Mark element in the PML region
+   void SetAttributes(Mesh *mesh_);
+
+   void SetOmega(double omega_) {omega = omega_;}
+
+   // PML complex stretching function
+   void StretchFunction(const Vector &x, vector<complex<double>> &dxs, double omega);
+};
+
+class ToroidPML
+{
+private:
+   Mesh *mesh;
+
+   Vector zlim, zpml_thickness; // range in axial direction
+   Vector rlim, rpml_thickness; // range in radial direction
+   Vector alim, apml_thickness; // range in azimuthal direction
+
+   // Integer Array identifying elements in the pml
+   // 0: in the pml, 1: not in the pml
+   Array<int> elems;
+
+   double GetAngle(const double x, const double y);
+
+   // Compute Domain and Computational Domain Boundaries
+   void SetBoundaries();
+
+   bool zstretch = false;
+   bool rstretch = false;
+   bool astretch = false;
+
+public:
+   // Constructor
+   ToroidPML(Mesh *mesh_);
+
+   int dim;
+   double omega;
+   // Return Computational Domain Boundary
+
+   // Return Domain Boundary
+   void GetDomainBdrs(Vector & zlim_, Vector & rlim_, Vector & alim_)
+   {
+      zlim_.SetSize(2); zlim_ = zlim;
+      rlim_.SetSize(2); rlim_ = rlim;
+      alim_.SetSize(2); alim_ = alim;
+   }
+
+   void SetPmlWidth(const Vector & zpml, const Vector & rpml, const Vector & apml)
+   {
+      MFEM_VERIFY(zpml.Size() == 2 , "Check zpml size");
+      MFEM_VERIFY(rpml.Size() == 2 , "Check rpml size");
+      MFEM_VERIFY(apml.Size() == 2 , "Check apml size");
+      zpml_thickness = zpml;
+      rpml_thickness = rpml;
+      apml_thickness = apml;
+   }
+
+   void SetPmlAxes(const bool zstretch_, 
+                   const bool rstretch_, 
+                   const bool astretch_ )
+   {
+      zstretch = zstretch_;
+      rstretch = rstretch_;
+      astretch = astretch_;
+   }
+
+   // // Return Marker list for elements
+   Array<int> * GetMarkedPMLElements() {return &elems;}
+
+   // Mark element in the PML region
+   void SetAttributes(Mesh *mesh_);
+
+   void SetOmega(double omega_) {omega = omega_;}
+
+   // PML complex stretching function
+   // void StretchFunction(const Vector &X, vector<complex<double>> &dxs, double omega);
+   void StretchFunction(const Vector &X, ComplexDenseMatrix & J, double omega);
+};
+
+class PmlCoefficient : public Coefficient
+{
+private:
+   CartesianPML * pml = nullptr;
+   double (*Function)(const Vector &, CartesianPML * );
+public:
+   PmlCoefficient(double (*F)(const Vector &, CartesianPML *), CartesianPML * pml_)
+      : pml(pml_), Function(F)
+   {}
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+      T.Transform(ip, transip);
+      return ((*Function)(transip, pml));
+   }
+};
+
+
+// This includes scalar coefficients
+class PmlMatrixCoefficient : public MatrixCoefficient
+{
+private:
+   CartesianPML * pml = nullptr;
+   void (*Function)(const Vector &, CartesianPML * , DenseMatrix &);
+public:
+   PmlMatrixCoefficient(int dim, void(*F)(const Vector &, CartesianPML *,
+                                          DenseMatrix &),
+                        CartesianPML * pml_)
+      : MatrixCoefficient(dim), pml(pml_), Function(F)
+   {}
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+      T.Transform(ip, transip);
+      K.SetSize(height, width);
+      (*Function)(transip, pml, K);
+   }
+};
+
+// Helmholtz pml Functions
+double pml_detJ_Re(const Vector & x, CartesianPML * pml);
+double pml_detJ_Im(const Vector & x, CartesianPML * pml);
+void pml_detJ_JT_J_inv_Re(const Vector & x, CartesianPML * pml , DenseMatrix & M);
+void pml_detJ_JT_J_inv_Im(const Vector & x, CartesianPML * pml , DenseMatrix & M);
+
+// Maxwell Pml functions
+void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M);
+void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M);
+void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M);
+void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M);
+void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M);

examples/maxwell-solver-dev/common/Utilities.cpp

@@ -0,0 +1,619 @@
+#include "Utilities.hpp"
+
+Sweep::Sweep(int dim_) : dim(dim_)
+{
+   nsweeps = pow(2,dim);
+   sweeps.resize(nsweeps);
+
+   for (int is = 0; is<nsweeps; is++)
+   {
+      sweeps[is].SetSize(dim);
+   }
+
+   switch(dim)
+   {
+      case 1: 
+         sweeps[0][0] =  1; 
+         sweeps[1][0] = -1; 
+         break;
+      case 2: 
+         sweeps[0][0] =  1; sweeps[0][1] =  1; 
+         sweeps[1][0] = -1; sweeps[1][1] =  1; 
+         sweeps[2][0] =  1; sweeps[2][1] = -1; 
+         sweeps[3][0] = -1; sweeps[3][1] = -1;
+         break;
+      default:
+         sweeps[0][0] =  1; sweeps[0][1] =  1; sweeps[0][2] =  1; 
+         sweeps[1][0] = -1; sweeps[1][1] =  1; sweeps[1][2] =  1; 
+         sweeps[2][0] =  1; sweeps[2][1] = -1; sweeps[2][2] =  1; 
+         sweeps[3][0] = -1; sweeps[3][1] = -1; sweeps[3][2] =  1;
+         sweeps[4][0] =  1; sweeps[4][1] =  1; sweeps[4][2] = -1; 
+         sweeps[5][0] = -1; sweeps[5][1] =  1; sweeps[5][2] = -1; 
+         sweeps[6][0] =  1; sweeps[6][1] = -1; sweeps[6][2] = -1; 
+         sweeps[7][0] = -1; sweeps[7][1] = -1; sweeps[7][2] = -1;
+         break;
+   }
+}
+
+Sweep::~Sweep()
+{
+   for (int i = 0; i<nsweeps; i++)
+   {
+      sweeps[i].DeleteAll();
+   }
+}
+
+
+double CutOffFncn(const Vector &x, const Vector & pmin, const Vector & pmax, const Array2D<double> & h_)
+{
+   int dim = pmin.Size();
+   Vector h0(dim);
+   Vector h1(dim);
+   for (int i=0; i<dim; i++)
+   {
+      h0(i) = h_[i][0];
+      h1(i) = h_[i][1];
+   }
+   Vector x0(dim);
+   Vector x1(dim);
+   x0 = pmin; x0+=h0;
+   x1 = pmax; x1-=h1;
+
+   double f = 1.0;
+   for (int i = 0; i<dim; i++)
+   {
+      double val = 1.0;
+      if( x(i) >= pmax(i) || x(i) <= pmin(i))
+      {
+         val = 0.0;
+      }  
+      else if (x(i) < pmax(i) && x(i) >= x1(i))
+      {
+         if(h1(i) != 0.0)
+            // val = (x(i)-pmax(i))/(x1(i)-pmax(i)); 
+            val = pow((x(i)-pmax(i))/(x1(i)-pmax(i)),1.0); 
+      }
+      else if (x(i) > pmin(i) && x(i) <= x0(i))
+      {
+         if (h0(i) != 0.0)
+            // val = (x(i)-pmin(i))/(x0(i)-pmin(i)); 
+            val = pow((x(i)-pmin(i))/(x0(i)-pmin(i)),1.0); 
+      }
+
+      if (h0(i) == 0 && x(i) <= x1(i))
+      {
+         val = 1.0;
+      }
+      if (h1(i) == 0 && x(i) >= x0(i))
+      {
+         val = 1.0;
+      }
+      f *= val;
+   }
+   return f;
+}
+
+double ChiFncn(const Vector &x, const Vector & pmin, const Vector & pmax, const Array2D<double> & h_)
+{
+   int dim = pmin.Size();
+   Vector h0(dim);
+   Vector h1(dim);
+   for (int i=0; i<dim; i++)
+   {
+      h0(i) = h_[i][0];
+      h1(i) = h_[i][1];
+   }
+   Vector x0(dim);
+   Vector x1(dim);
+   x0 = pmin; x0+=h0;
+   x1 = pmax; x1-=h1;
+
+   double f = 1.0;
+   for (int i = 0; i<dim; i++)
+   {
+      double val = 1.0;
+      if( x(i) >= pmax(i) || x(i) <= pmin(i))
+      {
+         val = 0.0;
+      }  
+      else if (x(i) < pmax(i) && x(i) >= x1(i))
+      {
+         if(h1(i) != 0.0)
+            val = (x(i)-pmax(i))/(x1(i)-pmax(i)); 
+            // This function has to be changed to smth more reasonable
+            // val = pow((x(i)-pmax(i))/(x1(i)-pmax(i)),100.0); 
+      }
+      else if (x(i) > pmin(i) && x(i) <= x0(i))
+      {
+         if (h0(i) != 0.0)
+            val = (x(i)-pmin(i))/(x0(i)-pmin(i)); 
+            // val = pow((x(i)-pmin(i))/(x0(i)-pmin(i)),100.0); 
+      }
+
+      if (h0(i) == 0 && x(i) <= x1(i))
+      {
+         val = 1.0;
+      }
+      if (h1(i) == 0 && x(i) >= x0(i))
+      {
+         val = 1.0;
+      }
+      f *= val;
+   }
+   return f;
+}
+
+
+DofMap::DofMap(FiniteElementSpace * fes , MeshPartition * partition) 
+{
+   const FiniteElementCollection * fec = fes->FEColl();
+   nrpatch = partition->nrpatch;
+
+   fespaces.SetSize(nrpatch);
+
+   Dof2GlobalDof.resize(nrpatch);
+
+   for (int ip=0; ip<nrpatch; ++ip)
+   {
+      // create finite element spaces for each patch 
+      fespaces[ip] = new FiniteElementSpace(partition->patch_mesh[ip],fec);
+
+      // construct the patch tdof to global tdof map
+      int nrdof = fespaces[ip]->GetTrueVSize();
+      Dof2GlobalDof[ip].SetSize(2*nrdof);
+
+      // loop through the elements in the patch
+      for (int iel = 0; iel<partition->element_map[ip].Size(); ++iel)
+      {
+         // index in the global mesh
+         int iel_idx = partition->element_map[ip][iel];
+         // get the dofs of this element
+         Array<int> ElemDofs;
+         Array<int> GlobalElemDofs;
+         fespaces[ip]->GetElementDofs(iel,ElemDofs);
+         fes->GetElementDofs(iel_idx,GlobalElemDofs);
+         // the sizes have to match
+         MFEM_VERIFY(ElemDofs.Size() == GlobalElemDofs.Size(),
+                     "Size inconsistency");
+         // loop through the dofs and take into account the signs;
+         int ndof = ElemDofs.Size();
+         for (int i = 0; i<ndof; ++i)
+         {
+            int pdof_ = ElemDofs[i];
+            int gdof_ = GlobalElemDofs[i];
+            int pdof = (pdof_ >= 0) ? pdof_ : abs(pdof_) - 1;
+            int gdof = (gdof_ >= 0) ? gdof_ : abs(gdof_) - 1;
+            Dof2GlobalDof[ip][pdof] = gdof;
+            Dof2GlobalDof[ip][pdof+nrdof] = gdof+fes->GetTrueVSize();
+         }
+      }
+   }
+}
+
+DofMap::DofMap(FiniteElementSpace * fes , MeshPartition * partition, int nrlayers) 
+{
+
+   nx = partition->nxyz[0];
+   ny = partition->nxyz[1];
+   nz = partition->nxyz[2];
+
+   int partition_kind = partition->partition_kind;
+   // Mesh * mesh = fespace->GetMesh();
+   const FiniteElementCollection * fec = fes->FEColl();
+   nrpatch = partition->nrpatch;
+
+   fespaces.SetSize(nrpatch);
+   PmlMeshes.SetSize(nrpatch);
+   // Extend patch meshes to include pml
+
+   for  (int ip = 0; ip<nrpatch; ip++)
+   {
+      int k = ip/(nx*ny);
+      int j = (ip-k*nx*ny)/nx;
+      int i = (ip-k*nx*ny)%nx;
+      
+      Array<int> directions;
+      if (i > 0)
+      {
+         for (int i=0; i<nrlayers; i++)
+         {
+            directions.Append(-1);
+         }
+      }
+      if (j > 0)
+      {
+         for (int i=0; i<nrlayers; i++)
+         {
+            directions.Append(-2);
+         }
+      }
+      if (k > 0)
+      {
+         for (int i=0; i<nrlayers; i++)
+         {
+            directions.Append(-3);
+         }
+      }
+      if (i < nx-1)
+      {
+         for (int i=0; i<nrlayers; i++)
+         {
+            if (partition_kind == 3 || partition_kind == 2) directions.Append(1);
+         }
+      }
+      if (j < ny-1)
+      {
+         for (int i=0; i<nrlayers; i++)
+         {
+            if (partition_kind == 3 || partition_kind == 2) directions.Append(2);
+         }
+      }
+      if (k < nz-1)
+      {
+         for (int i=0; i<nrlayers; i++)
+         {
+            if (partition_kind == 3 || partition_kind == 2) directions.Append(1);
+         }
+      }
+      PmlMeshes[ip] = ExtendMesh(partition->patch_mesh[ip],directions);
+   }
+
+   // Save PML_meshes
+   string meshpath;
+   string solpath;
+   if (partition_kind == 3 || partition_kind == 2) 
+   {
+      meshpath = "output/mesh_ovlp_pml.";
+      solpath = "output/sol_ovlp_pml.";
+   }
+   else if (partition_kind == 4)
+   {
+      meshpath = "output/mesh_novlp_pml.";
+      solpath = "output/sol_novlp_pml.";
+   }
+   else
+   {
+      MFEM_ABORT("This partition kind not supported yet");
+   }
+   
+   // SaveMeshPartition(PmlMeshes, meshpath, solpath);
+
+   PmlFespaces.SetSize(nrpatch);
+   Dof2GlobalDof.resize(nrpatch);
+   Dof2PmlDof.resize(nrpatch);
+
+   for (int ip=0; ip<nrpatch; ++ip)
+   {
+      // create finite element spaces for each patch 
+      fespaces[ip] = new FiniteElementSpace(partition->patch_mesh[ip],fec);
+      PmlFespaces[ip] = new FiniteElementSpace(PmlMeshes[ip],fec);
+
+      // construct the patch tdof to global tdof map
+      int nrdof = fespaces[ip]->GetTrueVSize();
+      Dof2GlobalDof[ip].SetSize(2*nrdof);
+      Dof2PmlDof[ip].SetSize(2*nrdof);
+
+      // build dof maps between patch and extended patch
+      // loop through the patch elements and constract the dof map
+      // The same elements in the extended mesh have the same ordering (but not the dofs)
+
+      // loop through the elements in the patch
+      for (int iel = 0; iel<partition->element_map[ip].Size(); ++iel)
+      {
+         // index in the global mesh
+         int iel_idx = partition->element_map[ip][iel];
+         // get the dofs of this element
+         Array<int> ElemDofs;
+         Array<int> PmlElemDofs;
+         Array<int> GlobalElemDofs;
+         fespaces[ip]->GetElementDofs(iel,ElemDofs);
+         PmlFespaces[ip]->GetElementDofs(iel,PmlElemDofs);
+         fes->GetElementDofs(iel_idx,GlobalElemDofs);
+         // the sizes have to match
+         MFEM_VERIFY(ElemDofs.Size() == GlobalElemDofs.Size(),
+                     "Size inconsistency");
+         MFEM_VERIFY(ElemDofs.Size() == PmlElemDofs.Size(),
+                     "Size inconsistency");            
+         // loop through the dofs and take into account the signs;
+         int ndof = ElemDofs.Size();
+         for (int i = 0; i<ndof; ++i)
+         {
+            int pdof_ = ElemDofs[i];
+            int gdof_ = GlobalElemDofs[i];
+            int pmldof_ = PmlElemDofs[i];
+            int pdof = (pdof_ >= 0) ? pdof_ : abs(pdof_) - 1;
+            int gdof = (gdof_ >= 0) ? gdof_ : abs(gdof_) - 1;
+            int pmldof = (pmldof_ >= 0) ? pmldof_ : abs(pmldof_) - 1;
+
+            Dof2GlobalDof[ip][pdof] = gdof;
+            Dof2GlobalDof[ip][pdof+nrdof] = gdof+fes->GetTrueVSize();
+            Dof2PmlDof[ip][pdof] = pmldof;
+            Dof2PmlDof[ip][pdof+nrdof] = pmldof+PmlFespaces[ip]->GetTrueVSize();
+         }
+      }
+   }
+}
+
+
+LocalDofMap::LocalDofMap(const FiniteElementCollection * fec_, MeshPartition * part1_, 
+               MeshPartition * part2_):fec(fec_), part1(part1_), part2(part2_)
+{
+   // Each overlapping patch has 2 non-overlapping subdomains
+   // Thre are n non-overlapping and and n-1 overlapping subdomains
+   int nrpatch = part2->nrpatch;
+   MFEM_VERIFY(part1->nrpatch-1 == part2->nrpatch, "Check number of subdomains");
+
+   cout << "Constructing local dof maps" << endl; 
+   map1.resize(nrpatch);
+   map2.resize(nrpatch);
+   for (int ip=0; ip<nrpatch; ip++)
+   {
+      // Get the 3 meshes involved
+      Mesh * mesh = part2->patch_mesh[ip];
+      Mesh * mesh1 = part1->patch_mesh[ip];
+      Mesh * mesh2 = part1->patch_mesh[ip+1];
+
+      // Define the fespaces
+      FiniteElementSpace fespace(mesh, fec);
+      FiniteElementSpace fespace1(mesh1, fec);
+      FiniteElementSpace fespace2(mesh2, fec);
+
+      int ndof1 = fespace1.GetTrueVSize();
+      int ndof2 = fespace2.GetTrueVSize();
+
+      map1[ip].SetSize(2*ndof1); // times 2 because it's complex
+      map2[ip].SetSize(2*ndof2); // times 2 because it's complex
+
+      // loop through the elements in the patches
+      // map 1 is constructed by the first half of elements
+      // map 2 is constructed by the second half of elements
+
+      for (int iel = 0; iel<part1->element_map[ip].Size(); ++iel)
+      {
+         // index in the overlapping mesh
+         int iel_idx = iel;
+         Array<int> ElemDofs;
+         Array<int> GlobalElemDofs;
+         fespace1.GetElementDofs(iel,ElemDofs);
+         fespace.GetElementDofs(iel_idx,GlobalElemDofs);
+         // the sizes have to match
+         MFEM_VERIFY(ElemDofs.Size() == GlobalElemDofs.Size(),
+                     "Size inconsistency");
+         // loop through the dofs and take into account the signs;
+         int ndof = ElemDofs.Size();
+         for (int i = 0; i<ndof; ++i)
+         {
+            int pdof_ = ElemDofs[i];
+            int gdof_ = GlobalElemDofs[i];
+            int pdof = (pdof_ >= 0) ? pdof_ : abs(pdof_) - 1;
+            int gdof = (gdof_ >= 0) ? gdof_ : abs(gdof_) - 1;
+            map1[ip][pdof] = gdof;
+            map1[ip][pdof+ndof1] = gdof+fespace.GetTrueVSize();
+         }
+      }
+      for (int iel = 0; iel<part1->element_map[ip+1].Size(); ++iel)
+      {
+         // index in the overlapping mesh
+         int k = part1->element_map[ip].Size();
+         int iel_idx = iel+k;
+         Array<int> ElemDofs;
+         Array<int> GlobalElemDofs;
+         fespace2.GetElementDofs(iel,ElemDofs);
+         fespace.GetElementDofs(iel_idx,GlobalElemDofs);
+         // the sizes have to match
+         MFEM_VERIFY(ElemDofs.Size() == GlobalElemDofs.Size(),
+                     "Size inconsistency");
+         // loop through the dofs and take into account the signs;
+         int ndof = ElemDofs.Size();
+         for (int i = 0; i<ndof; ++i)
+         {
+            int pdof_ = ElemDofs[i];
+            int gdof_ = GlobalElemDofs[i];
+            int pdof = (pdof_ >= 0) ? pdof_ : abs(pdof_) - 1;
+            int gdof = (gdof_ >= 0) ? gdof_ : abs(gdof_) - 1;
+            map2[ip][pdof] = gdof;
+            map2[ip][pdof+ndof2] = gdof+fespace.GetTrueVSize();
+         }
+      }
+   }
+};
+
+
+NeighborDofMaps::NeighborDofMaps(MeshPartition * part_, FiniteElementSpace * fes_,
+                                DofMap * dmap_, 
+                                int ovlp_layers_) : part(part_), fes(fes_), 
+                                dmap(dmap_),
+                                ovlp_layers(ovlp_layers_)
+{
+
+   nrsubdomains = part->nrpatch;
+   nxyz.SetSize(3);
+   mesh = fes->GetMesh();
+   dim = mesh->Dimension();
+   for (int d=0; d<3; d++)  nxyz[d] = part->nxyz[d];
+   MarkOvlpElements();
+   ComputeNeighborDofMaps();
+}
+
+void NeighborDofMaps::MarkOvlpElements()
+{
+   // Lists of elements
+   // x,y,z = +/- 1 ovlp 
+   OvlpElems.resize(nrsubdomains);
+
+   for (int ip = 0; ip<nrsubdomains; ip++)
+   {
+      int i0,j0,k0;
+      Getijk(ip,i0,j0,k0);
+      int ijk[dim]; ijk[0] = i0; ijk[1]=j0;
+      if (dim==3) ijk[2] = k0;
+
+      FiniteElementSpace * sub_fes = dmap->fespaces[ip];
+      Mesh * sub_mesh = sub_fes->GetMesh();
+      // OvlpElems[ip].resize(2*dim);
+      OvlpElems[ip].resize(pow(3,dim)); 
+
+      Vector pmin, pmax;
+      sub_mesh->GetBoundingBox(pmin,pmax);
+      double h = part->MeshSize;
+      // Loop through elements
+      for (int iel=0; iel<sub_mesh->GetNE(); iel++)
+      {
+         // Get element center
+         Vector center(dim);
+         int geom = sub_mesh->GetElementBaseGeometry(iel);
+         ElementTransformation * tr = sub_mesh->GetElementTransformation(iel);
+         tr->Transform(Geometries.GetCenter(geom),center);
+
+         // loop through dimensions   
+         Array<bool> pos(dim); pos = 0;
+         Array<bool> neg(dim); neg = 0;
+
+         for (int d=0;d<dim; d++)
+         {
+            if (ijk[d]>0 && center[d] < pmin[d]+2.0*h*ovlp_layers)
+            {
+               neg[d] = true;
+            }
+
+            if (ijk[d]<nxyz[d]-1 && center[d] > pmax[d]-2.0*h*ovlp_layers) 
+            {
+               pos[d] = true;
+            }
+         }
+         SetElementToOverlap(ip,iel,neg,pos);        
+      }   
+   }
+}
+
+void NeighborDofMaps::ComputeNeighborDofMaps()
+{
+   OvlpDofMaps.resize(nrsubdomains);
+
+   // Array<UniqueIndexGen * > Gen(nrsubdomains);
+   // // construct unique number generator for the elements of a patch
+   // for (int ip = 0; ip<nrsubdomains; ip++)
+   // {
+   //    Gen[ip] = new UniqueIndexGen;
+   //    // register the elements
+   //    int nel = part->element_map[ip].Size();
+   //    for (int iel=0; iel<nel; iel++)
+   //    {
+   //       int iel_idx = part->element_map[ip][iel];
+   //       Gen[ip]->Set(iel_idx);
+   //    }
+   // }
+
+   // construct dof maps
+   int nrneighbors = pow(3,dim); // including its self
+
+   for (int ip0 = 0; ip0<nrsubdomains; ip0++)
+   {
+      OvlpDofMaps[ip0].resize(nrneighbors);
+
+      FiniteElementSpace * fes0 = dmap->fespaces[ip0];
+      int tdofs0 = fes0->GetTrueVSize();
+      Array<int> marker0(tdofs0); marker0 = 0;
+      int i0, j0, k0;
+      Array<int> ijk(dim);
+      Getijk(ip0, i0,j0,k0);
+
+      int kbeg = (dim == 2) ? 0 : -1;
+      int kend = (dim == 2) ? 1 :  2;
+      for (int k=kbeg; k<kend; k++)
+      {
+         int k1 = k0 + k;
+         if (k1 <0 || k1>=nxyz[2]) continue;
+         int kk = (dim == 2) ? -1 : k;
+         for (int j=-1; j<2; j++)
+         {
+            int j1 = j0 + j;
+            if (j1 <0 || j1>=nxyz[1]) continue;
+            for (int i=-1; i<2; i++)
+            {
+               int i1 = i0 + i;
+               if (i1 <0 || i1>=nxyz[0]) continue;
+     
+               Array<int> ip0list;  marker0 = 0;
+               int directionId = GetDirectionId(i,j,kk);
+
+               Array<int> Elems = OvlpElems[ip0][directionId];
+               int nel = Elems.Size();
+
+               for (int iel = 0; iel<nel; ++iel)
+               {
+                  int iel0 = Elems[iel];   
+                  Array<int> ElemDofs0;
+
+                  fes0->GetElementDofs(iel0,ElemDofs0);
+                  int ndof = ElemDofs0.Size();
+                  // since the elements are added to the subdomain meshes 
+                  // in the same ordered fashion (as they come from the
+                  // original mesh) then the ordering of elements in each 
+                  // subdomain is the same. Hence the dof ovlp lists
+                  // can be computed for each subdomain independendly
+                  for (int l = 0; l<ndof; ++l)
+                  {
+                     int dof0_ = ElemDofs0[l];
+                     int dof0 = (dof0_ >= 0) ? dof0_ : abs(dof0_) - 1;
+                     if (!marker0[dof0])
+                     {
+                        ip0list.Append(dof0); // dofs of ip0 in ovlp
+                        marker0[dof0] = 1;
+                     }
+                  }
+               }
+               
+               OvlpDofMaps[ip0][directionId].Append(ip0list);
+               int tsize = fes0->GetTrueVSize();
+               // Imaginary part
+               for (int l=0;l<ip0list.Size(); l++) { ip0list[l] += tsize; }
+               OvlpDofMaps[ip0][directionId].Append(ip0list);
+            }
+         }
+      }
+   }
+}
+
+
+void NeighborDofMaps::GetNeighborDofMap(const int ip,
+                                        const Array<int> & directions, 
+                                        Array<int> & dofmap)
+{
+   int k = (dim == 2) ? -1 : directions[2];
+   int directionid = GetDirectionId(directions[0],directions[1],k);
+   dofmap = OvlpDofMaps[ip][directionid];
+}
+
+
+void NeighborDofMaps::SetElementToOverlap(int ip, int iel,
+                        const Array<bool> & neg, 
+                        const Array<bool> & pos)
+{
+   int kbeg = (dim == 2) ? 0 : -1;
+   int kend = (dim == 2) ? 0 :  1;
+   for (int k = kbeg; k<=kend; k++)
+   {
+      if (dim == 3)
+      {
+         if (k == -1 && !neg[2]) continue;   
+         if (k ==  1 && !pos[2]) continue;   
+      }   
+      for (int j = -1; j<=1; j++)
+      {
+         if (j== -1 && !neg[1]) continue;   
+         if (j==  1 && !pos[1]) continue;   
+         for (int i = -1; i<=1; i++)
+         {
+            // cases to skip 
+            if (i==-1 && !neg[0]) continue;   
+            if (i== 1 && !pos[0]) continue;   
+
+            if (i==0 && j==0 && k == 0) continue;
+            int kk  = (dim==2)?-1 : k;  
+            int DirId = GetDirectionId(i,j,kk);
+            OvlpElems[ip][DirId].Append(iel);
+         }
+      }
+   }
+}

examples/maxwell-solver-dev/common/Utilities.hpp

@@ -0,0 +1,179 @@
+#pragma once
+#include "MeshPartition.hpp"
+#include "complex_linalg.hpp"
+
+struct UniqueIndexGen
+{
+   int counter = 0;
+   std::unordered_map<int,int> idx;
+
+   void Set(int i)
+   {
+      std::unordered_map<int,int>::iterator f = idx.find(i);
+      if (f == idx.end())
+      {
+         idx[i] = counter;
+         counter++;
+      }
+   }
+
+   int Get(int i)
+   {
+      std::unordered_map<int,int>::iterator f = idx.find(i);
+      if (f == idx.end())
+      {
+         return -1;
+      }
+      else
+      {
+         return (*f).second;
+      }
+   }
+   void Reset()
+   {
+      counter = 0;
+      idx.clear();
+   }
+};
+
+
+struct Sweep
+{
+private:
+   int dim;
+   std::vector<Array<int>> sweeps;
+public:   
+   int nsweeps;
+   Sweep(int dim_);
+   ~Sweep();
+   void GetSweep(const int i, Array<int> & sweep)
+   {
+      MFEM_VERIFY(i<nsweeps, "Sweep number out of bounds");
+      sweep.SetSize(dim);
+      sweep = sweeps[i];
+   }
+
+
+};
+
+
+
+// Function coefficient that takes the bounding box of the mesh as an input
+class CutOffFnCoefficient : public Coefficient
+{
+private:
+   double (*Function)(const Vector &, const Vector &, const Vector &, const Array2D<double> &);
+   Vector pmin, pmax;
+   Array2D<double> h; // specify the with of the cutoff function (h in each direction)
+   
+
+public:
+   CutOffFnCoefficient(double (*F)(const Vector &, const Vector &, const Vector &, const Array2D<double> &), 
+                             const Vector & pmin_, const Vector & pmax_, Array2D<double> & h_)
+      : Function(F), pmin(pmin_), pmax(pmax_), h(h_)
+   {}
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+      T.Transform(ip, transip);
+      return ((*Function)(transip, pmin, pmax, h));
+   }
+};
+
+double CutOffFncn(const Vector &x, const Vector & pmin, 
+const Vector & pmax, const Array2D<double> & h_);
+double ChiFncn(const Vector &x, const Vector & pmin, 
+const Vector & pmax, const Array2D<double> & h_);
+
+class DofMap // Constructs dof maps for a given partition
+{
+public:
+   int nrpatch, nx, ny, nz;
+   vector<Array<int>> Dof2GlobalDof;
+   vector<Array<int>> Dof2PmlDof;
+   Array<Mesh *> PmlMeshes;
+   Array<FiniteElementSpace *> fespaces;
+   Array<FiniteElementSpace *> PmlFespaces;
+   // constructor
+   // Non PML constructor dof map
+   DofMap(FiniteElementSpace * fes, MeshPartition * partition);
+   // PML
+   DofMap(FiniteElementSpace * fes , MeshPartition * partition, int nrlayers);
+   ~DofMap(){};
+};
+
+
+
+class LocalDofMap // Constructs dof mapbetween two partitions
+{
+   const FiniteElementCollection *fec=nullptr;
+   MeshPartition * part1=nullptr;
+   MeshPartition * part2=nullptr;
+public:
+   int nrpatch, nx, ny, nz;
+   vector<Array<int>> map1;
+   vector<Array<int>> map2;
+   // constructor
+   LocalDofMap(const FiniteElementCollection * fec_, MeshPartition * part1_, 
+               MeshPartition * part2_);
+   ~LocalDofMap();
+};
+
+
+struct NeighborDofMaps
+{
+private:
+   int dim;
+   MeshPartition * part = nullptr;
+   FiniteElementSpace * fes = nullptr;
+   Mesh * mesh = nullptr;
+   std::vector<std::vector<Array<int>>> OvlpElems;
+   std::vector<std::vector<Array<int>>> OvlpDofMaps;
+
+   DofMap * dmap = nullptr;
+   int nrsubdomains = 0;
+   int ovlp_layers = 0;
+   Array<int> nxyz;
+   void SetElementToOverlap(int ip, int iel, 
+                           const Array<bool> & neg, 
+                           const Array<bool> & pos);
+
+   void MarkOvlpElements();
+   void ComputeNeighborDofMaps();
+
+   void Getijk(int ip, int & i, int & j, int & k) const
+   {
+      k = ip/(nxyz[0]*nxyz[1]);
+      j = (ip-k*nxyz[0]*nxyz[1])/nxyz[0];
+      i = (ip-k*nxyz[0]*nxyz[1])%nxyz[0];
+   }
+
+   int GetPatchId(const Array<int> & ijk) const
+   {
+      int d=ijk.Size();
+      int z = (d==2)? 0 : ijk[2];
+      return part->subdomains(ijk[0],ijk[1],z);
+   }
+      int GetDirectionId(int i, int j, int k=-1) 
+   {
+      int n = 3;
+      return (k+1)*n*n + (j+1)*n + i+1; 
+   }
+   void GetDirections(const int id, int & i, int & j, int & k) 
+   {
+      int n = 3;
+      k = id/(n*n) - 1;
+      j = (id-(k+1)*n*n)/n - 1;
+      i = (id-(k+1)*n*n)%n - 1;
+   }
+
+public:   
+   NeighborDofMaps(MeshPartition * part_, 
+                   FiniteElementSpace * fes_, 
+                   DofMap * dmap_,
+                   int ovlp_layers_);
+
+   void GetNeighborDofMap(const int ip, const Array<int> & directions,
+                          Array<int> & dofmap);
+};

examples/maxwell-solver-dev/common/complex_linalg.cpp

@@ -0,0 +1,358 @@
+
+#include "../../../linalg/kernels.hpp"
+#include "complex_linalg.hpp"
+
+
+ComplexDenseMatrix::ComplexDenseMatrix(){}
+
+ComplexDenseMatrix::ComplexDenseMatrix(int s) 
+{
+   MFEM_ASSERT(s >= 0, "invalid ComplexDenseMatrix size: " << s);
+   height = s;
+   width  = s;
+   if (s > 0)
+   {
+      data = new complex<double>[s*s];
+      *this = 0.0; // init with zeroes
+   }
+}
+
+
+ComplexDenseMatrix::ComplexDenseMatrix(int m, int n) 
+{
+   MFEM_VERIFY(m >= 0 && n >= 0,
+               "invalid DenseMatrix size: " << m << " x " << n);
+   const int s = m*n;
+   height = m;
+   width  = n;
+   if (s > 0)
+   {
+      data = new complex<double>[s];
+      *this = 0.0; // init with zeroes
+   }
+}
+
+void ComplexDenseMatrix::SetSize(int h, int w)
+{
+   MFEM_VERIFY(h >= 0 && w >= 0,
+               "invalid ComplexDenseMatrix size: " << h << " x " << w);
+   if (Height() == h && Width() == w)
+   {
+      return;
+   }
+   height = h;
+   width = w;
+   const int hw = h*w;
+   delete data;
+   data = new complex<double>[hw];
+   *this = 0.0; // init with zeroes
+}
+
+ComplexDenseMatrix &ComplexDenseMatrix::operator=(double c)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] = c;
+   }
+   return *this;
+}
+
+ComplexDenseMatrix &ComplexDenseMatrix::operator=(complex<double> c)
+{
+   const int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] = c;
+   }
+   return *this;
+}
+
+
+std::complex<double> ComplexDenseMatrix::Det() const
+{
+   MFEM_ASSERT(Height() == Width() && Height() > 0,
+               "The matrix must be square and "
+               << "sized larger than zero to compute the determinant."
+               << "  Height() = " << Height()
+               << ", Width() = " << Width());
+
+   switch (Height())
+   {
+      case 1:
+         return data[0];
+
+      case 2:
+         return data[0] * data[3] - data[1] * data[2];
+
+      case 3:
+      {
+         const complex<double> *d = data;
+         return
+            d[0] * (d[4] * d[8] - d[5] * d[7]) +
+            d[3] * (d[2] * d[7] - d[1] * d[8]) +
+            d[6] * (d[1] * d[5] - d[2] * d[4]);
+      }
+      default:
+      {
+         MFEM_ABORT("dim>3 not supported yet");         
+         return 0;
+      }
+   }
+}
+
+DenseMatrix * ComplexDenseMatrix::real() const
+{
+   DenseMatrix * Ar = new DenseMatrix(height,width);
+   double * data = Ar->Data();
+   complex<double> * zdata = this->data;
+   for (int s = 0; s<height*width; s++)
+   {
+         data[s] = zdata[s].real();
+   }
+   return Ar;
+}
+DenseMatrix * ComplexDenseMatrix::imag() const
+{
+   DenseMatrix * Ai = new DenseMatrix(height,width);
+   double * data = Ai->Data();
+   complex<double> * zdata = this->data;
+   for (int s = 0; s<height*width; s++)
+   {
+         data[s] = zdata[s].imag();
+   }
+   return Ai;
+}
+
+void ComplexDenseMatrix::GetReal(DenseMatrix & Ar) 
+{
+   MFEM_ASSERT(Ar.Height() == height && Ar.Width() == width, "Incompatible dimensions");
+   double * data = Ar.Data();
+   complex<double> * zdata = this->data;
+   for (int s = 0; s<height*width; s++)
+   {
+         data[s] = zdata[s].real();
+   }
+}
+
+void ComplexDenseMatrix::GetImag(DenseMatrix & Ai) 
+{
+   double * data = Ai.Data();
+   complex<double> * zdata = this->data;
+   for (int s = 0; s<height*width; s++)
+   {
+         data[s] = zdata[s].imag();
+   }
+}
+
+ComplexDenseMatrix &ComplexDenseMatrix::operator=(const ComplexDenseMatrix &m)
+{
+   SetSize(m.height, m.width);
+
+   const int hw = height * width;
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] = m.data[i];
+   }
+   return *this;
+}
+
+ComplexDenseMatrix &ComplexDenseMatrix::operator+=(const complex<double> *m)
+{
+   const int hw = Height()*Width();
+   for (int i = 0; i < hw; i++)
+   {
+      data[i] += m[i];
+   }
+   return *this;
+}
+
+ComplexDenseMatrix &ComplexDenseMatrix::operator+=(const ComplexDenseMatrix &m)
+{
+   MFEM_ASSERT(Height() == m.Height() && Width() == m.Width(),
+               "incompatible matrix sizes.");
+   return *this += m.GetData();
+}
+
+ComplexDenseMatrix &ComplexDenseMatrix::operator-=(const ComplexDenseMatrix &m)
+{
+   int s = Height()*Width();
+   complex<double> * mdata = m.GetData();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] -= mdata[s];
+   }
+   return *this;
+}
+
+ComplexDenseMatrix &ComplexDenseMatrix::operator*=(complex<double> c)
+{
+   int s = Height()*Width();
+   for (int i = 0; i < s; i++)
+   {
+      data[i] *= c;
+   }
+   return *this;
+}
+
+
+void ComplexDenseMatrix::Print(std::ostream &out, int width_) const
+{
+   // save current output flags
+   ios::fmtflags old_flags = out.flags();
+   // output flags = scientific + show sign
+   out << setiosflags(ios::scientific | ios::showpos);
+   for (int i = 0; i < height; i++)
+   {
+      out << "[row " << i << "]\n";
+      for (int j = 0; j < width; j++)
+      {
+         out << (*this)(i,j);
+         if (j+1 == width || (j+1) % width_ == 0)
+         {
+            out << '\n';
+         }
+         else
+         {
+            out << ' ';
+         }
+      }
+   }
+   // reset output flags to original values
+   out.flags(old_flags);
+}
+
+void ComplexDenseMatrix::PrintMatlab(std::ostream &out) const
+{
+   // save current output flags
+   // ios::fmtflags old_flags = out.flags();
+   // output flags = scientific + show sign
+   // out << setiosflags(ios::scientific | ios::showpos);
+   for (int i = 0; i < height; i++)
+   {
+      for (int j = 0; j < width; j++)
+      {
+         out << (*this)(i,j);
+         out << ' ';
+      }
+      out << "\n";
+   }
+   // reset output flags to original values
+   // out.flags(old_flags);
+}
+
+ComplexDenseMatrixInverse::ComplexDenseMatrixInverse(const ComplexDenseMatrix & A) : ComplexDenseMatrix(A.Height())
+{
+   MFEM_VERIFY(A.Height() == A.Width(), "The matrix is not square");
+   MFEM_VERIFY(A.Height() < 4, "dim > 3 is not supported yet");
+
+   std::complex<double> detA = A.Det();
+   MFEM_VERIFY(abs(A.Det())>1e-14, "The given matrix is singular");
+
+   std::complex<double> * d = this->Data();
+   std::complex<double> *dA = A.GetData();
+   switch (A.Height())
+   {
+   case 1:
+      d[0] = 1.0/dA[0];
+      break;
+   case 2:
+      d[0] =  1.0/detA * dA[3];
+      d[1] = -1.0/detA * dA[1];
+      d[2] = -1.0/detA * dA[2];
+      d[3] =  1.0/detA * dA[0];
+      break;   
+   case 3:
+      d[0] =   1.0/detA*(dA[4]*dA[8] - dA[5]*dA[7]);
+      d[1] =  -1.0/detA*(dA[1]*dA[8] - dA[2]*dA[7]);
+      d[2] =   1.0/detA*(dA[1]*dA[5] - dA[2]*dA[4]);
+      d[3] =  -1.0/detA*(dA[3]*dA[8] - dA[5]*dA[6]);
+      d[4] =   1.0/detA*(dA[0]*dA[8] - dA[2]*dA[6]);
+      d[5] =  -1.0/detA*(dA[0]*dA[5] - dA[2]*dA[3]);
+      d[6] =   1.0/detA*(dA[3]*dA[7] - dA[4]*dA[6]);
+      d[7] =  -1.0/detA*(dA[0]*dA[7] - dA[1]*dA[6]);
+      d[8] =   1.0/detA*(dA[0]*dA[4] - dA[1]*dA[3]);
+      break;      
+   default:
+      // Should be unreachable 
+      break;
+   }
+}
+
+
+
+/// Matrix matrix multiplication.  A = B * C.
+void Mult(const ComplexDenseMatrix &b, const ComplexDenseMatrix &c, ComplexDenseMatrix &a)
+{
+   MFEM_ASSERT(a.Height() == b.Height() && a.Width() == c.Width() &&
+               b.Width() == c.Height(), "incompatible dimensions");
+
+   const int ah = a.Height();
+   const int aw = a.Width();
+   const int bw = b.Width();
+   complex<double> *ad = a.Data();
+   const complex<double> *bd = b.Data();
+   const complex<double> *cd = c.Data();
+   kernels::Mult(ah,aw,bw,bd,cd,ad);               
+}
+
+/// Multiply the transpose of a matrix A with a matrix B:   At*B
+void MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B, ComplexDenseMatrix &AtB)
+{
+   MFEM_ASSERT(A.Width() == AtB.Height() && B.Width() == AtB.Width() &&
+               A.Height() == B.Height(), "incompatible dimensions");
+   const int ah = A.Height();
+   const int aw = A.Width();
+   const int bw = B.Width();
+   const complex<double> *ad = A.Data();
+   const complex<double> *bd = B.Data();
+   complex<double> *cd = AtB.Data();
+
+   for (int j = 0; j < bw; j++)
+   {
+      const complex<double> *ap = ad;
+      for (int i = 0; i < aw; i++)
+      {
+         complex<double> d = 0.0;
+         for (int k = 0; k < ah; k++)
+         {
+            d += ap[k] * bd[k];
+         }
+         *(cd++) = d;
+         ap += ah;
+      }
+      bd += ah;
+   }            
+}
+
+/// Multiply the conjugate transpose of a matrix A with a matrix B:   At*B
+void MultAhB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B, ComplexDenseMatrix &AtB)
+{
+   MFEM_ASSERT(A.Width() == AtB.Height() && B.Width() == AtB.Width() &&
+               A.Height() == B.Height(), "incompatible dimensions");
+      MFEM_ASSERT(A.Width() == AtB.Height() && B.Width() == AtB.Width() &&
+               A.Height() == B.Height(), "incompatible dimensions");
+   const int ah = A.Height();
+   const int aw = A.Width();
+   const int bw = B.Width();
+   const complex<double> *ad = A.Data();
+   const complex<double> *bd = B.Data();
+   complex<double> *cd = AtB.Data();
+
+   for (int j = 0; j < bw; j++)
+   {
+      const complex<double> *ap = ad;
+      for (int i = 0; i < aw; i++)
+      {
+         complex<double> d = 0.0;
+         for (int k = 0; k < ah; k++)
+         {
+            d += conj(ap[k]) * bd[k];
+         }
+         *(cd++) = d;
+         ap += ah;
+      }
+      bd += ah;
+   }                    
+
+}

examples/maxwell-solver-dev/common/complex_linalg.hpp

@@ -0,0 +1,98 @@
+#pragma once
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ComplexDenseMatrix 
+{
+private:
+   std::complex<double> * data = nullptr;
+   int height = 0;
+   int width  = 0;
+public:
+   ComplexDenseMatrix();
+   
+   /// Creates square matrix of size s.
+   explicit ComplexDenseMatrix(int s);
+
+   /// Creates rectangular matrix of size m x n.
+   ComplexDenseMatrix(int m, int n);
+
+   /// Change the size of the DenseMatrix to s x s.
+   void SetSize(int s) { SetSize(s, s); }
+
+   /// Change the size of the DenseMatrix to h x w.
+   void SetSize(int h, int w);
+
+   /// Returns the matrix data array.
+   inline complex<double> *Data() const
+   { return const_cast<complex<double>*>((const complex<double>*)data);}
+
+   /// Returns the matrix data array.
+   inline complex<double> *GetData() const { return Data(); }
+
+   /// Returns reference to a_{ij}.
+   inline complex<double> &operator()(int i, int j);
+   inline const complex<double> &operator()(int i, int j) const;
+
+   inline int Height() const { return height; }
+   inline int Width() const { return width; }
+
+   /// Sets the matrix elements equal to constant c
+   ComplexDenseMatrix &operator=(std::complex<double> c);
+   ComplexDenseMatrix &operator=(double c);
+
+   /// Sets the matrix size and elements equal to those of m
+   ComplexDenseMatrix &operator=(const ComplexDenseMatrix &m);
+   ComplexDenseMatrix &operator+=(const complex<double> *m);
+   ComplexDenseMatrix &operator+=(const ComplexDenseMatrix &m);
+   ComplexDenseMatrix &operator-=(const ComplexDenseMatrix &m);
+   ComplexDenseMatrix &operator*=(complex<double> c);
+
+   /// Calculates the determinant of the matrix
+   /// (for 2x2, 3x3)
+   std::complex<double> Det() const;
+
+   virtual void Print(std::ostream &out = mfem::out, int width_ = 4) const;
+   virtual void PrintMatlab(std::ostream &out = mfem::out) const;
+
+   DenseMatrix * real() const;
+   DenseMatrix * imag() const;
+
+   void GetReal(DenseMatrix & Ar);
+   void GetImag(DenseMatrix & Ai);
+};
+
+inline complex<double> &ComplexDenseMatrix::operator()(int i, int j)
+{
+   MFEM_VERIFY(data && i >= 0 && i < height && j >= 0 && j < width, "");
+   // return data[i*width+j];
+   return data[j*height+i];
+}
+
+inline const complex<double> &ComplexDenseMatrix::operator()(int i, int j) const
+{
+   MFEM_VERIFY(data && i >= 0 && i < height && j >= 0 && j < width, "");
+   // return data[i*width+j];
+   return data[j*height+i];
+}
+
+
+class ComplexDenseMatrixInverse : public ComplexDenseMatrix
+{
+private:
+public:
+   ComplexDenseMatrixInverse(const ComplexDenseMatrix & );   
+};
+
+/// Matrix matrix multiplication.  A = B * C.
+void Mult(const ComplexDenseMatrix &b, const ComplexDenseMatrix &c, ComplexDenseMatrix &a);
+
+/// Multiply the transpose of a matrix A with a matrix B:   At*B
+void MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B, ComplexDenseMatrix &AtB);
+
+/// Multiply the conjugate transpose of a matrix A with a matrix B:   At*B
+void MultAhB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B, ComplexDenseMatrix &AtB);

examples/maxwell-solver-dev/helmholtzp.cpp

@@ -0,0 +1,404 @@
+//
+// Compile with: make helmholtzp
+//
+// Sample runs:   mpirun -np 4 ./helmholtzp -nd 2 -nx 4 -ny 4 -sr 3 -pr 3 -k 16.0 -o 2
+//                mpirun -np 4 ./helmholtzp -nd 3 -nx 2 -ny 2 -nz 2 -sr 3 -pr 1 -k 2.0 -o 2
+//
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "ParDST/ParDST.hpp"
+
+using namespace std;
+using namespace mfem;
+
+// Exact solution and r.h.s., see below for implementation.
+double f_exact_Re(const Vector &x);
+double f_exact_Im(const Vector &x);
+
+double wavespeed(const Vector &x);
+
+double funccoeff_re(const Vector & x);
+double funccoeff_im(const Vector & x);
+
+
+int dim;
+double omega;
+int sol = 1;
+double length = 1.0;
+double pml_length = 0.25;
+Array2D<double>comp_bdr;
+
+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.
+   // finite element order of approximation
+   int order = 1;
+   bool visualization = 1;
+   // number of wavelengths
+   double k = 0.5;
+   // number of serial refinements
+   int ser_ref_levels = 1;
+   // number of parallel refinements
+   int par_ref_levels = 2;
+   // dimension
+   int nd = 2;
+   int nx=2;
+   int ny=2;
+   int nz=2;
+   bool herm_conv = true;
+
+   // optional command line inputs
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&nd, "-nd", "--dim","Problem space dimension");
+   args.AddOption(&nx, "-nx", "--nx","Number of subdomains in x direction");
+   args.AddOption(&ny, "-ny", "--ny","Number of subdomains in y direction");
+   args.AddOption(&nz, "-nz", "--nz","Number of subdomains in z direction");
+   args.AddOption(&sol, "-sol", "--exact",
+                  "Exact solution flag - 0:polynomial, 1: plane wave, -1: unknown exact");
+   args.AddOption(&k, "-k", "--wavelengths",
+                  "Number of wavelengths.");
+   args.AddOption(&pml_length, "-pml_length", "--pml_length",
+                  "Length of the PML region in each direction");
+   args.AddOption(&length, "-length", "--length",
+                  "length of the domain in each direction.");
+   args.AddOption(&ser_ref_levels, "-sr", "--ser_ref_levels",
+                  "Number of Serial Refinements.");
+   args.AddOption(&par_ref_levels, "-pr", "--par_ref_levels",
+                  "Number of Parallel Refinements.");                  
+   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
+                  "--no-hermitian", "Use convention for Hermitian operators.");                  
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   // check if the inputs are correct
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   // Angular frequency
+   omega = 2.0 * M_PI * k;
+
+   // 3. Read the mesh from the given mesh file.
+   Mesh *mesh;
+
+   if (nd == 2)
+   {
+      mesh = new Mesh(1, 1, Element::QUADRILATERAL, true, length, length, false);
+   }
+   else
+   {
+      mesh = new Mesh(1, 1, 1, Element::HEXAHEDRON, true, length, length, length,false);
+   }
+
+   // 3. Executing uniform h-refinement
+   dim = mesh->Dimension();
+   for (int i = 0; i < ser_ref_levels; i++ )
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define a parallel mesh by a partitioning of the serial mesh.
+   // ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   int nprocs;
+   int nprocsx;
+   int nprocsy;
+   int nprocsz;
+   if (dim == 2)
+   {
+      nprocs = sqrt(num_procs); 
+      // nprocsx = nprocs;
+      // nprocsy = nprocs;
+      nprocsx = 1;
+      nprocsy = num_procs;
+      nprocsz = 1;
+   }    
+   else
+   {
+      nprocs = cbrt(num_procs); 
+      // nprocsx = nprocs;
+      // nprocsy = nprocs;
+      // nprocsz = nprocs;
+      nprocsx = 1;
+      if (nz != 1)
+      {
+         nprocsy = sqrt(num_procs);
+         nprocsz = nprocsy;
+      }
+      else
+      {
+         nprocsy = num_procs;
+         nprocsz = 1;
+      }
+   }
+   // MFEM_VERIFY(nprocs*nprocs == num_procs, "Check MPI partitioning");
+   // int nxyz[3] = {num_procs,1,1};
+   // int nxyz[3] = {nprocs,nprocs,1};
+   // int nxyz[3] = {1,num_procs,1};
+   
+   int nxyz[3] = {nprocsx,nprocsy,nprocsz};
+   // int nxyz[3] = {num_procs,1,1};
+   int * part = mesh->CartesianPartitioning(nxyz);
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD,*mesh,part);
+   // ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD,*mesh);
+   delete [] part;
+   delete mesh;
+
+   for (int l = 0; l < par_ref_levels; l++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+
+   double hl = GetUniformMeshElementSize(pmesh);
+   int nrlayers = 3;
+
+   Array2D<double> lengths(dim,2);
+   lengths = hl*nrlayers;
+   // lengths[0][1] = 0.0;
+   // lengths[1][1] = 0.0;
+   // lengths[1][0] = 0.0;
+   // lengths[0][0] = 0.0;
+   CartesianPML pml(pmesh,lengths);
+   pml.SetOmega(omega);
+   comp_bdr.SetSize(dim,2);
+   comp_bdr = pml.GetCompDomainBdr(); 
+
+   // 6. Define a finite element space on the mesh.
+   FiniteElementCollection *fec = new H1_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+   // 6. Set up the linear form (Real and Imaginary part)
+   FunctionCoefficient f_Re(f_exact_Re);
+   FunctionCoefficient f_Im(f_exact_Im);
+
+   // 8. Setup Complex Operator convention
+   ComplexOperator::Convention conv =
+      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
+
+   // ParLinearForm *b_Re(new ParLinearForm);
+   ParComplexLinearForm b(fespace, conv);
+   b.AddDomainIntegrator(new DomainLFIntegrator(f_Re),
+                         new DomainLFIntegrator(f_Im));
+   b.real().Vector::operator=(0.0);
+   b.imag().Vector::operator=(0.0);
+   b.Assemble();
+
+   // 7. Set up the bilinear form (Real and Imaginary part)
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient sigma(-pow(omega, 2));
+
+   FunctionCoefficient ws(wavespeed);
+
+   PmlMatrixCoefficient c1_re(dim,pml_detJ_JT_J_inv_Re,&pml);
+   PmlMatrixCoefficient c1_im(dim,pml_detJ_JT_J_inv_Im,&pml);
+
+   PmlCoefficient detJ_re(pml_detJ_Re,&pml);
+   PmlCoefficient detJ_im(pml_detJ_Im,&pml);
+
+   ProductCoefficient c2_re0(sigma, detJ_re);
+   ProductCoefficient c2_im0(sigma, detJ_im);
+
+   ProductCoefficient c2_re(c2_re0, ws);
+   ProductCoefficient c2_im(c2_im0, ws);
+
+   ParSesquilinearForm a(fespace,conv);
+   a.AddDomainIntegrator(new DiffusionIntegrator(c1_re),
+                         new DiffusionIntegrator(c1_im));
+   a.AddDomainIntegrator(new MassIntegrator(c2_re),
+                         new MassIntegrator(c2_im));
+   a.Assemble();
+   a.Finalize();
+
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+   ess_bdr = 1;
+   fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+   // Solution grid function
+   ParComplexGridFunction p_gf(fespace); p_gf = 0.0;
+   OperatorHandle Ah;
+   Vector X, B;
+
+   a.FormLinearSystem(ess_tdof_list, p_gf, b, Ah, X, B);
+   {
+      StopWatch chrono;
+      chrono.Clear();
+      chrono.Start();
+      ParDST S(&a,lengths,omega, &ws,nrlayers,nx,ny,nz);
+      chrono.Stop();
+      double t1 = chrono.RealTime();
+
+      chrono.Clear();
+      chrono.Start();
+      // X = 0.0;
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPreconditioner(S);
+      gmres.SetOperator(*Ah);
+      gmres.SetRelTol(1e-6);
+      gmres.SetMaxIter(20);
+      gmres.SetPrintLevel(1);
+      gmres.Mult(B, X);
+      chrono.Stop();
+
+      double t2 = chrono.RealTime();
+      
+      MPI_Barrier(MPI_COMM_WORLD);
+
+
+      cout << " myid: " << myid 
+           << ", setup time: " << t1
+           << ", solution time: " << t2 << endl; 
+
+
+      a.RecoverFEMSolution(X,B,p_gf);
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         string keys;
+         if (dim ==2 )
+         {
+            keys = "keys mrRljc\n";
+         }
+         else
+         {
+            keys = "keys mc\n";
+         }
+         socketstream sol_sock_re(vishost, visport);
+         sol_sock_re.precision(8);
+         sol_sock_re << "parallel " << num_procs << " " << myid << "\n"
+                     << "solution\n" << *pmesh << p_gf.real() << keys 
+                     << "window_title 'Numerical Pressure: Real Part' " << flush;                     
+
+         socketstream sol_sock_im(vishost, visport);
+         sol_sock_im.precision(8);
+         sol_sock_im << "parallel " << num_procs << " " << myid << "\n"
+                     << "solution\n" << *pmesh << p_gf.imag() << keys 
+                     << "window_title 'Numerical Pressure: Imag Part' " << flush;                     
+      }
+   }
+
+   delete fespace;
+   delete fec;
+	delete pmesh;
+   MPI_Finalize();
+   return 0;
+}
+
+
+double f_exact_Re(const Vector &x)
+{
+   
+   int nrsources = (dim == 2) ? 4 : 8;
+   Vector x0(nrsources);
+   Vector y0(nrsources);
+   Vector z0(nrsources);
+   x0(0) = 0.25; y0(0) = 0.25; z0(0) = 0.25;
+   x0(1) = 0.75; y0(1) = 0.25; z0(1) = 0.25;
+   x0(2) = 0.25; y0(2) = 0.75; z0(2) = 0.25;
+   x0(3) = 0.75; y0(3) = 0.75; z0(3) = 0.25;
+   if (dim == 3)
+   {
+      x0(4) = 0.25; y0(4) = 0.25; z0(4) = 0.75;
+      x0(5) = 0.75; y0(5) = 0.25; z0(5) = 0.75;
+      x0(6) = 0.25; y0(6) = 0.75; z0(6) = 0.75;
+      x0(7) = 0.75; y0(7) = 0.75; z0(7) = 0.75;
+   }
+  
+   double n = 4.0*omega/M_PI;
+   double coeff = 16.0*omega*omega/M_PI/M_PI/M_PI;
+
+   double f_re = 0.0;
+   // for (int i = 0; i<1; i++)
+   for (int i = 0; i<nrsources; i++)
+   {
+      double beta = pow(x0(i)-x(0),2) + pow(y0(i)-x(1),2);
+      if (dim == 3) { beta += pow(z0(i)-x(2),2); }
+      double alpha = -pow(n,2) * beta;
+      f_re += coeff*exp(alpha);
+   }
+
+   bool in_pml = false;
+   for (int i = 0; i<dim; i++)
+   {
+      if (x(i)<=comp_bdr(i,0) || x(i)>=comp_bdr(i,1))
+      {
+         in_pml = true;
+         break;
+      }
+   }
+   if (in_pml) f_re = 0.0;
+
+   return f_re;
+
+}
+double f_exact_Im(const Vector &x)
+{
+   double f_im;
+   f_im = 0.0;
+   return f_im;
+}
+
+double wavespeed(const Vector &x)
+{
+   double ws;
+   ws = 1.0;
+   return ws;
+}
+
+double funccoeff_re(const Vector & x)
+{
+   return sin(3*M_PI*(x.Sum()));
+}
+
+double funccoeff_im(const Vector & x)
+{
+   return cos(10*M_PI*(x.Sum()));
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

examples/maxwell-solver-dev/makefile

@@ -0,0 +1,67 @@
+# 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 = $(if $(MFEM_DIR:../../..=),$(MFEM_DIR)/examples/maxwell-solver-dev/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = 
+PAR_EXAMPLES = helmholtzp maxwellp
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean
+.PRECIOUS: %.o
+
+COMMON_O= common/PML.o common/MeshPartition.o \
+			 common/Utilities.o common/complex_linalg.o\
+			 ParDST/ParDST.o ParDST/DofMapsDST.o
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(EXAMPLES)
+
+# Rules for building the EXAMPLES
+
+%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(COMMON_O) $(MFEM_LIBS)
+
+# Rules for compiling miniapp dependencies
+$(COMMON_O) $($(EXAMPLES)): \
+%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -f DST/*.o
+	rm -f ParDST/*.o
+	rm -f common/*.o
+	rm -f DST2D/*.o
+	rm -rf *.dSYM *.TVD.*breakpoints
+	rm output/*
+
+

examples/maxwell-solver-dev/maxwellp.cpp

@@ -0,0 +1,535 @@
+//
+// Compile with: make maxwellp
+//
+// Sample runs:   mpirun -np 4 ./maxwellp -nd 2 -nx 4 -ny 4 -sr 3 -pr 3 -k 16.0 -o 2
+//                mpirun -np 4 ./maxwellp -nd 3 -nx 2 -ny 2 -nz 2 -sr 3 -pr 1 -k 2.0 -o 2
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "ParDST/ParDST.hpp"
+
+using namespace std;
+using namespace mfem;
+  
+void source_re(const Vector &x, Vector & f);
+void source_im(const Vector &x, Vector & f);
+void exact_re(const Vector & x, Vector & E);
+void exact_im(const Vector & x, Vector & E);
+void maxwell_solution(const Vector & x, double E[], double curl2E[]);
+double wavespeed(const Vector &x);
+void Mwavespeed(const Vector & x, DenseMatrix & M);
+
+void ess_data_func(const Vector & x, Vector & E);
+
+
+double mu = 1.0;
+double epsilon = 1.0;
+double omega;
+int dim;
+double length = 1.0;
+Array2D<double> comp_bdr;
+Array2D<double> domain_bdr;
+bool exact_known = false;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+   int order = 1;
+   // number of serial refinements
+   int ser_ref_levels = 1;
+   // number of parallel refinements
+   int par_ref_levels = 2;
+   // number of wavelengths
+   double k = 5.0; // 
+   bool herm_conv = true;
+   bool visualization = 1;
+   int nd=2;
+   int nx=2;
+   int ny=2;
+   int nz=2;
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&nd, "-nd", "--dim","Problem space dimension");
+   args.AddOption(&nx, "-nx", "--nx","Number of subdomains in x direction");
+   args.AddOption(&ny, "-ny", "--ny","Number of subdomains in y direction");
+   args.AddOption(&nz, "-nz", "--nz","Number of subdomains in z direction");               
+   args.AddOption(&ser_ref_levels, "-sr", "--ser_ref_levels",
+                  "Number of Serial Refinements.");
+   args.AddOption(&par_ref_levels, "-pr", "--par_ref_levels",
+                  "Number of Parallel Refinements."); 
+   args.AddOption(&mu, "-mu", "--permeability",
+                  "Permeability of free space (or 1/(spring constant)).");
+   args.AddOption(&epsilon, "-eps", "--permittivity",
+                  "Permittivity of free space (or mass constant).");
+   args.AddOption(&k, "-k", "--wavelengths",
+                  "Number of wavelengths.");
+   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
+                  "--no-hermitian", "Use convention for Hermitian operators.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+
+   // check if the inputs are correct
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // Angular frequency
+   omega = 2.0 * M_PI * k;
+
+   Mesh *mesh;
+
+
+   int nel = 1;
+   if (nd == 2)
+   {
+      mesh = new Mesh(nel, nel, Element::QUADRILATERAL, true, length, length, false);
+   }
+   else
+   {
+      mesh = new Mesh(nel, nel, nel, Element::HEXAHEDRON, true, length, length, length,false);
+   }
+
+   dim = mesh->Dimension();
+
+   // 4. Refine the mesh to increase the resolution.
+   for (int l = 0; l < ser_ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define a parallel mesh by a partitioning of the serial mesh.
+   // ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   int nprocs;
+   int nprocsx;
+   int nprocsy;
+   int nprocsz;
+   if (dim == 2)
+   {
+      nprocs = sqrt(num_procs); 
+      nprocsx = nprocs;
+      nprocsy = nprocs;
+      nprocsz = 1;
+   }    
+   else
+   {
+      nprocs = cbrt(num_procs); 
+      nprocsx = nprocs;
+      nprocsy = nprocs;
+      nprocsz = nprocs;
+   }
+   int nxyz[3] = {nprocsx,nprocsy,nprocsz};
+   int * part = mesh->CartesianPartitioning(nxyz);
+   // ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD,*mesh,part);
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD,*mesh);
+   delete [] part;
+
+   
+   delete mesh;
+
+   for (int l = 0; l < par_ref_levels; l++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+
+   double hl = GetUniformMeshElementSize(pmesh);
+   int nrlayers = 3;
+   Array2D<double> lengths(dim,2);
+   lengths = hl*nrlayers;
+   // lengths[0][1] = 0.0;
+   // lengths[1][1] = 0.0;
+   // lengths[1][0] = 0.0;
+   // lengths[0][0] = 0.0;
+   if (exact_known) lengths = 0.0;
+   // CartesianPML pml(mesh,lengths);
+   CartesianPML pml(pmesh,lengths);
+   pml.SetOmega(omega);
+   comp_bdr.SetSize(dim,2);
+   comp_bdr = pml.GetCompDomainBdr(); 
+
+
+   // 6. Define a finite element space on the mesh. Here we use the Nedelec
+   //    finite elements of the specified order.
+   FiniteElementCollection *fec = new ND_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 7. Determine the list of true essential boundary dofs. In this example,
+   //    the boundary conditions are defined based on the specific mesh and the
+   //    problem type.
+   Array<int> ess_tdof_list;
+   if (pmesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 8. Setup Complex Operator convention
+   ComplexOperator::Convention conv =
+      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
+
+   // 9. Set up the linear form b(.) which corresponds to the right-hand side of
+   //    the FEM linear system.
+   VectorFunctionCoefficient f_re(dim, source_re);
+   VectorFunctionCoefficient f_im(dim, source_re);
+   ParComplexLinearForm b(fespace, conv);
+   b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f_re),
+                         new VectorFEDomainLFIntegrator(f_im));
+   b.Vector::operator=(0.0);
+   b.Assemble();
+
+   // 10. Define the solution vector x as a complex finite element grid function
+   //     corresponding to fespace.
+   ParComplexGridFunction x(fespace);
+   x = 0.0;
+   // VectorFunctionCoefficient done(dim,ess_data_func);
+   // x.ProjectCoefficient(done,done);
+   VectorFunctionCoefficient E_re(dim,exact_re);
+   VectorFunctionCoefficient E_im(dim,exact_re);
+   if (exact_known)
+   {
+      x.ProjectCoefficient(E_re,E_re);
+   }
+   // 11. Set up the sesquilinear form a(.,.)
+   //
+   //       1/mu (1/det(J) J^T J Curl E, Curl F)
+   //        - omega^2 * epsilon (det(J) * (J^T J)^-1 * E, F)
+   //
+   FunctionCoefficient ws(wavespeed);
+
+   // MatrixFunctionCoefficient Mws(dim,Mwavespeed);
+
+   // DenseMatrix M(dim); M = 0.0;
+   // M(0,0) = -pow(omega, 2);
+   // M(1,1) = -pow(omega, 2);
+   // M(2,2) = -pow(omega, 2);
+   // MatrixConstantCoefficient Momeg(M);
+   MatrixFunctionCoefficient eps_func(dim,Mwavespeed);
+
+   ConstantCoefficient omeg(-pow(omega, 2));
+   int cdim = (dim == 2) ? 1 : dim;
+   PmlMatrixCoefficient pml_c1_Re(cdim,detJ_inv_JT_J_Re, &pml);
+   PmlMatrixCoefficient pml_c1_Im(cdim,detJ_inv_JT_J_Im, &pml);
+
+   PmlMatrixCoefficient pml_c2_Re(dim, detJ_JT_J_inv_Re,&pml);
+   PmlMatrixCoefficient pml_c2_Im(dim, detJ_JT_J_inv_Im,&pml);
+   ScalarMatrixProductCoefficient c2_Re0(omeg,pml_c2_Re);
+   ScalarMatrixProductCoefficient c2_Im0(omeg,pml_c2_Im);
+   // ScalarMatrixProductCoefficient c2_Re(ws,c2_Re0);
+   // ScalarMatrixProductCoefficient c2_Im(ws,c2_Im0);
+
+   MatrixMatrixProductCoefficient c2_Re(c2_Re0,eps_func);
+   MatrixMatrixProductCoefficient c2_Im(c2_Im0,eps_func);
+
+   // MatrixMatrixProductCoefficient c2_Re0(Momeg,pml_c2_Re);
+   // MatrixMatrixProductCoefficient c2_Im0(Momeg,pml_c2_Im);
+   // MatrixMatrixProductCoefficient c2_Re(Mws,c2_Re0);
+   // MatrixMatrixProductCoefficient c2_Im(Mws,c2_Im0);
+
+
+   ParSesquilinearForm a(fespace, conv);
+   a.AddDomainIntegrator(new CurlCurlIntegrator(pml_c1_Re),
+                         new CurlCurlIntegrator(pml_c1_Im));
+   a.AddDomainIntegrator(new VectorFEMassIntegrator(c2_Re),
+                         new VectorFEMassIntegrator(c2_Im));
+
+   a.Assemble(0);
+
+   OperatorHandle Ah;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, x, b, Ah, X, B);
+
+   ComplexSparseMatrix * Ac = Ah.As<ComplexSparseMatrix>();
+   StopWatch chrono;
+
+
+   chrono.Clear();
+   chrono.Start();
+   ParDST * S = new ParDST(&a,lengths, omega, &ws, nrlayers, nx, ny, nz);
+   chrono.Stop();
+   double t1 = chrono.RealTime();
+
+   chrono.Clear();
+   chrono.Start();
+   // X = 0.0;
+	GMRESSolver gmres(MPI_COMM_WORLD);
+	// gmres.iterative_mode = true;
+   gmres.SetPreconditioner(*S);
+	gmres.SetOperator(*Ac);
+	gmres.SetRelTol(1e-8);
+	gmres.SetMaxIter(100);
+	gmres.SetPrintLevel(1);
+	gmres.Mult(B, X);
+   delete S;
+   chrono.Stop();
+   double t2 = chrono.RealTime();
+
+   MPI_Barrier(MPI_COMM_WORLD);
+
+
+   cout << " myid: " << myid 
+         << ", setup time: " << t1
+         << ", solution time: " << t2 << endl; 
+
+   a.RecoverFEMSolution(X, b, x);
+
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      string keys;
+      if (dim ==2 )
+      {
+         keys = "keys mrRljc\n";
+      }
+      else
+      {
+         keys = "keys mc\n";
+      }
+      socketstream sol_sock_re(vishost, visport);
+      sol_sock_re.precision(8);
+      sol_sock_re << "parallel " << num_procs << " " << myid << "\n"
+                  << "solution\n" << *pmesh << x.real() << keys 
+                  << "window_title 'E: Real Part' " << flush;                     
+
+      socketstream sol_sock_im(vishost, visport);
+      sol_sock_im.precision(8);
+      sol_sock_im << "parallel " << num_procs << " " << myid << "\n"
+                  << "solution\n" << *pmesh << x.imag() << keys 
+                  << "window_title 'E: Imag Part' " << flush;   
+
+      {
+         ParGridFunction x_t(fespace);
+         x_t = x.real();
+
+         socketstream sol_sock(vishost, visport);
+         sol_sock.precision(8);
+         sol_sock << "parallel " << num_procs << " " << myid << "\n"
+                  << "solution\n" << *pmesh << x_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 num_frames = 32;
+         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), x.real(), sin(2.0*M_PI*t), x.imag(), x_t);
+            sol_sock << "parallel " << num_procs << " " << myid << "\n";
+            sol_sock << "solution\n" << *pmesh << x_t
+                     << "window_title '" << oss.str() << "'" << flush;
+            i++;
+         }
+      }
+   }
+
+   // 18. Free the used memory.
+   delete fespace;
+   delete fec;
+   delete pmesh;
+   MPI_Finalize();
+   return 0;
+}
+
+void source_re(const Vector &x, Vector &f)
+{
+   f = 0.0;
+   if (exact_known)
+   {
+      double E[3], curl2E[3];
+      maxwell_solution(x, E, curl2E);
+      // curl ( curl E) +/- omega^2 E = f
+      double coeff = -omega * omega;
+      f(0) = curl2E[0] + coeff * E[0];
+      f(1) = curl2E[1] + coeff * E[1];
+      if (dim == 2)
+      {
+         if (x.Size() == 3) {f(2)=0.0;}
+      }
+      else
+      {
+         f(2) = curl2E[2] + coeff * E[2];
+      }
+   }
+   else
+   {
+      int nrsources = (dim == 2) ? 4 : 8;
+      Vector x0(nrsources);
+      Vector y0(nrsources);
+      Vector z0(nrsources);
+      x0(0) = 0.25; y0(0) = 0.25; z0(0) = 0.25;
+      x0(1) = 0.75; y0(1) = 0.25; z0(1) = 0.25;
+      x0(2) = 0.25; y0(2) = 0.75; z0(2) = 0.25;
+      x0(3) = 0.75; y0(3) = 0.75; z0(3) = 0.25;
+      if (dim == 3)
+      {
+         x0(4) = 0.25; y0(4) = 0.25; z0(4) = 0.75;
+         x0(5) = 0.75; y0(5) = 0.25; z0(5) = 0.75;
+         x0(6) = 0.25; y0(6) = 0.75; z0(6) = 0.75;
+         x0(7) = 0.75; y0(7) = 0.75; z0(7) = 0.75;
+      }
+  
+      double n = 4.0*omega/M_PI;
+      double coeff = 16.0*omega*omega/M_PI/M_PI/M_PI;
+
+      for (int i = 0; i<nrsources; i++)
+      {
+         double beta = pow(x0(i)-x(0),2) + pow(y0(i)-x(1),2);
+         if (dim == 3) { beta += pow(z0(i)-x(2),2); }
+         double alpha = -pow(n,2) * beta;
+         f[0] += coeff*exp(alpha);
+      }
+
+      bool in_pml = false;
+      for (int i = 0; i<dim; i++)
+      {
+         if (x(i)<=comp_bdr(i,0) || x(i)>=comp_bdr(i,1))
+         {
+            in_pml = true;
+            break;
+         }
+      }
+      if (in_pml) f = 0.0;
+   }
+}
+
+void source_im(const Vector &x, Vector &f)
+{
+   f = 0.0;
+}
+
+double wavespeed(const Vector &x)
+{
+   double ws;
+   ws = 1.0;
+   return ws;
+}
+
+void Mwavespeed(const Vector & x, DenseMatrix & M)
+{
+   M = 0.0;
+   M(0,0) = 1.0;
+   M(1,1) = 1.0;
+   // M(2,2) = 4.0*x(0)-1.0;
+   if (dim == 3) M(2,2) = 1.0;
+}
+
+
+void exact_re(const Vector & x, Vector & E)
+{
+   double curl2E[3];
+   maxwell_solution(x, E, curl2E);
+}
+void exact_im(const Vector & x, Vector & E)
+{
+   // double curl2E[3];
+   // maxwell_solution(x, E, curl2E);
+   E = 0.0;
+}
+void maxwell_solution(const Vector & x, double E[], double curl2E[])
+{
+   // point source
+   if (dim == 2)
+   {
+      // shift to avoid singularity
+      double x0 = x(0) + 0.1;
+      double x1 = x(1) + 0.1;
+      //
+      double r = sqrt(x0 * x0 + x1 * x1);
+
+      E[0] = cos(omega * r);
+      E[1] = 0.0;
+
+      double r_x = x0 / r;
+      double r_y = x1 / r;
+      double r_xy = -(r_x / r) * r_y;
+      double r_yx = r_xy;
+      double r_yy = (1.0 / r) * (1.0 - r_y * r_y);
+
+      curl2E[0] = omega * ((r_yy ) * sin(omega * r) + (omega * r_y * r_y) * cos(omega * r));
+      curl2E[1] = -omega * (r_yx * sin(omega * r) + omega * r_y * r_x * cos(omega * r));
+      curl2E[2] = 0.0;
+   }
+   else
+   {
+   // shift to avoid singularity
+      double x0 = x(0) + 0.1;
+      double x1 = x(1) + 0.1;
+      double x2 = x(2) + 0.1;
+      //
+      double r = sqrt(x0 * x0 + x1 * x1 + x2 * x2);
+
+      E[0] = cos(omega * r);
+      E[1] = 0.0;
+      E[2] = 0.0;
+
+      double r_x = x0 / r;
+      double r_y = x1 / r;
+      double r_z = x2 / r;
+      double r_xy = -(r_x / r) * r_y;
+      double r_xz = -(r_x / r) * r_z;
+      double r_yx = r_xy;
+      double r_yy = (1.0 / r) * (1.0 - r_y * r_y);
+      double r_zx = r_xz;
+      double r_zz = (1.0 / r) * (1.0 - r_z * r_z);
+
+      curl2E[0] = omega * ((r_yy + r_zz) * sin(omega * r) +
+                           (omega * r_y * r_y + omega * r_z * r_z) * cos(omega * r));
+      curl2E[1] = -omega * (r_yx * sin(omega * r) + omega * r_y * r_x * cos(omega * r));
+      curl2E[2] = -omega * (r_zx * sin(omega * r) + omega * r_z * r_x * cos(omega * r));
+   }
+}
+
+
+void ess_data_func(const Vector & x, Vector & E)
+{
+   E = 0.0;
+   if (x(1)==0.0) E[0] = sin(x(0)+x(1));
+
+   bool in_pml = false;
+   for (int i = 0; i<dim; i++)
+   {
+      if (x(i)<comp_bdr(i,0) || x(i)>comp_bdr(i,1))
+      {
+         in_pml = true;
+         break;
+      }
+   }
+   if (in_pml) E = 0.0;
+
+}

fem/coefficient.cpp

@@ -629,6 +629,27 @@ void MatrixVectorProductCoefficient::Eval(Vector &V, ElementTransformation &T,
    ma.Mult(vb, V);
 }
 
+
+MatrixMatrixProductCoefficient::MatrixMatrixProductCoefficient(MatrixCoefficient &A,
+                                           MatrixCoefficient &B)
+   : MatrixCoefficient(A.GetHeight(), A.GetWidth()),
+     a(&A), b(&B),
+     ma(A.GetHeight(), A.GetWidth()), 
+     mb(B.GetHeight(), B.GetWidth())
+{
+   MFEM_ASSERT(A.GetWidth() == B.GetHeight(),
+               "MatrixMatrixProductCoefficient:  "
+               "Arguments must have the same dimensions.");
+}
+
+void MatrixMatrixProductCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
+                                          const IntegrationPoint &ip)
+{
+   a->Eval(ma, T, ip);
+   b->Eval(mb, T, ip);
+   Mult(ma, mb, M);
+}
+
 void IdentityMatrixCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
                                      const IntegrationPoint &ip)
 {

fem/coefficient.hpp

@@ -1470,6 +1470,42 @@ public:
                      const IntegrationPoint &ip);
 };
 
+/** @brief Matrix coefficient defined as a product of two 
+    matrix coefficients */
+class MatrixMatrixProductCoefficient : public MatrixCoefficient
+{
+private:
+   MatrixCoefficient * a;
+   MatrixCoefficient * b;
+
+   mutable DenseMatrix ma;
+   mutable DenseMatrix mb;
+
+public:
+   /// Constructor with two coefficients.  Result is A*B.
+   MatrixMatrixProductCoefficient(MatrixCoefficient &A, MatrixCoefficient &B);
+
+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Reset the vector coefficient
+   void SetBCoef(MatrixCoefficient &B) { b = &B; }
+   /// Return the vector coefficient
+   MatrixCoefficient * GetBCoef() const { return b; }
+
+   /// Evaluate the vector coefficient at @a ip.
+   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
+/// Convenient alias for the MatrixVectorProductCoefficient
+typedef MatrixMatrixProductCoefficient MatMatCoefficient;
+
+
+
+
 /// Matrix coefficient defined as the linear combination of two matrices
 class MatrixSumCoefficient : public MatrixCoefficient
 {