MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex5.cpp
Go to the documentation of this file.
1// MFEM Example 5
2//
3// Compile with: make ex5
4//
5// Sample runs: ex5 -m ../data/square-disc.mesh
6// ex5 -m ../data/star.mesh
7// ex5 -m ../data/star.mesh -pa
8// ex5 -m ../data/beam-tet.mesh
9// ex5 -m ../data/beam-hex.mesh
10// ex5 -m ../data/beam-hex.mesh -pa
11// ex5 -m ../data/escher.mesh
12// ex5 -m ../data/fichera.mesh
13//
14// Device sample runs:
15// ex5 -m ../data/star.mesh -pa -d cuda
16// ex5 -m ../data/star.mesh -pa -d raja-cuda
17// ex5 -m ../data/star.mesh -pa -d raja-omp
18// ex5 -m ../data/beam-hex.mesh -pa -d cuda
19//
20// Description: This example code solves a simple 2D/3D mixed Darcy problem
21// corresponding to the saddle point system
22//
23// k*u + grad p = f
24// - div u = g
25//
26// with natural boundary condition -p = <given pressure>.
27// Here, we use a given exact solution (u,p) and compute the
28// corresponding r.h.s. (f,g). We discretize with Raviart-Thomas
29// finite elements (velocity u) and piecewise discontinuous
30// polynomials (pressure p).
31//
32// The example demonstrates the use of the BlockOperator class, as
33// well as the collective saving of several grid functions in
34// VisIt (visit.llnl.gov) and ParaView (paraview.org) formats.
35//
36// We recommend viewing examples 1-4 before viewing this example.
37
38#include "mfem.hpp"
39#include <fstream>
40#include <iostream>
41#include <algorithm>
42
43using namespace std;
44using namespace mfem;
45
46// Define the analytical solution and forcing terms / boundary conditions
47void uFun_ex(const Vector & x, Vector & u);
48real_t pFun_ex(const Vector & x);
49void fFun(const Vector & x, Vector & f);
50real_t gFun(const Vector & x);
51real_t f_natural(const Vector & x);
52
53int main(int argc, char *argv[])
54{
55 StopWatch chrono;
56
57 // 1. Parse command-line options.
58 const char *mesh_file = "../data/star.mesh";
59 int order = 1;
60 bool pa = false;
61 const char *device_config = "cpu";
62 bool visualization = 1;
63
64 OptionsParser args(argc, argv);
65 args.AddOption(&mesh_file, "-m", "--mesh",
66 "Mesh file to use.");
67 args.AddOption(&order, "-o", "--order",
68 "Finite element order (polynomial degree).");
69 args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
70 "--no-partial-assembly", "Enable Partial Assembly.");
71 args.AddOption(&device_config, "-d", "--device",
72 "Device configuration string, see Device::Configure().");
73 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
74 "--no-visualization",
75 "Enable or disable GLVis visualization.");
76 args.Parse();
77 if (!args.Good())
78 {
79 args.PrintUsage(cout);
80 return 1;
81 }
82 args.PrintOptions(cout);
83
84 // 2. Enable hardware devices such as GPUs, and programming models such as
85 // CUDA, OCCA, RAJA and OpenMP based on command line options.
86 Device device(device_config);
87 device.Print();
88
89 // 3. Read the mesh from the given mesh file. We can handle triangular,
90 // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
91 // the same code.
92 Mesh *mesh = new Mesh(mesh_file, 1, 1);
93 int dim = mesh->Dimension();
94
95 // 4. Refine the mesh to increase the resolution. In this example we do
96 // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
97 // largest number that gives a final mesh with no more than 10,000
98 // elements.
99 {
100 int ref_levels =
101 (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
102 for (int l = 0; l < ref_levels; l++)
103 {
104 mesh->UniformRefinement();
105 }
106 }
107
108 // 5. Define a finite element space on the mesh. Here we use the
109 // Raviart-Thomas finite elements of the specified order.
110 FiniteElementCollection *hdiv_coll(new RT_FECollection(order, dim));
111 FiniteElementCollection *l2_coll(new L2_FECollection(order, dim));
112
113 FiniteElementSpace *R_space = new FiniteElementSpace(mesh, hdiv_coll);
114 FiniteElementSpace *W_space = new FiniteElementSpace(mesh, l2_coll);
115
116 // 6. Define the BlockStructure of the problem, i.e. define the array of
117 // offsets for each variable. The last component of the Array is the sum
118 // of the dimensions of each block.
119 Array<int> block_offsets(3); // number of variables + 1
120 block_offsets[0] = 0;
121 block_offsets[1] = R_space->GetVSize();
122 block_offsets[2] = W_space->GetVSize();
123 block_offsets.PartialSum();
124
125 std::cout << "***********************************************************\n";
126 std::cout << "dim(R) = " << block_offsets[1] - block_offsets[0] << "\n";
127 std::cout << "dim(W) = " << block_offsets[2] - block_offsets[1] << "\n";
128 std::cout << "dim(R+W) = " << block_offsets.Last() << "\n";
129 std::cout << "***********************************************************\n";
130
131 // 7. Define the coefficients, analytical solution, and rhs of the PDE.
132 ConstantCoefficient k(1.0);
133
137
140
141 // 8. Allocate memory (x, rhs) for the analytical solution and the right hand
142 // side. Define the GridFunction u,p for the finite element solution and
143 // linear forms fform and gform for the right hand side. The data
144 // allocated by x and rhs are passed as a reference to the grid functions
145 // (u,p) and the linear forms (fform, gform).
146 MemoryType mt = device.GetMemoryType();
147 BlockVector x(block_offsets, mt), rhs(block_offsets, mt);
148
149 LinearForm *fform(new LinearForm);
150 fform->Update(R_space, rhs.GetBlock(0), 0);
153 fform->Assemble();
154 fform->SyncAliasMemory(rhs);
155
156 LinearForm *gform(new LinearForm);
157 gform->Update(W_space, rhs.GetBlock(1), 0);
159 gform->Assemble();
160 gform->SyncAliasMemory(rhs);
161
162 // 9. Assemble the finite element matrices for the Darcy operator
163 //
164 // D = [ M B^T ]
165 // [ B 0 ]
166 // where:
167 //
168 // M = \int_\Omega k u_h \cdot v_h d\Omega u_h, v_h \in R_h
169 // B = -\int_\Omega \div u_h q_h d\Omega u_h \in R_h, q_h \in W_h
170 BilinearForm *mVarf(new BilinearForm(R_space));
171 MixedBilinearForm *bVarf(new MixedBilinearForm(R_space, W_space));
172
173 if (pa) { mVarf->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
175 mVarf->Assemble();
176 if (!pa) { mVarf->Finalize(); }
177
178 if (pa) { bVarf->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
180 bVarf->Assemble();
181 if (!pa) { bVarf->Finalize(); }
182
183 BlockOperator darcyOp(block_offsets);
184
185 TransposeOperator *Bt = NULL;
186
187 if (pa)
188 {
189 Bt = new TransposeOperator(bVarf);
190
191 darcyOp.SetBlock(0,0, mVarf);
192 darcyOp.SetBlock(0,1, Bt, -1.0);
193 darcyOp.SetBlock(1,0, bVarf, -1.0);
194 }
195 else
196 {
197 SparseMatrix &M(mVarf->SpMat());
198 SparseMatrix &B(bVarf->SpMat());
199 B *= -1.;
200 Bt = new TransposeOperator(&B);
201
202 darcyOp.SetBlock(0,0, &M);
203 darcyOp.SetBlock(0,1, Bt);
204 darcyOp.SetBlock(1,0, &B);
205 }
206
207 // 10. Construct the operators for preconditioner
208 //
209 // P = [ diag(M) 0 ]
210 // [ 0 B diag(M)^-1 B^T ]
211 //
212 // Here we use Symmetric Gauss-Seidel to approximate the inverse of the
213 // pressure Schur Complement
214 SparseMatrix *MinvBt = NULL;
215 Vector Md(mVarf->Height());
216
217 BlockDiagonalPreconditioner darcyPrec(block_offsets);
218 Solver *invM, *invS;
219 SparseMatrix *S = NULL;
220
221 if (pa)
222 {
223 mVarf->AssembleDiagonal(Md);
224 auto Md_host = Md.HostRead();
225 Vector invMd(mVarf->Height());
226 for (int i=0; i<mVarf->Height(); ++i)
227 {
228 invMd(i) = 1.0 / Md_host[i];
229 }
230
231 Vector BMBt_diag(bVarf->Height());
233
234 Array<int> ess_tdof_list; // empty
235
236 invM = new OperatorJacobiSmoother(Md, ess_tdof_list);
237 invS = new OperatorJacobiSmoother(BMBt_diag, ess_tdof_list);
238 }
239 else
240 {
241 SparseMatrix &M(mVarf->SpMat());
242 M.GetDiag(Md);
244
245 SparseMatrix &B(bVarf->SpMat());
246 MinvBt = Transpose(B);
247
248 for (int i = 0; i < Md.Size(); i++)
249 {
250 MinvBt->ScaleRow(i, 1./Md(i));
251 }
252
253 S = Mult(B, *MinvBt);
254
255 invM = new DSmoother(M);
256
257#ifndef MFEM_USE_SUITESPARSE
258 invS = new GSSmoother(*S);
259#else
260 invS = new UMFPackSolver(*S);
261#endif
262 }
263
264 invM->iterative_mode = false;
265 invS->iterative_mode = false;
266
267 darcyPrec.SetDiagonalBlock(0, invM);
268 darcyPrec.SetDiagonalBlock(1, invS);
269
270 // 11. Solve the linear system with MINRES.
271 // Check the norm of the unpreconditioned residual.
272 int maxIter(1000);
273 real_t rtol(1.e-6);
274 real_t atol(1.e-10);
275
276 chrono.Clear();
277 chrono.Start();
278 MINRESSolver solver;
279 solver.SetAbsTol(atol);
280 solver.SetRelTol(rtol);
281 solver.SetMaxIter(maxIter);
282 solver.SetOperator(darcyOp);
283 solver.SetPreconditioner(darcyPrec);
284 solver.SetPrintLevel(1);
285 x = 0.0;
286 solver.Mult(rhs, x);
287 if (device.IsEnabled()) { x.HostRead(); }
288 chrono.Stop();
289
290 if (solver.GetConverged())
291 {
292 std::cout << "MINRES converged in " << solver.GetNumIterations()
293 << " iterations with a residual norm of "
294 << solver.GetFinalNorm() << ".\n";
295 }
296 else
297 {
298 std::cout << "MINRES did not converge in " << solver.GetNumIterations()
299 << " iterations. Residual norm is " << solver.GetFinalNorm()
300 << ".\n";
301 }
302 std::cout << "MINRES solver took " << chrono.RealTime() << "s.\n";
303
304 // 12. Create the grid functions u and p. Compute the L2 error norms.
306 u.MakeRef(R_space, x.GetBlock(0), 0);
307 p.MakeRef(W_space, x.GetBlock(1), 0);
308
309 int order_quad = max(2, 2*order+1);
311 for (int i=0; i < Geometry::NumGeom; ++i)
312 {
313 irs[i] = &(IntRules.Get(i, order_quad));
314 }
315
316 real_t err_u = u.ComputeL2Error(ucoeff, irs);
317 real_t norm_u = ComputeLpNorm(2., ucoeff, *mesh, irs);
318 real_t err_p = p.ComputeL2Error(pcoeff, irs);
319 real_t norm_p = ComputeLpNorm(2., pcoeff, *mesh, irs);
320
321 std::cout << "|| u_h - u_ex || / || u_ex || = " << err_u / norm_u << "\n";
322 std::cout << "|| p_h - p_ex || / || p_ex || = " << err_p / norm_p << "\n";
323
324 // 13. Save the mesh and the solution. This output can be viewed later using
325 // GLVis: "glvis -m ex5.mesh -g sol_u.gf" or "glvis -m ex5.mesh -g
326 // sol_p.gf".
327 {
328 ofstream mesh_ofs("ex5.mesh");
329 mesh_ofs.precision(8);
330 mesh->Print(mesh_ofs);
331
332 ofstream u_ofs("sol_u.gf");
333 u_ofs.precision(8);
334 u.Save(u_ofs);
335
336 ofstream p_ofs("sol_p.gf");
337 p_ofs.precision(8);
338 p.Save(p_ofs);
339 }
340
341 // 14. Save data in the VisIt format
342 VisItDataCollection visit_dc("Example5", mesh);
343 visit_dc.RegisterField("velocity", &u);
344 visit_dc.RegisterField("pressure", &p);
345 visit_dc.Save();
346
347 // 15. Save data in the ParaView format
348 ParaViewDataCollection paraview_dc("Example5", mesh);
349 paraview_dc.SetPrefixPath("ParaView");
350 paraview_dc.SetLevelsOfDetail(order);
351 paraview_dc.SetCycle(0);
352 paraview_dc.SetDataFormat(VTKFormat::BINARY);
353 paraview_dc.SetHighOrderOutput(true);
354 paraview_dc.SetTime(0.0); // set the time
355 paraview_dc.RegisterField("velocity",&u);
356 paraview_dc.RegisterField("pressure",&p);
357 paraview_dc.Save();
358
359 // 16. Send the solution by socket to a GLVis server.
360 if (visualization)
361 {
362 char vishost[] = "localhost";
363 int visport = 19916;
365 u_sock.precision(8);
366 u_sock << "solution\n" << *mesh << u << "window_title 'Velocity'" << endl;
368 p_sock.precision(8);
369 p_sock << "solution\n" << *mesh << p << "window_title 'Pressure'" << endl;
370 }
371
372 // 17. Free the used memory.
373 delete fform;
374 delete gform;
375 delete invM;
376 delete invS;
377 delete S;
378 delete Bt;
379 delete MinvBt;
380 delete mVarf;
381 delete bVarf;
382 delete W_space;
383 delete R_space;
384 delete l2_coll;
385 delete hdiv_coll;
386 delete mesh;
387
388 return 0;
389}
390
391
392void uFun_ex(const Vector & x, Vector & u)
393{
394 real_t xi(x(0));
395 real_t yi(x(1));
396 real_t zi(0.0);
397 if (x.Size() == 3)
398 {
399 zi = x(2);
400 }
401
402 u(0) = - exp(xi)*sin(yi)*cos(zi);
403 u(1) = - exp(xi)*cos(yi)*cos(zi);
404
405 if (x.Size() == 3)
406 {
407 u(2) = exp(xi)*sin(yi)*sin(zi);
408 }
409}
410
411// Change if needed
413{
414 real_t xi(x(0));
415 real_t yi(x(1));
416 real_t zi(0.0);
417
418 if (x.Size() == 3)
419 {
420 zi = x(2);
421 }
422
423 return exp(xi)*sin(yi)*cos(zi);
424}
425
426void fFun(const Vector & x, Vector & f)
427{
428 f = 0.0;
429}
430
432{
433 if (x.Size() == 3)
434 {
435 return -pFun_ex(x);
436 }
437 else
438 {
439 return 0;
440 }
441}
442
444{
445 return (-pFun_ex(x));
446}
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
T & Last()
Return the last element in the array.
Definition array.hpp:802
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
virtual void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag. Note that diag is a tdof Vector.
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
A class to handle Block diagonal preconditioners in a matrix-free implementation.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
A class to handle Block systems in a matrix-free implementation.
void SetBlock(int iRow, int iCol, Operator *op, real_t c=1.0)
Add a block op in the block-entry (iblock, jblock).
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
A coefficient that is constant across space and time.
Data type for scaled Jacobi-type smoother of sparse matrix.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition device.hpp:278
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:286
static bool IsEnabled()
Return true if any backend other than Backend::CPU is enabled.
Definition device.hpp:247
Class for domain integration .
Definition lininteg.hpp:109
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
A general function coefficient.
Data type for Gauss-Seidel smoother of sparse matrix.
static const int NumGeom
Definition geom.hpp:42
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void MakeRef(FiniteElementSpace *f, real_t *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition gridfunc.cpp:203
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
real_t GetFinalNorm() const
Returns the final residual norm after termination of the solver during the last call to Mult().
Definition solvers.hpp:275
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition solvers.hpp:260
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
bool GetConverged() const
Returns true if the last call to Mult() converged successfully.
Definition solvers.hpp:262
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Vector with associated FE space and LinearFormIntegrators.
Adds new Domain Integrator. Assumes ownership of lfi.
Adds new Boundary Integrator. Assumes ownership of lfi.
void Update()
Update the object according to the associated FE space fes.
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
MINRES method.
Definition solvers.hpp:628
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the MINRES method.
Definition solvers.cpp:1616
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:1595
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.hpp:640
Mesh data type.
Definition mesh.hpp:56
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY.
void Assemble(int skip_zeros=1)
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACY.
void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const
Assemble the diagonal of ADA^T into diag, where A is this mixed bilinear form and D is a diagonal.
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
Adds a domain integrator. Assumes ownership of bfi.
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition solvers.hpp:313
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Helper class for ParaView visualization data.
void SetHighOrderOutput(bool high_order_output_)
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
void SetDataFormat(VTKFormat fmt)
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Base class for solvers.
Definition operator.hpp:683
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition operator.hpp:686
Data type sparse matrix.
Definition sparsemat.hpp:51
void GetDiag(Vector &d) const
Returns the Diagonal of A.
void ScaleRow(const int row, const real_t scale)
Timing object.
Definition tic_toc.hpp:36
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition tic_toc.cpp:432
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition tic_toc.cpp:411
void Stop()
Stop the stopwatch.
Definition tic_toc.cpp:422
void Clear()
Clear the elapsed time on the stopwatch and restart it if it's running.
Definition tic_toc.cpp:406
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition operator.hpp:750
Direct sparse solver using UMFPACK.
Definition solvers.hpp:1096
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition lininteg.hpp:347
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
virtual const real_t * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:478
void SyncAliasMemory(const Vector &v) const
Update the alias memory location of the vector to match v.
Definition vector.hpp:259
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:494
Data collection with VisIt I/O routines.
virtual void Save() override
Save the collection and a VisIt root file.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
int dim
Definition ex24.cpp:53
real_t pFun_ex(const Vector &x)
Definition ex5.cpp:412
void uFun_ex(const Vector &x, Vector &u)
Definition ex5.cpp:392
real_t gFun(const Vector &x)
Definition ex5.cpp:431
real_t f_natural(const Vector &x)
Definition ex5.cpp:443
void fFun(const Vector &x, Vector &f)
Definition ex5.cpp:426
int main()
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:476
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition table.cpp:414
real_t ComputeLpNorm(real_t p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
Compute the Lp norm of a function f. .
float real_t
Definition config.hpp:43
MemoryType
Memory types supported by MFEM.
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30