49 Data(
double x_,
double val_) {x=x_; val=val_;};
52 inline bool operator==(
const Data& d1,
const Data& d2) {
return (d1.x == d2.x); };
53 inline bool operator <(
const Data& d1,
const Data& d2) {
return (d1.x < d2.x); };
60 #ifndef MFEM_THREAD_SAFE 67 Diffusion2Integrator() { Q = NULL; }
82 #ifdef MFEM_THREAD_SAFE 95 if (el.
Space() == FunctionSpace::Pk)
104 if (el.
Space() == FunctionSpace::rQk)
126 w *= Q->
Eval(Trans, ip);
129 for (
int jj = 0; jj < nd; jj++)
131 for (
int ii = 0; ii < nd; ii++)
133 elmat(ii, jj) += w*shape(ii)*laplace(jj);
141 int main(
int argc,
char *argv[])
144 const char *mesh_file =
"../../data/star.mesh";
145 const char *per_file =
"none";
146 const char *ref_file =
"";
150 bool static_cond =
false;
151 bool visualization = 1;
160 args.
AddOption(&mesh_file,
"-m",
"--mesh",
161 "Mesh file to use.");
162 args.
AddOption(&ref_levels,
"-r",
"--refine",
163 "Number of times to refine the mesh uniformly, -1 for auto.");
165 "Periodic BCS file.");
166 args.
AddOption(&ref_file,
"-rf",
"--ref-file",
167 "File with refinement data");
168 args.
AddOption(&master,
"-pm",
"--master",
169 "Master boundaries for periodic BCs");
171 "Slave boundaries for periodic BCs");
173 "Finite element order (polynomial degree) or -1 for" 174 " isoparametric space.");
175 args.
AddOption(&ibp,
"-ibp",
"--ibp",
"-no-ibp",
177 "Selects the standard weak form (IBP) or the nonstandard (NO-IBP).");
178 args.
AddOption(&strongBC,
"-sbc",
"--strong-bc",
"-wbc",
180 "Selects strong or weak enforcement of Dirichlet BCs.");
182 "Sets the SIPG penalty parameters, should be positive." 183 " Negative values are replaced with (order+1)^2.");
184 args.
AddOption(&static_cond,
"-sc",
"--static-condensation",
"-no-sc",
185 "--no-static-condensation",
"Enable static condensation.");
186 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
187 "--no-visualization",
188 "Enable or disable GLVis visualization.");
189 args.
AddOption(&lod,
"-lod",
"--level-of-detail",
190 "Refinement level for 1D solution output (0 means no output).");
197 if (!strongBC & (
kappa < 0))
206 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
215 if (mesh->
NURBSext && (strlen(ref_file) != 0))
223 (int)floor(log(5000./mesh->
GetNE())/log(2.)/
dim);
226 for (
int l = 0; l < ref_levels; l++)
246 if (order.
Size() == 1)
253 if (order.
Size() != nkv ) {
mfem_error(
"Wrong number of orders set."); }
258 in.open(per_file, std::ifstream::in);
265 master.
Load(in, psize);
266 slave.
Load(in, psize);
273 else if (order[0] == -1)
279 cout <<
"Using isoparametric FEs: " << fec->
Name() << endl;
283 cout <<
"Mesh does not have FEs --> Assume order 1.\n";
290 if (order.
Size() > 1) { cout <<
"Wrong number of orders set, needs one.\n"; }
296 cout <<
"Number of finite element unknowns: " 303 cout <<
"No integration by parts requires a NURBS mesh."<< endl;
308 cout <<
"No integration by parts requires a NURBS mesh, with only 1 patch."<<
310 cout <<
"A C_1 discretisation is required."<< endl;
311 cout <<
"Currently only C_0 multipatch coupling implemented."<< endl;
316 cout <<
"No integration by parts requires at least quadratic NURBS."<< endl;
317 cout <<
"A C_1 discretisation is required."<< endl;
340 for (
int i = 0; i < master.
Size(); i++)
342 ess_bdr[master[i]-1] = 0;
343 ess_bdr[slave[i]-1] = 0;
357 b->AddBdrFaceIntegrator(
377 a->AddDomainIntegrator(
new Diffusion2Integrator(one));
389 if (static_cond) {
a->EnableStaticCondensation(); }
394 a->FormLinearSystem(ess_tdof_list, x, *
b, A, X, B);
396 cout <<
"Size of linear system: " << A.
Height() << endl;
398 #ifndef MFEM_USE_SUITESPARSE 402 PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
406 umf_solver.
Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
408 umf_solver.
Mult(B, X);
412 a->RecoverFEMSolution(X, *
b, x);
417 ofstream mesh_ofs(
"refined.mesh");
418 mesh_ofs.precision(8);
419 mesh->
Print(mesh_ofs);
420 ofstream sol_ofs(
"sol.gf");
421 sol_ofs.precision(8);
432 sol_sock.precision(8);
433 sol_sock <<
"solution\n" << *mesh << x << flush;
448 for (
int i = 0; i < mesh->
GetNE(); i++)
457 for (
int j = 0; j < vals.
Size(); j++)
459 sol.push_back(Data(coords[j],vals[j]));
464 ofstream sol_ofs(
"solution.dat");
465 for (std::list<Data>::iterator d = sol.begin(); d != sol.end(); ++d)
467 sol_ofs<<d->x <<
"\t"<<d->val<<endl;
482 if (own_fec) {
delete fec; }
Abstract class for all finite elements.
Arbitrary order non-uniform rational B-splines (NURBS) finite elements.
Class for domain integration L(v) := (f, v)
int GetNPoints() const
Returns the number of the points in the integration rule.
void Load(std::istream &in, int fmt=0)
Read an Array from the stream in using format fmt. The format fmt can be:
Class for an integration rule - an Array of IntegrationPoint.
Class for grid function - Vector with associated FE space.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
A coefficient that is constant across space and time.
void PrintOptions(std::ostream &out) const
Print the options.
Geometry::Type GetElementBaseGeometry(int i) const
virtual void PrintInfo(std::ostream &os=mfem::out)
In serial, this method calls PrintCharacteristics(). In parallel, additional information about the pa...
int Dimension() const
Dimension of the reference space used within the elements.
void SetSize(int s)
Resize the vector to size s.
int Space() const
Returns the type of FunctionSpace on the element.
void PrintUsage(std::ostream &out) const
Print the usage message.
virtual void CalcPhysLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in reference space at the gi...
int Size() const
Returns the size of the vector.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Data type dense matrix using column-major storage.
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
IntegrationRules RefinedIntRules(1, Quadrature1D::GaussLegendre)
A global object with all refined integration rules.
bool Good() const
Return true if the command line options were parsed successfully.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
GeometryRefiner GlobGeometryRefiner
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Data collection with VisIt I/O routines.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
int main(int argc, char *argv[])
virtual const char * Name() const
void Print(std::ostream &out=mfem::out, int width=4) const
Prints array to stream with width elements per row.
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
bool operator==(const Data &d1, const Data &d2)
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
double Control[UMFPACK_CONTROL]
int GetDim() const
Returns the reference space dimension for the finite element.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
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...
virtual void Save() override
Save the collection and a VisIt root file.
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.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
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().
void RefineNURBSFromFile(std::string ref_file)
int GetNE() const
Returns number of elements.
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Class for integration point with weight.
bool operator<(const Data &d1, const Data &d2)
int Size() const
Return the logical size of the array.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
virtual void Print(std::ostream &os=mfem::out) const
Arbitrary order H1-conforming (continuous) finite elements.
void GetNodes(Vector &node_coord) const
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.