38 Table face_dof, dof_face;
57 int *I = dof_dof.
GetI();
58 int *J = dof_dof.
GetJ();
113 MFEM_ABORT(
"the assembly level has already been set!");
134 MFEM_ABORT(
"BilinearForm: unknown assembly level");
143 MFEM_WARNING(
"Static condensation not supported for this assembly level");
149 bool symmetric =
false;
150 bool block_diagonal =
false;
167 MFEM_WARNING(
"Hybridization not supported for this assembly level");
198 "invalid matrix A dimensions: "
200 MFEM_ASSERT(A.
Finalized(),
"matrix A must be Finalized");
290 domain_integs[0]->AssembleElementMatrix(fe, *eltrans, elmat);
331 tr = mesh -> GetFaceElementTransformations (i);
335 if (tr->Elem2No >= 0)
359 if (tr->Elem2No >= 0)
373 tr = mesh -> GetBdrFaceTransformations (i);
377 fe1 =
fes -> GetFE (tr -> Elem1No);
478#ifdef MFEM_USE_LEGACY_OPENMP
479 int free_element_matrices = 0;
483 free_element_matrices = 1;
495 "invalid element marker for domain integrator #"
496 << k <<
", counting from zero");
501 MFEM_VERIFY(
fes->
GetNURBSext(),
"Patchwise integration requires a "
502 <<
"NURBS FE space");
507 for (
int i = 0; i <
fes -> GetNE(); i++)
515 elmat_p = &(*element_matrices)(i);
532 if (elmat.
Size() == 0)
542 if (elmat.
Size() == 0)
575 bool vdofsSet =
false;
591 for (
int r=0; r<spmat->
Height(); ++r)
593 spmat->
GetRow(r, cols, srow);
594 for (
int i=0; i<cols.
Size(); ++i)
596 cols[i] =
vdofs[cols[i]];
622 MFEM_ASSERT(bdr_marker.
Size() == bdr_attr_marker.
Size(),
623 "invalid boundary marker for boundary integrator #"
624 << k <<
", counting from zero");
625 for (
int i = 0; i < bdr_attr_marker.
Size(); i++)
627 bdr_attr_marker[i] |= bdr_marker[i];
631 for (
int i = 0; i <
fes -> GetNBE(); i++)
634 if (bdr_attr_marker[bdr_attr-1] == 0) {
continue; }
637 doftrans =
fes -> GetBdrElementVDofs (i,
vdofs);
638 eltrans =
fes -> GetBdrElementTransformation (i);
683 for (
int i = 0; i < nfaces; i++)
685 tr = mesh -> GetInteriorFaceTransformations (i);
688 fes -> GetElementVDofs (tr -> Elem1No,
vdofs);
689 fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
694 AssembleFaceMatrix(*
fes->
GetFE(tr->Elem1No),
720 MFEM_ASSERT(bdr_marker.
Size() == bdr_attr_marker.
Size(),
721 "invalid boundary marker for boundary face integrator #"
722 << k <<
", counting from zero");
723 for (
int i = 0; i < bdr_attr_marker.
Size(); i++)
725 bdr_attr_marker[i] |= bdr_marker[i];
729 for (
int i = 0; i <
fes -> GetNBE(); i++)
732 if (bdr_attr_marker[bdr_attr-1] == 0) {
continue; }
734 tr = mesh -> GetBdrFaceTransformations (i);
737 fes -> GetElementVDofs (tr -> Elem1No,
vdofs);
738 fe1 =
fes -> GetFE (tr -> Elem1No);
757#ifdef MFEM_USE_LEGACY_OPENMP
758 if (free_element_matrices)
772 MFEM_ASSERT(
mat,
"the BilinearForm is not assembled");
803 "Vector for holding diagonal has wrong size!");
807 MFEM_ASSERT(
mat,
"the BilinearForm is not assembled!");
809 "BilinearForm::ConformingAssemble() is not called!");
816 ext->AssembleDiagonal(diag);
825 ext->AssembleDiagonal(local_diag);
831 Vector &B,
int copy_interior)
846 ext->FormLinearSystem(ess_tdof_list, x,
b, A, X, B, copy_interior);
917 const int remove_zeros = 0;
923 ext->FormSystemMatrix(ess_tdof_list, A);
948 const int remove_zeros = 0;
967 ext->RecoverFEMSolution(X,
b, x);
1046#ifdef MFEM_USE_LEGACY_OPENMP
1047 #pragma omp parallel for private(tmp,eltrans)
1049 for (
int i = 0; i < num_elements; i++)
1052 num_dofs_per_el, num_dofs_per_el);
1056 mfem_error(
"BilinearForm::ComputeElementMatrices:"
1057 " all elements must have same number of dofs");
1061 domain_integs[0]->AssembleElementMatrix(fe, eltrans, elmat);
1135 for (
int i = 0; i < vdofs_.
Size(); i++)
1137 int vdof = vdofs_[i];
1140 mat -> EliminateRowCol (vdof, sol(vdof), rhs, dpolicy);
1144 mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, dpolicy);
1158 for (
int i = 0; i < vdofs_.
Size(); i++)
1160 int vdof = vdofs_[i];
1163 mat -> EliminateRowCol (vdof, *
mat_e, dpolicy);
1167 mat -> EliminateRowCol (-1-vdof, *
mat_e, dpolicy);
1176 MFEM_ASSERT(ess_dofs.
Size() ==
height,
"incorrect dof Array size");
1177 MFEM_ASSERT(sol.
Size() ==
height,
"incorrect sol Vector size");
1178 MFEM_ASSERT(rhs.
Size() ==
height,
"incorrect rhs Vector size");
1180 for (
int i = 0; i < ess_dofs.
Size(); i++)
1181 if (ess_dofs[i] < 0)
1183 mat -> EliminateRowCol (i, sol(i), rhs, dpolicy);
1191 "incorrect dof Array size: " << ess_dofs.
Size() <<
' ' <<
height);
1193 for (
int i = 0; i < ess_dofs.
Size(); i++)
1194 if (ess_dofs[i] < 0)
1196 mat -> EliminateRowCol (i, dpolicy);
1204 "incorrect dof Array size: " << ess_dofs.
Size() <<
' ' <<
height);
1206 for (
int i = 0; i < ess_dofs.
Size(); i++)
1207 if (ess_dofs[i] < 0)
1209 mat -> EliminateRowColDiag (i, value);
1236 ext->MultTranspose(x, y);
1249 if (nfes && nfes !=
fes)
1282 if (
ext) {
ext->Update(); }
1310 :
Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1324 :
Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1352 MFEM_ABORT(
"the assembly level has already been set!");
1364 MFEM_ABORT(
"Element assembly not supported yet... stay tuned!");
1371 MFEM_ABORT(
"Matrix-free action not supported yet... stay tuned!");
1375 MFEM_ABORT(
"Unknown assembly level");
1381 return (*
mat)(i, j);
1386 return (*
mat)(i, j);
1400 ext->AddMult(x, y,
a);
1419 ext->AddMultTranspose(x, y,
a);
1431 MFEM_WARNING(
"MixedBilinearForm::Inverse not possible with this "
1453 "MixedBilinearForm::GetBlocks: both trial and test spaces "
1454 "must use Ordering::byNODES!");
1552 "invalid element marker for domain integrator #"
1553 << k <<
", counting from zero");
1557 for (
int i = 0; i <
test_fes -> GetNE(); i++)
1562 eltrans =
test_fes -> GetElementTransformation (i);
1577 if (ran_dof_trans || dom_dof_trans)
1590 bdr_attr_marker = 0;
1595 bdr_attr_marker = 1;
1599 MFEM_ASSERT(bdr_marker.
Size() == bdr_attr_marker.
Size(),
1600 "invalid boundary marker for boundary integrator #"
1601 << k <<
", counting from zero");
1602 for (
int i = 0; i < bdr_attr_marker.
Size(); i++)
1604 bdr_attr_marker[i] |= bdr_marker[i];
1608 for (
int i = 0; i <
test_fes -> GetNBE(); i++)
1611 if (bdr_attr_marker[bdr_attr-1] == 0) {
continue; }
1615 eltrans =
test_fes -> GetBdrElementTransformation (i);
1629 if (ran_dof_trans || dom_dof_trans)
1641 const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
1644 for (
int i = 0; i < nfaces; i++)
1667 trial_fe2 = trial_fe1;
1668 test_fe2 = test_fe1;
1685 const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
1690 bdr_attr_marker = 0;
1695 bdr_attr_marker = 1;
1699 MFEM_ASSERT(bdr_marker.
Size() == bdr_attr_marker.
Size(),
1700 "invalid boundary marker for boundary face integrator #"
1701 << k <<
", counting from zero");
1702 for (
int i = 0; i < bdr_attr_marker.
Size(); i++)
1704 bdr_attr_marker[i] |= bdr_marker[i];
1708 for (
int i = 0; i <
trial_fes -> GetNBE(); i++)
1711 if (bdr_attr_marker[bdr_attr-1] == 0) {
continue; }
1713 ftr = mesh -> GetBdrFaceTransformations (i);
1723 trial_fe2 = trial_fe1;
1724 test_fe2 = test_fe1;
1746 for (
int i = 0; i < nfaces; i++)
1764 test_fe2 = test_fe1;
1784 bdr_attr_marker = 0;
1789 bdr_attr_marker = 1;
1793 MFEM_ASSERT(bdr_marker.
Size() == bdr_attr_marker.
Size(),
1794 "invalid boundary marker for boundary trace face"
1795 "integrator #" << k <<
", counting from zero");
1796 for (
int i = 0; i < bdr_attr_marker.
Size(); i++)
1798 bdr_attr_marker[i] |= bdr_marker[i];
1802 for (
int i = 0; i <
trial_fes -> GetNBE(); i++)
1805 if (bdr_attr_marker[bdr_attr-1] == 0) {
continue; }
1818 test_fe2 = test_fe1;
1842 "Vector for holding diagonal has wrong size!");
1844 "Vector for holding diagonal has wrong size!");
1850 P_trial->
Mult(D, local_D);
1855 ext->AssembleDiagonal_ADAt(local_D, local_diag);
1860 ext->AssembleDiagonal_ADAt(local_D, diag);
1868 ext->AssembleDiagonal_ADAt(D, local_diag);
1873 ext->AssembleDiagonal_ADAt(D, diag);
1879 MFEM_ABORT(
"Not implemented. Maybe assemble your bilinear form into a "
1880 "matrix and use SparseMatrix functions?");
1888 MFEM_WARNING(
"Conforming assemble not supported for this assembly level!");
1925 domain_integs[0]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
1929 domain_integs[k]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
1939 elmat.
SetSize(te_dofs, tr_dofs);
1952 boundary_integs[0]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
1956 boundary_integs[k]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
1966 elmat.
SetSize(te_dofs, tr_dofs);
1976 MFEM_ASSERT(ftr,
"No associated face transformations.");
1978 const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
1992 trial_fe2 = trial_fe1;
1993 test_fe2 = test_fe1;
2019 elmat.
SetSize(te_dofs, tr_dofs);
2029 MFEM_ASSERT(ftr,
"No associated boundary face.");
2031 const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
2038 trial_fe2 = trial_fe1;
2039 test_fe2 = test_fe1;
2059 elmat.
SetSize(te_dofs, tr_dofs);
2069 MFEM_ASSERT(ftr,
"No associated face transformation.");
2084 test_fe2 = test_fe1;
2107 elmat.
SetSize(te_dofs, tr_face_dofs);
2118 MFEM_ASSERT(ftr,
"No associated boundary face.");
2127 test_fe2 = test_fe1;
2147 elmat.
SetSize(te_dofs, tr_face_dofs);
2211 trial_vdofs_marker);
2224 trial_vdofs_marker);
2247 for (i = 0; i <
test_fes -> GetNBE(); i++)
2248 if (bdr_attr_is_ess[
test_fes -> GetBdrAttribute (i)-1])
2250 test_fes -> GetBdrElementVDofs (i, te_vdofs);
2251 for (j = 0; j < te_vdofs.
Size(); j++)
2253 if ( (k = te_vdofs[j]) < 0 )
2257 mat -> EliminateRow (k);
2264 for (
int i=0; i<test_vdofs_.
Size(); ++i)
2278 ext->FormRectangularSystemOperator(trial_tdof_list, test_tdof_list, A);
2287 if (test_P && trial_P)
2321 ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
2350 if (
ext) {
ext->Update(); }
2378 MFEM_ABORT(
"the assembly level has already been set!");
2388 MFEM_ABORT(
"Element assembly not supported yet... stay tuned!");
2394 MFEM_ABORT(
"Matrix-free action not supported yet... stay tuned!");
2397 MFEM_ABORT(
"Unknown assembly level");
2429 "invalid element marker for domain integrator #"
2430 << k <<
", counting from zero");
2454 if (ran_dof_trans || dom_dof_trans)
2465 for (
int i = 0; i < nfaces; i++)
Dynamic 2D array using row-major layout.
void SetSize(int m, int n)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
int Size() const
Return the logical size of the array.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Square Operator for imposing essential boundary conditions using only the action, Mult(),...
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate "essential boundary condition" values specified in x from the given right-hand side b.
Data type dense matrix using column-major storage.
real_t * GetData() const
Returns the matrix data array.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
int Size() const
For backward compatibility define Size to be synonym of Width()
Rank 3 tensor (array of matrices)
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::FULL.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element,...
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Convert an array of indices (list) to a Boolean marker array where all indices in the list are marked...
const NURBSExtension * GetNURBSext() const
virtual const Operator * GetProlongationMatrix() const
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Ordering::Type GetOrdering() const
Return the ordering method.
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
int GetNE() const
Returns number of elements in the mesh.
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
const SparseMatrix * GetConformingProlongation() const
Mesh * GetMesh() const
Returns the mesh.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
int GetVDim() const
Returns the vector dimension of the finite element space.
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Mark degrees of freedom associated with boundary elements with the specified boundary attributes (mar...
void GetPatchVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom in vdofs for NURBS patch i.
Abstract class for all finite elements.
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Auxiliary class Hybridization, used to implement BilinearForm hybridization.
Abstract data type for matrix inverse.
Abstract data type matrix.
Class used by MFEM to store pointers to host and/or device memory.
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
NURBSExtension * NURBSext
Optional NURBS mesh extension.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
int GetAttribute(int i) const
Return the attribute of element i.
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
int GetBdrElementFaceIndex(int be_idx) const
Return the local face (codimension-1) index for the given boundary element index.
ElementTransformation * GetFaceTransformation(int FaceNo)
Returns a pointer to the transformation defining the given face element.
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
See GetFaceElementTransformations().
Table * GetFaceToElementTable() const
Array< int > attributes
A list of all unique element attributes used by the Mesh.
int GetNP() const
Return the number of patches.
Pointer to an Operator of a specified type.
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
int width
Dimension of the input / number of columns in the matrix.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int height
Dimension of the output / number of rows in the matrix.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
@ DIAG_ONE
Set the diagonal value to one.
@ DIAG_KEEP
Keep the diagonal value.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
void InitTVectors(const Operator *Po, const Operator *Ri, const Operator *Pi, Vector &x, Vector &b, Vector &X, Vector &B) const
Initializes memory for true vectors of linear system.
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 ...
Partial assembly extension for DiscreteLinearOperator.
int GetRow(const int row, Array< int > &cols, Vector &srow) const override
Extract all column indices and values from a given row.
void GetDiag(Vector &d) const
Returns the Diagonal of A.
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
bool ColumnsAreSorted() const
Returns whether or not the columns are sorted.
void MultTranspose(const Vector &x, Vector &y) const override
Multiply a vector with the transposed matrix. y = At * x.
bool Finalized() const
Returns whether or not CSR format has been finalized.
void AbsMultTranspose(const Vector &x, Vector &y) const
y = |At| * x, using entry-wise absolute values of the transpose of matrix A
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const override
y += At * x (default) or y += a * At * x
void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const override
y += A * x (default) or y += a * A * x
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication.
void EliminateRow(int row, const real_t sol, Vector &rhs)
Eliminates a column from the transpose matrix.
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
void Finalize(int skip_zeros=1) override
Finalize the matrix initialization, switching the storage format from LIL to CSR.
void EliminateCols(const Array< int > &cols, const Vector *x=NULL, Vector *b=NULL)
Eliminate all columns i for which cols[i] != 0.
int * GetJ()
Return the array J.
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
int * GetI()
Return the array I.
void LoseData()
Call this if data has been stolen.
void SortRows()
Sort the column (TYPE II) indices in each row.
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
void SyncMemory(const Vector &v) const
Update the memory location of the vector to match v.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
void SetSubVectorComplement(const Array< int > &dofs, const real_t val)
Set all vector entries NOT in the dofs Array to the given val.
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
void mfem_error(const char *msg)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
void TransformDual(const DofTransformation *ran_dof_trans, const DofTransformation *dom_dof_trans, DenseMatrix &elmat)
void TransformPrimal(const DofTransformation *ran_dof_trans, const DofTransformation *dom_dof_trans, DenseMatrix &elmat)
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
bool IsIdentityProlongation(const Operator *P)
@ NATIVE
Native ordering as defined by the FiniteElement.
std::function< real_t(const Vector &)> f(real_t mass_coeff)
SparseMatrix * TransposeMult(const SparseMatrix &A, const SparseMatrix &B)
C = A^T B.
real_t p(const Vector &x, real_t t)