15 #include "../general/device.hpp"
37 Table face_dof, dof_face;
56 int *I = dof_dof.
GetI();
57 int *J = dof_dof.
GetJ();
119 MFEM_ABORT(
"the assembly level has already been set!");
149 MFEM_WARNING(
"Static condensation not supported for this assembly level");
155 bool symmetric =
false;
156 bool block_diagonal =
false;
175 MFEM_WARNING(
"Hybridization not supported for this assembly level");
202 "invalid matrix A dimensions: "
204 MFEM_ASSERT(A.
Finalized(),
"matrix A must be Finalized");
293 domain_integs[0]->AssembleElementMatrix(fe, *eltrans, elmat);
402 #ifdef MFEM_USE_LEGACY_OPENMP
403 int free_element_matrices = 0;
407 free_element_matrices = 1;
419 "invalid element marker for domain integrator #"
420 << k <<
", counting from zero");
424 for (
int i = 0; i <
fes -> GetNE(); i++)
430 elmat_p = &(*element_matrices)(i);
443 if (elmat.
Size() == 0)
453 if (elmat.
Size() == 0)
491 MFEM_ASSERT(bdr_marker.
Size() == bdr_attr_marker.
Size(),
492 "invalid boundary marker for boundary integrator #"
493 << k <<
", counting from zero");
494 for (
int i = 0; i < bdr_attr_marker.
Size(); i++)
496 bdr_attr_marker[i] |= bdr_marker[i];
500 for (
int i = 0; i <
fes -> GetNBE(); i++)
503 if (bdr_attr_marker[bdr_attr-1] == 0) {
continue; }
506 fes -> GetBdrElementVDofs (i,
vdofs);
507 eltrans =
fes -> GetBdrElementTransformation (i);
547 for (
int i = 0; i < nfaces; i++)
549 tr = mesh -> GetInteriorFaceTransformations (i);
552 fes -> GetElementVDofs (tr -> Elem1No,
vdofs);
553 fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
584 MFEM_ASSERT(bdr_marker.
Size() == bdr_attr_marker.
Size(),
585 "invalid boundary marker for boundary face integrator #"
586 << k <<
", counting from zero");
587 for (
int i = 0; i < bdr_attr_marker.
Size(); i++)
589 bdr_attr_marker[i] |= bdr_marker[i];
593 for (
int i = 0; i <
fes -> GetNBE(); i++)
596 if (bdr_attr_marker[bdr_attr-1] == 0) {
continue; }
598 tr = mesh -> GetBdrFaceTransformations (i);
601 fes -> GetElementVDofs (tr -> Elem1No,
vdofs);
602 fe1 =
fes -> GetFE (tr -> Elem1No);
621 #ifdef MFEM_USE_LEGACY_OPENMP
622 if (free_element_matrices)
636 MFEM_ASSERT(
mat,
"the BilinearForm is not assembled");
667 "Vector for holding diagonal has wrong size!");
671 MFEM_ASSERT(
mat,
"the BilinearForm is not assembled!");
672 MFEM_ASSERT(cP ==
nullptr ||
mat->
Height() == cP->Width(),
673 "BilinearForm::ConformingAssemble() is not called!");
688 Vector local_diag(cP->Height());
690 cP->AbsMultTranspose(local_diag, diag);
695 Vector &B,
int copy_interior)
791 const int remove_zeros = 0;
879 #ifdef MFEM_USE_LEGACY_OPENMP
880 #pragma omp parallel for private(tmp,eltrans)
882 for (
int i = 0; i < num_elements; i++)
885 num_dofs_per_el, num_dofs_per_el);
888 if (num_dofs_per_el != fe.GetDof()*
fes->
GetVDim())
889 mfem_error(
"BilinearForm::ComputeElementMatrices:"
890 " all elements must have same number of dofs");
961 for (
int i = 0; i < vdofs.
Size(); i++)
966 mat -> EliminateRowCol (vdof, sol(vdof), rhs, dpolicy);
970 mat -> EliminateRowCol (-1-vdof, sol(-1-vdof), rhs, dpolicy);
983 for (
int i = 0; i < vdofs.
Size(); i++)
988 mat -> EliminateRowCol (vdof, *
mat_e, dpolicy);
992 mat -> EliminateRowCol (-1-vdof, *
mat_e, dpolicy);
1001 MFEM_ASSERT(ess_dofs.
Size() ==
height,
"incorrect dof Array size");
1002 MFEM_ASSERT(sol.
Size() ==
height,
"incorrect sol Vector size");
1003 MFEM_ASSERT(rhs.
Size() ==
height,
"incorrect rhs Vector size");
1005 for (
int i = 0; i < ess_dofs.
Size(); i++)
1006 if (ess_dofs[i] < 0)
1008 mat -> EliminateRowCol (i, sol(i), rhs, dpolicy);
1015 MFEM_ASSERT(ess_dofs.
Size() ==
height,
"incorrect dof Array size");
1017 for (
int i = 0; i < ess_dofs.
Size(); i++)
1018 if (ess_dofs[i] < 0)
1020 mat -> EliminateRowCol (i, dpolicy);
1027 MFEM_ASSERT(ess_dofs.
Size() ==
height,
"incorrect dof Array size");
1029 for (
int i = 0; i < ess_dofs.
Size(); i++)
1030 if (ess_dofs[i] < 0)
1032 mat -> EliminateRowColDiag (i, value);
1059 if (nfes && nfes !=
fes)
1127 :
Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1141 :
Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
1167 MFEM_ABORT(
"the assembly level has already been set!");
1179 mfem_error(
"Element assembly not supported yet... stay tuned!");
1186 mfem_error(
"Matrix-free action not supported yet... stay tuned!");
1196 return (*
mat)(i, j);
1201 return (*
mat)(i, j);
1211 const double a)
const
1230 const double a)
const
1246 MFEM_WARNING(
"MixedBilinearForm::Inverse not possible with this "
1268 "MixedBilinearForm::GetBlocks: both trial and test spaces "
1269 "must use Ordering::byNODES!");
1334 for (
int i = 0; i <
test_fes -> GetNE(); i++)
1336 trial_fes -> GetElementVDofs (i, tr_vdofs);
1337 test_fes -> GetElementVDofs (i, te_vdofs);
1338 eltrans =
test_fes -> GetElementTransformation (i);
1344 mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1354 bdr_attr_marker = 0;
1359 bdr_attr_marker = 1;
1363 MFEM_ASSERT(bdr_marker.
Size() == bdr_attr_marker.
Size(),
1364 "invalid boundary marker for boundary integrator #"
1365 << k <<
", counting from zero");
1366 for (
int i = 0; i < bdr_attr_marker.
Size(); i++)
1368 bdr_attr_marker[i] |= bdr_marker[i];
1372 for (
int i = 0; i <
test_fes -> GetNBE(); i++)
1375 if (bdr_attr_marker[bdr_attr-1] == 0) {
continue; }
1377 trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1378 test_fes -> GetBdrElementVDofs (i, te_vdofs);
1379 eltrans =
test_fes -> GetBdrElementTransformation (i);
1388 mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
1400 for (
int i = 0; i < nfaces; i++)
1410 te_vdofs.
Append(te_vdofs2);
1418 test_fe2 = test_fe1;
1423 *test_fe2, *ftr, elemmat);
1438 bdr_attr_marker = 0;
1443 bdr_attr_marker = 1;
1447 MFEM_ASSERT(bdr_marker.
Size() == bdr_attr_marker.
Size(),
1448 "invalid boundary marker for boundary trace face"
1449 "integrator #" << k <<
", counting from zero");
1450 for (
int i = 0; i < bdr_attr_marker.
Size(); i++)
1452 bdr_attr_marker[i] |= bdr_marker[i];
1456 for (
int i = 0; i <
trial_fes -> GetNBE(); i++)
1459 if (bdr_attr_marker[bdr_attr-1] == 0) {
continue; }
1471 test_fe2 = test_fe1;
1495 "Vector for holding diagonal has wrong size!");
1497 "Vector for holding diagonal has wrong size!");
1502 Vector local_D(P_trial->Height());
1503 P_trial->Mult(D, local_D);
1532 MFEM_ABORT(
"Not implemented. Maybe assemble your bilinear form into a "
1533 "matrix and use SparseMatrix functions?");
1541 MFEM_WARNING(
"Conforming assemble not supported for this assembly level!");
1577 domain_integs[0]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
1581 domain_integs[k]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
1602 boundary_integs[0]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
1606 boundary_integs[k]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
1665 for (i = 0; i <
trial_fes -> GetNBE(); i++)
1666 if (bdr_attr_is_ess[
trial_fes -> GetBdrAttribute (i)-1])
1668 trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
1669 for (j = 0; j < tr_vdofs.
Size(); j++)
1671 if ( (k = tr_vdofs[j]) < 0 )
1678 mat -> EliminateCols (cols_marker, &sol, &rhs);
1684 mat -> EliminateCols (marked_vdofs, &sol, &rhs);
1692 for (i = 0; i <
test_fes -> GetNBE(); i++)
1693 if (bdr_attr_is_ess[
test_fes -> GetBdrAttribute (i)-1])
1695 test_fes -> GetBdrElementVDofs (i, te_vdofs);
1696 for (j = 0; j < te_vdofs.
Size(); j++)
1698 if ( (k = te_vdofs[j]) < 0 )
1702 mat -> EliminateRow (k);
1731 Array<int> ess_trial_tdof_marker, ess_test_tdof_marker;
1733 ess_trial_tdof_marker);
1735 ess_test_tdof_marker);
1740 for (
int i=0; i<test_tdof_list.
Size(); ++i)
1811 MFEM_ABORT(
"the assembly level has already been set!");
1821 mfem_error(
"Element assembly not supported yet... stay tuned!");
1827 mfem_error(
"Matrix-free action not supported yet... stay tuned!");
1862 domain_integs[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
1866 domain_integs[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
1877 for (
int i = 0; i < nfaces; i++)
Abstract class for all finite elements.
int Size() const
For backward compatibility define Size to be synonym of Width()
Ordering::Type GetOrdering() const
Return the ordering method.
int Size() const
Return the logical size of the array.
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
void SetConstraintIntegrator(BilinearFormIntegrator *c_integ)
bool ColumnsAreSorted() const
Returns whether or not the columns are sorted.
bool ReducesTrueVSize() const
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
void SortRows()
Sort the column (TYPE II) indices in each row.
void SetSize(int s)
Resize the vector to size s.
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Pointer to an Operator of a specified type.
int * GetJ()
Return the array J.
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...
Data type dense matrix using column-major storage.
void AddMult(const Vector &x, Vector &y, const double a=1.0) const
y += A * x (default) or y += a * A * x
int Size() const
Returns the size of the vector.
void ReduceSystem(Vector &x, Vector &b, Vector &X, Vector &B, int copy_interior=0) const
Set the reduced solution X and r.h.s B vectors from the full linear system solution x and r...
int * GetI()
Return the array I.
void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
y += At * x (default) or y += a * At * x
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, Matrix::DiagonalPolicy dpolicy)
Eliminate the given reduced true dofs from the Schur complement matrix S.
Abstract data type for matrix inverse.
void LoseData()
Call this if data has been stolen.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
void ComputeSolution(const Vector &b, const Vector &sc_sol, Vector &sol) const
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::FULL.
void Reset()
Destroy the current hybridization matrix while preserving the computed constraint matrix and the set ...
void GetBlocks(Array2D< SparseMatrix * > &blocks) 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 ...
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
int GetNE() const
Returns number of elements in the mesh.
void SetSize(int m, int n)
void Finalize()
Finalize the construction of the hybridized matrix.
void AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
Assemble the boundary element matrix A into the hybridized system matrix.
double f(const Vector &xvec)
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
SparseMatrix & GetMatrix()
Return the serial hybridized matrix.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Mesh * GetMesh() const
Returns the mesh.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
void GetFaceTransformation(int i, IsoparametricTransformation *FTr)
Returns the transformation defining the given face element in a user-defined variable.
void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const
Abstract data type matrix.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
int GetVDim() const
Returns vector dimension.
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
void SyncMemory(const Vector &v) const
Update the memory location of the vector to match v.
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Dynamic 2D array using row-major layout.
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the dofs Array to the given val.
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
bool Finalized() const
Returns whether or not CSR format has been finalized.
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
bool IsIdentityProlongation(const Operator *P)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void AssembleMatrix(int el, const DenseMatrix &A)
Assemble the element matrix A into the hybridized system matrix.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
void Finalize()
Finalize the construction of the Schur complement matrix.
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.
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
void AssembleBdrMatrix(int el, const DenseMatrix &elmat)
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
int height
Dimension of the output / number of rows in the matrix.
void EliminateRow(int row, const double sol, Vector &rhs)
Eliminates a column from the transpose matrix.
Table * GetFaceToElementTable() const
void GetDiag(Vector &d) const
Returns the Diagonal of A.
void AssembleMatrix(int el, const DenseMatrix &elmat)
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix...
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
void EliminateCols(const Array< int > &cols, const Vector *x=NULL, Vector *b=NULL)
Eliminate all columns i for which cols[i] != 0.
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th face element (2D and 3D).
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...
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
SparseMatrix & GetMatrix()
Return the serial Schur complement matrix.
void ReduceRHS(const Vector &b, Vector &b_r) const
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element...
void Init(const Array< int > &ess_tdof_list)
Prepare the Hybridization object for assembly.
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th boundary element.
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
void Init(bool symmetric, bool block_diagonal)
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
Rank 3 tensor (array of matrices)
int GetAttribute(int i) const
Return the attribute of element i.
Array< int > attributes
A list of all unique element attributes used by the Mesh.
int width
Dimension of the input / number of columns in the matrix.
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Partial assembly extension for DiscreteLinearOperator.