15 #include "../general/forall.hpp" 43 trial_fes(
a->FESpace()),
44 test_fes(
a->FESpace())
54 const int integratorCount = integrators.
Size();
55 for (
int i = 0; i < integratorCount; ++i)
57 integrators[i]->AssembleMF(*
a->
FESpace());
65 const int iSz = integrators.
Size();
69 for (
int i = 0; i < iSz; ++i)
71 integrators[i]->AssembleDiagonalMF(
localY);
88 for (
int i = 0; i < iSz; ++i)
90 integrators[i]->AssembleDiagonalMF(y);
130 const int iSz = integrators.
Size();
135 for (
int i = 0; i < iSz; ++i)
137 integrators[i]->AddMultMF(x, y);
144 for (
int i = 0; i < iSz; ++i)
152 const int iFISz = intFaceIntegrators.
Size();
159 for (
int i = 0; i < iFISz; ++i)
168 const int bFISz = bdrFaceIntegrators.
Size();
175 for (
int i = 0; i < bFISz; ++i)
187 const int iSz = integrators.
Size();
192 for (
int i = 0; i < iSz; ++i)
202 for (
int i = 0; i < iSz; ++i)
204 integrators[i]->AddMultTransposeMF(x, y);
209 const int iFISz = intFaceIntegrators.
Size();
216 for (
int i = 0; i < iFISz; ++i)
225 const int bFISz = bdrFaceIntegrators.
Size();
232 for (
int i = 0; i < bFISz; ++i)
244 trial_fes(
a->FESpace()),
245 test_fes(
a->FESpace())
295 const int integratorCount = integrators.
Size();
296 for (
int i = 0; i < integratorCount; ++i)
298 integrators[i]->AssemblePA(*
a->
FESpace());
301 MFEM_VERIFY(
a->
GetBBFI()->Size() == 0,
302 "Partial assembly does not support AddBoundaryIntegrator yet.");
305 const int intFaceIntegratorCount = intFaceIntegrators.Size();
306 for (
int i = 0; i < intFaceIntegratorCount; ++i)
308 intFaceIntegrators[i]->AssemblePAInteriorFaces(*
a->
FESpace());
312 const int boundFaceIntegratorCount = bdrFaceIntegrators.
Size();
313 for (
int i = 0; i < boundFaceIntegratorCount; ++i)
315 bdrFaceIntegrators[i]->AssemblePABoundaryFaces(*
a->
FESpace());
323 const int iSz = integrators.
Size();
327 for (
int i = 0; i < iSz; ++i)
329 integrators[i]->AssembleDiagonalPA(
localY);
346 for (
int i = 0; i < iSz; ++i)
348 integrators[i]->AssembleDiagonalPA(y);
388 const int iSz = integrators.
Size();
393 for (
int i = 0; i < iSz; ++i)
395 integrators[i]->AddMultPA(x, y);
402 for (
int i = 0; i < iSz; ++i)
410 const int iFISz = intFaceIntegrators.
Size();
417 for (
int i = 0; i < iFISz; ++i)
426 const int bFISz = bdrFaceIntegrators.
Size();
433 for (
int i = 0; i < bFISz; ++i)
445 const int iSz = integrators.
Size();
450 for (
int i = 0; i < iSz; ++i)
460 for (
int i = 0; i < iSz; ++i)
462 integrators[i]->AddMultTransposePA(x, y);
467 const int iFISz = intFaceIntegrators.
Size();
474 for (
int i = 0; i < iFISz; ++i)
483 const int bFISz = bdrFaceIntegrators.
Size();
490 for (
int i = 0; i < bFISz; ++i)
502 factorize_face_terms(false)
521 const int integratorCount = integrators.
Size();
522 if ( integratorCount == 0 )
526 for (
int i = 0; i < integratorCount; ++i)
535 MFEM_VERIFY(
a->
GetBBFI()->Size() == 0,
536 "Element assembly does not support AddBoundaryIntegrator yet.");
539 const int intFaceIntegratorCount = intFaceIntegrators.Size();
540 if (intFaceIntegratorCount>0)
546 for (
int i = 0; i < intFaceIntegratorCount; ++i)
548 intFaceIntegrators[i]->AssembleEAInteriorFaces(*
a->
FESpace(),
555 const int boundFaceIntegratorCount = bdrFaceIntegrators.
Size();
556 if (boundFaceIntegratorCount>0)
562 for (
int i = 0; i < boundFaceIntegratorCount; ++i)
599 MFEM_FORALL(glob_j,
ne*NDOFS,
601 const int e = glob_j/NDOFS;
602 const int j = glob_j%NDOFS;
604 for (
int i = 0; i < NDOFS; i++)
606 res += A(i, j, e)*X(i, e);
619 const int iFISz = intFaceIntegrators.
Size();
634 MFEM_FORALL(glob_j,
nf_int*NDOFS,
636 const int f = glob_j/NDOFS;
637 const int j = glob_j%NDOFS;
639 for (
int i = 0; i < NDOFS; i++)
641 res += A_int(i, j, 0,
f)*X(i, 0,
f);
645 for (
int i = 0; i < NDOFS; i++)
647 res += A_int(i, j, 1,
f)*X(i, 1,
f);
653 MFEM_FORALL(glob_j,
nf_int*NDOFS,
655 const int f = glob_j/NDOFS;
656 const int j = glob_j%NDOFS;
658 for (
int i = 0; i < NDOFS; i++)
660 res += A_ext(i, j, 0,
f)*X(i, 0,
f);
664 for (
int i = 0; i < NDOFS; i++)
666 res += A_ext(i, j, 1,
f)*X(i, 1,
f);
677 const int bFISz = bdrFaceIntegrators.
Size();
690 MFEM_FORALL(glob_j,
nf_bdr*NDOFS,
692 const int f = glob_j/NDOFS;
693 const int j = glob_j%NDOFS;
695 for (
int i = 0; i < NDOFS; i++)
697 res += A(i, j,
f)*X(i,
f);
727 MFEM_FORALL(glob_j,
ne*NDOFS,
729 const int e = glob_j/NDOFS;
730 const int j = glob_j%NDOFS;
732 for (
int i = 0; i < NDOFS; i++)
734 res += A(j, i, e)*X(i, e);
747 const int iFISz = intFaceIntegrators.
Size();
762 MFEM_FORALL(glob_j,
nf_int*NDOFS,
764 const int f = glob_j/NDOFS;
765 const int j = glob_j%NDOFS;
767 for (
int i = 0; i < NDOFS; i++)
769 res += A_int(j, i, 0,
f)*X(i, 0,
f);
773 for (
int i = 0; i < NDOFS; i++)
775 res += A_int(j, i, 1,
f)*X(i, 1,
f);
781 MFEM_FORALL(glob_j,
nf_int*NDOFS,
783 const int f = glob_j/NDOFS;
784 const int j = glob_j%NDOFS;
786 for (
int i = 0; i < NDOFS; i++)
788 res += A_ext(j, i, 1,
f)*X(i, 0,
f);
792 for (
int i = 0; i < NDOFS; i++)
794 res += A_ext(j, i, 0,
f)*X(i, 1,
f);
805 const int bFISz = bdrFaceIntegrators.
Size();
818 MFEM_FORALL(glob_j,
nf_bdr*NDOFS,
820 const int f = glob_j/NDOFS;
821 const int j = glob_j%NDOFS;
823 for (
int i = 0; i < NDOFS; i++)
825 res += A(j, i,
f)*X(i,
f);
856 bool keep_nbr_block =
false;
870 keep_nbr_block =
true;
884 "Full Assembly not yet supported on NCMesh.");
892 for (
int i =
height; i > 0; i--)
917 "Full Assembly not yet supported on NCMesh.");
923 if (restF) { restF->
FillI(*mat, keep_nbr_block); }
927 for (
int i = 0; i <
height; i++)
929 const int nnz = h_I[i];
944 for (
int i =
height; i > 0; i--)
968 if (
auto pa = dynamic_cast<ParBilinearForm*>(
a) )
970 pa->ParallelRAP(*pa->mat, A);
983 "Only DiagonalPolicy::DIAG_ONE supported with" 984 " FABilinearFormExtension.");
986 if ( dynamic_cast<ParBilinearForm*>(
a) )
989 DiagonalPolicy::DIAG_ONE);
995 DiagonalPolicy::DIAG_ONE);
1022 if ( (pfes = dynamic_cast<const ParFiniteElementSpace*>(
test_fes)) )
1026 x_gf.
MakeRef(const_cast<ParFiniteElementSpace*>(pfes),
1027 const_cast<Vector&>(x),0);
1031 auto dg_x_ptr = dg_x.
Write();
1032 auto x_ptr = x.
Read();
1033 MFEM_FORALL(i,local_size,
1035 dg_x_ptr[i] = x_ptr[i];
1037 const int shared_size = shared_x.
Size();
1038 auto shared_x_ptr = shared_x.
Read();
1039 MFEM_FORALL(i,shared_size,
1041 dg_x_ptr[local_size+i] = shared_x_ptr[i];
1044 if ((pform = dynamic_cast<ParBilinearForm*>(
a)) && (pform->keep_nbr_block))
1046 mat->
Mult(dg_x, dg_y);
1048 auto dg_y_ptr = dg_y.
Read();
1050 MFEM_FORALL(i,local_size,
1052 y_ptr[i] += dg_y_ptr[i];
1083 if ( (pfes = dynamic_cast<const ParFiniteElementSpace*>(
test_fes)) )
1087 x_gf.
MakeRef(const_cast<ParFiniteElementSpace*>(pfes),
1088 const_cast<Vector&>(x),0);
1092 auto dg_x_ptr = dg_x.
Write();
1093 auto x_ptr = x.
Read();
1094 MFEM_FORALL(i,local_size,
1096 dg_x_ptr[i] = x_ptr[i];
1098 const int shared_size = shared_x.
Size();
1099 auto shared_x_ptr = shared_x.
Read();
1100 MFEM_FORALL(i,shared_size,
1102 dg_x_ptr[local_size+i] = shared_x_ptr[i];
1105 if ((pb = dynamic_cast<ParBilinearForm*>(
a)) && (pb->keep_nbr_block))
1109 auto dg_y_ptr = dg_y.
Read();
1111 MFEM_FORALL(i,local_size,
1113 y_ptr[i] += dg_y_ptr[i];
1142 :
Operator(form->Height(), form->Width()),
a(form)
1172 trial_fes(form->TrialFESpace()),
1173 test_fes(form->TestFESpace()),
1174 elem_restrict_trial(NULL),
1175 elem_restrict_test(NULL)
1183 const int integratorCount = integrators.
Size();
1184 for (
int i = 0; i < integratorCount; ++i)
1188 MFEM_VERIFY(
a->
GetBBFI()->Size() == 0,
1189 "Partial assembly does not support AddBoundaryIntegrator yet.");
1191 "Partial assembly does not support AddTraceFaceIntegrator yet.");
1193 "Partial assembly does not support AddBdrTraceFaceIntegrator yet.");
1250 const double c)
const 1253 if (elem_restrict_x)
1255 elem_restrict_x->
Mult(x, localX);
1272 if (elem_restrict_y)
1290 const double c)
const 1293 const int iSz = integrators.
Size();
1300 for (
int i = 0; i < iSz; ++i)
1322 const double c)
const 1325 const int iSz = integrators.
Size();
1332 for (
int i = 0; i < iSz; ++i)
1351 const int iSz = integrators.
Size();
1357 if (H1elem_restrict_trial)
1370 for (
int i = 0; i < iSz; ++i)
1378 integrators[i]->AssembleDiagonalPA_ADAt(D,
localTest);
1383 if (H1elem_restrict_test)
1396 for (
int i = 0; i < iSz; ++i)
1400 integrators[i]->AssembleDiagonalPA_ADAt(
localTrial, diag);
1404 integrators[i]->AssembleDiagonalPA_ADAt(D, diag);
1426 const int integratorCount = integrators.
Size();
1427 for (
int i = 0; i < integratorCount; ++i)
1445 mfem_error(
"A real ElementRestriction is required in this setting!");
1448 auto tm = test_multiplicity.
ReadWrite();
1449 MFEM_FORALL(i, test_multiplicity.
Size(),
1451 tm[i] = 1.0 / tm[i];
1459 const int iSz = integrators.
Size();
1466 for (
int i = 0; i < iSz; ++i)
1484 mfem_error(
"In this setting you need a real ElementRestriction!");
1492 const int iSz = integrators.
Size();
1498 MFEM_VERIFY(x.
Size() == test_multiplicity.
Size(),
"Input vector of wrong size");
1500 auto tm = test_multiplicity.
Read();
1501 MFEM_FORALL(i, x.
Size(),
1509 for (
int i = 0; i < iSz; ++i)
1523 mfem_error(
"Trial ElementRestriction not defined");
void FillI(SparseMatrix &mat) const
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
void FormRectangularSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Operator *&A)
Return in A a parallel (on truedofs) version of this rectangular operator (including constraints)...
void Mult(const Vector &x, Vector &y) const override=0
Extract the face degrees of freedom from x into y.
void ExchangeFaceNbrData()
void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const
void SetSize(int s)
Resize the vector to size s.
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Geometry::Type GetFaceGeometry(int i) const
Return the Geometry::Type associated with face i.
Pointer to an Operator of a specified type.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
int Size() const
Returns the size of the vector.
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Abstract parallel finite element space.
int GetFaceNbrVSize() const
virtual void FillJAndData(const Vector &fea_data, SparseMatrix &mat, const bool keep_nbr_block=false) const
Fill the J and Data arrays of the SparseMatrix corresponding to the sparsity pattern given by this L2...
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
void SortColumnIndices()
Sort the column indices corresponding to each row.
Memory< double > & GetMemoryData()
Memory< int > & GetMemoryI()
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Memory< int > & GetMemoryJ()
void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, Operator *&A, Vector &X, Vector &B)
Form a column-constrained linear system using a matrix-free approach.
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
virtual void AddFaceMatricesToElementMatrices(const Vector &fea_data, Vector &ea_data) const
This methods adds the DG face matrices to the element matrices.
void ExchangeFaceNbrData()
double f(const Vector &xvec)
PADiscreteLinearOperatorExtension(DiscreteLinearOperator *linop)
virtual void AddMultTransposeInPlace(Vector &x, Vector &y) const
Add the face degrees of freedom x to the element degrees of freedom y. Perform the same computation a...
Native ordering as defined by the FiniteElement.
Operator that extracts Face degrees of freedom for L2 spaces.
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
void MultLeftInverse(const Vector &x, Vector &y) const
void MultTransposeUnsigned(const Vector &x, Vector &y) const
Compute MultTranspose without applying signs based on DOF orientations.
virtual const Operator * GetOutputRestrictionTranspose() const
Transpose of GetOutputRestriction, directly available in this form to facilitate matrix-free RAP-type...
void FormRectangularSystemOperator(const Array< int > &, const Array< int > &, OperatorHandle &A)
Setup OperatorHandle A to contain constrained linear operator.
virtual void FillI(SparseMatrix &mat, const bool keep_nbr_block=false) const
Fill the I array of SparseMatrix corresponding to the sparsity pattern given by this L2FaceRestrictio...
void Assemble()
Partial assembly of all internal integrators.
const Operator * GetOutputRestrictionTranspose() const
Transpose of GetOutputRestriction, directly available in this form to facilitate matrix-free RAP-type...
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Mesh * GetMesh() const
Returns the mesh.
Operator that converts FiniteElementSpace L-vectors to E-vectors.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
bool DeviceCanUseCeed()
Function that determines if a CEED kernel should be used, based on the current mfem::Device configura...
void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const
void MultUnsigned(const Vector &x, Vector &y) const
Compute Mult without applying signs based on DOF orientations.
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, Operator *&A, Vector &X, Vector &B, int copy_interior=0)
Form a constrained linear system using a matrix-free approach.
int GetDof() const
Returns the number of degrees of freedom in the finite element.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
static bool Allows(unsigned long b_mask)
Return true if any of the backends in the backend mask, b_mask, are allowed.
void FillSparseMatrix(const Vector &mat_ea, SparseMatrix &mat) const
Fill a Sparse Matrix with Element Matrices.
Vector & Set(const double a, const Vector &x)
(*this) = a * x
int height
Dimension of the output / number of rows in the matrix.
int GetNE() const
Returns number of elements.
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Rectangular Operator for imposing essential boundary conditions on the input space using only the act...
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
void SyncAliasMemory(const Vector &v) const
Update the alias memory location of the vector to match v.
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
void FormSystemOperator(const Array< int > &ess_tdof_list, Operator *&A)
Return in A a parallel (on truedofs) version of this square operator.
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
void AddMult(const Vector &x, Vector &y, const double c=1.0) const
y += c*A*x
void OverrideSize(int height_, int width_)
Sets the height and width of the matrix.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Lexicographic ordering for tensor-product FiniteElements.
Operator * SetupRAP(const Operator *Pi, const Operator *Po)
Returns RAP Operator of this, using input/output Prolongation matrices Pi corresponds to "P"...
int Size() const
Return the logical size of the array.
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
Class for parallel grid function.
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering e_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
Wrapper for hypre's ParCSR matrix class.
Bitwise-OR of all CEED backends.
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Operator that converts L2 FiniteElementSpace L-vectors to E-vectors.
int width
Dimension of the input / number of columns in the matrix.
MemoryType GetMemoryType() const
Return a MemoryType that is currently valid. If both the host and the device pointers are currently v...
void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const
y += c*A^T*x
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.