124 #include "../common/mfem-common.hpp" 129 using namespace mfem;
167 static const char *enum_str[] =
173 "pml_plane_wave_scatter",
179 int main(
int argc,
char *argv[])
182 int myid = Mpi::WorldRank();
185 const char *mesh_file =
"../../data/inline-quad.mesh";
188 bool visualization =
true;
191 bool static_cond =
false;
196 bool with_pml =
false;
197 bool paraview =
false;
200 args.
AddOption(&mesh_file,
"-m",
"--mesh",
201 "Mesh file to use.");
203 "Finite element order (polynomial degree)");
204 args.
AddOption(&rnum,
"-rnum",
"--number-of-wavelengths",
205 "Number of wavelengths");
206 args.
AddOption(&iprob,
"-prob",
"--problem",
"Problem case" 207 " 0: plane wave, 1: Gaussian beam, 2: Generic PML," 208 " 3: Scattering of a Gaussian beam" 209 " 4: Scattering of a plane wave, 5: Point source");
210 args.
AddOption(&delta_order,
"-do",
"--delta-order",
211 "Order enrichment for DPG test space.");
212 args.
AddOption(&theta,
"-theta",
"--theta",
213 "Theta parameter for AMR");
214 args.
AddOption(&sr,
"-sref",
"--serial-ref",
215 "Number of parallel refinements.");
216 args.
AddOption(&pr,
"-pref",
"--parallel-ref",
217 "Number of parallel refinements.");
218 args.
AddOption(&static_cond,
"-sc",
"--static-condensation",
"-no-sc",
219 "--no-static-condensation",
"Enable static condensation.");
220 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
221 "--no-visualization",
222 "Enable or disable GLVis visualization.");
223 args.
AddOption(¶view,
"-paraview",
"--paraview",
"-no-paraview",
225 "Enable or disable ParaView visualization.");
236 if (iprob > 5) { iprob = 0; }
238 omega = 2.*M_PI*rnum;
243 if (
prob > 2) { mesh_file =
"meshes/scatter.mesh"; }
255 Mesh mesh(mesh_file, 1, 1);
257 for (
int i = 0; i<sr; i++)
262 MFEM_VERIFY(
dim > 1,
"Dimension = 1 is not supported in this example");
273 ParMesh pmesh(MPI_COMM_WORLD, mesh);
312 int test_order = order+delta_order;
320 trial_fes.
Append(hatp_fes);
321 trial_fes.
Append(hatu_fes);
345 negomeg_cf = &negomeg;
381 omeg2_abs_Jt_J_detJinv_2,attrPML);
390 TrialSpace::p_space,TestSpace::q_space);
393 nullptr,TrialSpace::u_space,TestSpace::q_space);
396 TrialSpace::p_space,TestSpace::v_space);
398 a->AddTrialIntegrator(
nullptr,
400 TrialSpace::u_space,TestSpace::v_space);
403 TrialSpace::hatp_space,TestSpace::v_space);
406 TrialSpace::hatu_space,TestSpace::q_space);
411 TestSpace::q_space, TestSpace::q_space);
414 TestSpace::q_space, TestSpace::q_space);
417 TestSpace::v_space, TestSpace::v_space);
420 TestSpace::v_space, TestSpace::v_space);
423 TestSpace::q_space, TestSpace::v_space);
425 a->AddTestIntegrator(
nullptr,
427 TestSpace::v_space, TestSpace::q_space);
430 TestSpace::v_space, TestSpace::v_space);
433 TestSpace::v_space, TestSpace::q_space);
436 TestSpace::q_space, TestSpace::v_space);
439 TestSpace::q_space, TestSpace::q_space);
452 TrialSpace::p_space,TestSpace::q_space);
460 TrialSpace::u_space,TestSpace::v_space);
465 omeg_Jt_J_detJinv_i_restr);
468 negomeg_Jt_J_detJinv_r_restr);
470 a->AddTestIntegrator(integ0_r, integ0_i,
471 TestSpace::q_space,TestSpace::v_space);
477 negomeg_Jt_J_detJinv_i_restr),
479 TestSpace::v_space,TestSpace::q_space);
483 omeg2_abs_Jt_J_detJinv_2_restr);
485 a->AddTestIntegrator(integ1,
nullptr,TestSpace::v_space,TestSpace::v_space);
491 TestSpace::v_space,TestSpace::q_space);
496 negomeg_detJ_i_restr),
498 TestSpace::q_space,TestSpace::v_space);
503 a->AddTestIntegrator(integ,
nullptr,
504 TestSpace::q_space,TestSpace::q_space);
532 std::cout <<
"\n Ref |" 537 std::cout <<
" L2 Error |" 540 std::cout <<
" Residual |" 542 <<
" PCG it |" << endl;
543 std::cout << std::string((
exact_known) ? 82 : 60,
'-')
570 if (static_cond) {
a->EnableStaticCondensation(); }
571 for (
int it = 0; it<=pr; it++)
590 for (
int j = 0; j < ess_tdof_list.
Size(); j++)
611 hatp_gf_i.ProjectBdrCoefficient(hatpex_i, ess_bdr);
616 a->FormLinearSystem(ess_tdof_list,x,Ah, X,B);
626 int skip = (static_cond) ? 0 : 2;
627 int k = (static_cond) ? 2 : 0;
628 for (
int i=0; i<num_blocks; i++)
630 tdof_offsets[i+1] = trial_fes[i+k]->GetTrueVSize();
631 tdof_offsets[num_blocks+i+1] = trial_fes[i+k]->GetTrueVSize();
636 for (
int i = 0; i<num_blocks; i++)
638 for (
int j = 0; j<num_blocks; j++)
640 blockA.
SetBlock(i,j,&BlockA_r->GetBlock(i,j));
642 blockA.
SetBlock(i+num_blocks,j+num_blocks,&BlockA_r->GetBlock(i,j));
654 BlockA_r->GetBlock(0,0));
658 BlockA_r->GetBlock(1,1));
668 BlockA_r->GetBlock(skip,skip));
677 dynamic_cast<HypreAMS*
>(solver_hatu)->SetPrintLevel(0);
684 dynamic_cast<HypreADS*
>(solver_hatu)->SetPrintLevel(0);
700 for (
int i = 0; i<num_blocks; i++)
707 a->RecoverFEMSolution(X,x);
709 Vector & residuals =
a->ComputeResidual(x);
711 double residual = residuals.
Norml2();
712 double maxresidual = residuals.
Max();
713 double globalresidual = residual * residual;
714 MPI_Allreduce(MPI_IN_PLACE,&maxresidual,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
715 MPI_Allreduce(MPI_IN_PLACE,&globalresidual,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
717 globalresidual = sqrt(globalresidual);
722 u_r.
MakeRef(u_fes,x, offsets[1]);
726 for (
int i = 0; i<trial_fes.
Size(); i++)
728 dofs += trial_fes[i]->GlobalTrueVSize();
731 double L2Error = 0.0;
732 double rate_err = 0.0;
747 L2Error = sqrt(p_err_r*p_err_r + p_err_i*p_err_i
748 +u_err_r*u_err_r + u_err_i*u_err_i);
750 rate_err = (it) ?
dim*log(err0/L2Error)/log((
double)dof0/dofs) : 0.0;
754 double rate_res = (it) ?
dim*log(res0/globalresidual)/log((
755 double)dof0/dofs) : 0.0;
757 res0 = globalresidual;
762 std::ios oldState(
nullptr);
763 oldState.copyfmt(std::cout);
764 std::cout << std::right << std::setw(5) << it <<
" | " 765 << std::setw(10) << dof0 <<
" | " 766 << std::setprecision(1) << std::fixed
767 << std::setw(4) << 2*rnum <<
" π | ";
770 std::cout << std::setprecision(3) << std::setw(10)
771 << std::scientific << err0 <<
" | " 772 << std::setprecision(2)
773 << std::setw(6) << std::fixed << rate_err <<
" | " ;
775 std::cout << std::setprecision(3)
776 << std::setw(10) << std::scientific << res0 <<
" | " 777 << std::setprecision(2)
778 << std::setw(6) << std::fixed << rate_res <<
" | " 779 << std::setw(6) << std::fixed << num_iter <<
" | " 781 std::cout.copyfmt(oldState);
786 const char * keys = (it == 0 &&
dim == 2) ?
"jRcml\n" :
nullptr;
790 "Numerical presure (real part)", 0, 0, 500, 500, keys);
792 "Numerical presure (imaginary part)", 501, 0, 500, 500, keys);
798 paraview_dc->
SetTime((
double)it);
810 for (
int iel = 0; iel<pmesh.
GetNE(); iel++)
812 if (residuals[iel] > theta * maxresidual)
814 elements_to_refine.
Append(iel);
824 for (
int i =0; i<trial_fes.
Size(); i++)
826 trial_fes[i]->Update(
false);
881 vector<complex<double>> grad;
883 for (
unsigned i = 0; i < grad.size(); i++)
885 grad_r[i] = grad[i].real();
892 vector<complex<double>> grad;
894 for (
unsigned i = 0; i < grad.size(); i++)
896 grad_i[i] = grad[i].imag();
967 complex<double> zi = complex<double>(0., 1.);
983 double alpha = (180+degrees) * M_PI/180.;
984 double sina = sin(
alpha);
985 double cosa = cos(
alpha);
988 double xprim=X(0) + shift;
989 double yprim=X(1) + shift;
991 double x = xprim*sina - yprim*cosa;
992 double y = xprim*cosa + yprim*sina;
994 double rl = 2.*M_PI/rk;
1000 double fact = rl/M_PI/(w0*w0);
1001 double aux = 1. + (fact*y)*(fact*y);
1003 double w = w0*sqrt(aux);
1005 double phi0 = atan(fact*y);
1007 double r = y + 1./y/(fact*fact);
1010 complex<double> ze = - x*x/(w*w) - zi*rk*y - zi * M_PI * x * x/rl/r +
1012 double pf = pow(2.0/M_PI/(w*w),0.25);
1019 double x = X(0)-0.5;
1020 double y = X(1)-0.5;
1021 double r = sqrt(x*x + y*y);
1023 complex<double> Ho = jn(0,
beta) + zi * yn(0,
beta);
1028 MFEM_ABORT(
"Should be unreachable");
1036 dp.resize(X.
Size());
1037 complex<double> zi = complex<double>(0., 1.);
1039 for (
int i = 0; i<X.
Size(); i++) { dp[i] = 0.0; }
1047 complex<double>
p = exp(
alpha);
1048 for (
int i = 0; i<X.
Size(); i++)
1050 dp[i] = zi *
beta *
p;
1058 double degrees = 45;
1059 double alpha = (180+degrees) * M_PI/180.;
1060 double sina = sin(
alpha);
1061 double cosa = cos(
alpha);
1064 double xprim=X(0) + shift;
1065 double yprim=X(1) + shift;
1067 double x = xprim*sina - yprim*cosa;
1068 double y = xprim*cosa + yprim*sina;
1069 double dxdxprim = sina, dxdyprim = -cosa;
1070 double dydxprim = cosa, dydyprim = sina;
1072 double rl = 2.*M_PI/rk;
1078 double fact = rl/M_PI/(w0*w0);
1079 double aux = 1. + (fact*y)*(fact*y);
1081 double w = w0*sqrt(aux);
1082 double dwdy = w0*fact*fact*y/sqrt(aux);
1084 double phi0 = atan(fact*y);
1085 double dphi0dy = cos(phi0)*cos(phi0)*fact;
1087 double r = y + 1./y/(fact*fact);
1088 double drdy = 1. - 1./(y*y)/(fact*fact);
1091 complex<double> ze = - x*x/(w*w) - zi*rk*y - zi * M_PI * x * x/rl/r +
1094 complex<double> zdedx = -2.*x/(w*w) - 2.*zi*M_PI*x/rl/r;
1095 complex<double> zdedy = 2.*x*x/(w*w*w)*dwdy - zi*rk + zi*M_PI*x*x/rl/
1096 (r*r)*drdy + zi*dphi0dy/2.;
1098 double pf = pow(2.0/M_PI/(w*w),0.25);
1099 double dpfdy = -pow(2./M_PI/(w*w),-0.75)/M_PI/(w*w*w)*dwdy;
1101 complex<double> zp = pf*exp(ze);
1102 complex<double> zdpdx = zp*zdedx;
1103 complex<double> zdpdy = dpfdy*exp(ze)+zp*zdedy;
1105 dp[0] = (zdpdx*dxdxprim + zdpdy*dydxprim);
1106 dp[1] = (zdpdx*dxdyprim + zdpdy*dydyprim);
1110 MFEM_ABORT(
"Should be unreachable");
1117 complex<double> zi = complex<double>(0., 1.);
1132 double degrees = 45;
1133 double alpha = (180+degrees) * M_PI/180.;
1134 double sina = sin(
alpha);
1135 double cosa = cos(
alpha);
1138 double xprim=X(0) + shift;
1139 double yprim=X(1) + shift;
1141 double x = xprim*sina - yprim*cosa;
1142 double y = xprim*cosa + yprim*sina;
1143 double dxdxprim = sina, dxdyprim = -cosa;
1144 double dydxprim = cosa, dydyprim = sina;
1146 double rl = 2.*M_PI/rk;
1152 double fact = rl/M_PI/(w0*w0);
1153 double aux = 1. + (fact*y)*(fact*y);
1155 double w = w0*sqrt(aux);
1156 double dwdy = w0*fact*fact*y/sqrt(aux);
1157 double d2wdydy = w0*fact*fact*(1. - (fact*y)*(fact*y)/aux)/sqrt(aux);
1159 double phi0 = atan(fact*y);
1160 double dphi0dy = cos(phi0)*cos(phi0)*fact;
1161 double d2phi0dydy = -2.*cos(phi0)*sin(phi0)*fact*dphi0dy;
1163 double r = y + 1./y/(fact*fact);
1164 double drdy = 1. - 1./(y*y)/(fact*fact);
1165 double d2rdydy = 2./(y*y*y)/(fact*fact);
1168 complex<double> ze = - x*x/(w*w) - zi*rk*y - zi * M_PI * x * x/rl/r +
1171 complex<double> zdedx = -2.*x/(w*w) - 2.*zi*M_PI*x/rl/r;
1172 complex<double> zdedy = 2.*x*x/(w*w*w)*dwdy - zi*rk + zi*M_PI*x*x/rl/
1173 (r*r)*drdy + zi*dphi0dy/2.;
1174 complex<double> zd2edxdx = -2./(w*w) - 2.*zi*M_PI/rl/r;
1175 complex<double> zd2edxdy = 4.*x/(w*w*w)*dwdy + 2.*zi*M_PI*x/rl/(r*r)*drdy;
1176 complex<double> zd2edydx = zd2edxdy;
1177 complex<double> zd2edydy = -6.*x*x/(w*w*w*w)*dwdy*dwdy + 2.*x*x/
1178 (w*w*w)*d2wdydy - 2.*zi*M_PI*x*x/rl/(r*r*r)*drdy*drdy
1179 + zi*M_PI*x*x/rl/(r*r)*d2rdydy + zi/2.*d2phi0dydy;
1181 double pf = pow(2.0/M_PI/(w*w),0.25);
1182 double dpfdy = -pow(2./M_PI/(w*w),-0.75)/M_PI/(w*w*w)*dwdy;
1183 double d2pfdydy = -1./M_PI*pow(2./M_PI,-0.75)*(-1.5*pow(w,-2.5)
1184 *dwdy*dwdy + pow(w,-1.5)*d2wdydy);
1187 complex<double> zp = pf*exp(ze);
1188 complex<double> zdpdx = zp*zdedx;
1189 complex<double> zdpdy = dpfdy*exp(ze)+zp*zdedy;
1190 complex<double> zd2pdxdx = zdpdx*zdedx + zp*zd2edxdx;
1191 complex<double> zd2pdxdy = zdpdy*zdedx + zp*zd2edxdy;
1192 complex<double> zd2pdydx = dpfdy*exp(ze)*zdedx + zdpdx*zdedy + zp*zd2edydx;
1193 complex<double> zd2pdydy = d2pfdydy*exp(ze) + dpfdy*exp(
1194 ze)*zdedy + zdpdy*zdedy + zp*zd2edydy;
1197 return (zd2pdxdx*dxdxprim + zd2pdydx*dydxprim)*dxdxprim
1198 + (zd2pdxdy*dxdxprim + zd2pdydy*dydxprim)*dydxprim
1199 + (zd2pdxdx*dxdyprim + zd2pdydx*dydyprim)*dxdyprim
1200 + (zd2pdxdy*dxdyprim + zd2pdydy*dydyprim)*dydyprim;
1204 MFEM_ABORT(
"Should be unreachable");
1215 for (
int i = 0; i <
dim; ++i)
1217 r += pow(x[i] - center[i], 2.);
1219 double n = 5.0 *
omega / M_PI;
1220 double coeff = pow(n, 2) / M_PI;
1221 double alpha = -pow(n, 2) * r;
Class for domain integration L(v) := (f, v)
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Conjugate gradient method.
virtual Operator & real()
Real or imaginary part accessor methods.
Class for an integration rule - an Array of IntegrationPoint.
The Auxiliary-space Maxwell Solver in hypre.
double source_function(const Vector &x)
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
The Auxiliary-space Divergence Solver in hypre.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void SetDataFormat(VTKFormat fmt)
A coefficient that is constant across space and time.
double rhs_func_r(const Vector &x)
void PrintOptions(std::ostream &out) const
Print the options.
int Dimension() const
Dimension of the reference space used within the elements.
double d2_exact_r(const Vector &x)
void SetSize(int s)
Resize the vector to size s.
Mimic the action of a complex operator using two real operators.
Helper class for ParaView visualization data.
int NumRowBlocks() const
Return the number of row blocks.
void hatu_exact_i(const Vector &X, Vector &hatu)
void PrintUsage(std::ostream &out) const
Print the usage message.
Pointer to an Operator of a specified type.
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
double divu_exact_r(const Vector &x)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
int Size() const
Returns the size of the vector.
complex< double > acoustics_solution_laplacian(const Vector &X)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void gradp_exact_i(const Vector &x, Vector &gradu)
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
void Jt_J_detJinv_i_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Matrix coefficient defined as a product of a scalar coefficient and a matrix coefficient.
void SetOmega(double omega_)
complex< double > acoustics_solution(const Vector &X)
(Q div u, div v) for RT elements
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
double detJ_i_function(const Vector &x, CartesianPML *pml)
Operator & GetDiagonalBlock(int iblock)
Return a reference to block i,i.
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
The BoomerAMG solver in hypre.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
int Append(const T &el)
Append element 'el' to array, resize if necessary.
void EnsureNCMesh(bool simplices_nonconforming=false)
Derived coefficient that takes the value of the parent coefficient for the active attributes and is z...
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
void SetPrintLevel(int print_level)
void abs_Jt_J_detJinv_2_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
void SetMaxIter(int max_it)
double Sum() const
Return the sum of the vector entries.
void SetAttributes(Mesh *mesh_, Array< int > *attrNonPML=nullptr, Array< int > *attrPML=nullptr)
Mark element in the PML region.
double abs_detJ_2_function(const Vector &x, CartesianPML *pml)
void u_exact_r(const Vector &x, Vector &u)
void SetHighOrderOutput(bool high_order_output_)
void VisualizeField(socketstream &sock, const char *vishost, int visport, ParGridFunction &gf, const char *title, int x=0, int y=0, int w=400, int h=400, bool vec=false)
void SetTime(double t)
Set physical time (for time-dependent simulations)
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
void gradp_exact_r(const Vector &x, Vector &gradu)
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
double d2_exact_i(const Vector &x)
A general vector function coefficient.
void u_exact_i(const Vector &x, Vector &u)
Geometry::Type GetElementGeometry(int i) const
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
double p(const Vector &x, double t)
void SetRelTol(double rtol)
double hatp_exact_r(const Vector &X)
double hatp_exact_i(const Vector &X)
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr)
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
double divu_exact_i(const Vector &x)
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
virtual Operator & imag()
double p_exact_i(const Vector &x)
double Max() const
Returns the maximal element of the vector.
Derived matrix coefficient that has the value of the parent matrix coefficient where it is active and...
int GetNE() const
Returns number of elements.
double detJ_r_function(const Vector &x, CartesianPML *pml)
PML stretching functions: See https://doi.org/10.1006/jcph.1994.1159.
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Class for setting up a simple Cartesian PML region.
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
double rhs_func_i(const Vector &x)
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
void Jt_J_detJinv_r_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
double Norml2() const
Returns the l2 norm of the vector.
Abstract class for hypre's solvers and preconditioners.
int Size() const
Return the logical size of the array.
T & Last()
Return the last element in the array.
void SetLevelsOfDetail(int levels_of_detail_)
void acoustics_solution_grad(const Vector &X, vector< complex< double >> &dp)
void Clear()
Clear the contents of the Mesh.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
A general function coefficient.
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
void SetSystemsOptions(int dim, bool order_bynodes=false)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
void hatu_exact_r(const Vector &X, Vector &hatu)
Arbitrary order H1-conforming (continuous) finite elements.
virtual void Save() override
double u(const Vector &xvec)
Class for parallel grid function.
Wrapper for hypre's ParCSR matrix class.
A class to handle Block systems in a matrix-free implementation.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Class for parallel meshes.
double p_exact_r(const Vector &x)
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
int main(int argc, char *argv[])
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Arbitrary order "L2-conforming" discontinuous finite elements.