mfem/mfem

Lightweight, general, scalable C++ library for finite element methods

                Finite Element Discretization Library
                               __
                   _ __ ___   / _|  ___  _ __ ___
                  | '_ ` _ \ | |_  / _ \| '_ ` _ \
                  | | | | | ||  _||  __/| | | | | |
                  |_| |_| |_||_|   \___||_| |_| |_|

                           https://mfem.org

MFEM is a modular parallel C++ library for finite element methods. Its goal is to enable high-performance scalable finite element discretization research and application development on a wide variety of platforms, ranging from laptops to supercomputers.

We welcome contributions and feedback from the community. Please see the file CONTRIBUTING.md for additional details about our development process.

• For building instructions, see the file INSTALL, or type "make help".

• Copyright and licensing information can be found in files LICENSE and NOTICE.

• The best starting point for new users interested in MFEM's features is to review the examples and miniapps at https://mfem.org/examples.

• Instructions for learning with Docker are in config/docker.

Conceptually, MFEM can be viewed as a finite element toolbox that provides the building blocks for developing finite element algorithms in a manner similar to that of MATLAB for linear algebra methods. In particular, MFEM provides support for arbitrary high-order H1-conforming, discontinuous (L2), H(div)-conforming, H(curl)-conforming and NURBS finite element spaces in 2D and 3D, as well as many bilinear, linear and nonlinear forms defined on them. It enables the quick prototyping of various finite element discretizations, including Galerkin methods, mixed finite elements, Discontinuous Galerkin (DG), isogeometric analysis, hybridization and Discontinuous Petrov-Galerkin (DPG) approaches.

MFEM includes classes for dealing with a wide range of mesh types: triangular, quadrilateral, tetrahedral and hexahedral, as well as surface and topologically periodical meshes. It has general support for mesh refinement, including local conforming and non-conforming (AMR) adaptive refinement. Arbitrary element transformations, allowing for high-order mesh elements with curved boundaries, are also supported.

When used as a "finite element to linear algebra translator", MFEM can take a problem described in terms of finite element-type objects, and produce the corresponding linear algebra vectors and fully or partially assembled operators, e.g. in the form of global sparse matrices or matrix-free operators. The library includes simple smoothers and Krylov solvers, such as PCG, MINRES and GMRES, as well as support for sequential sparse direct solvers from the SuiteSparse library. Nonlinear solvers (the Newton method), eigensolvers (LOBPCG), and several explicit and implicit Runge-Kutta time integrators are also available.

MFEM supports MPI-based parallelism throughout the library, and can readily be used as a scalable unstructured finite element problem generator. Starting with version 4.0, MFEM offers support for GPU acceleration, and programming models, such as CUDA, HIP, OCCA, RAJA and OpenMP. MFEM-based applications require minimal changes to switch from a serial to a highly-performant MPI-parallel version of the code, where they can take advantage of the integrated linear solvers from the hypre library. Comprehensive support for other external packages, e.g. PETSc, SUNDIALS and libCEED is also included, giving access to additional linear and nonlinear solvers, preconditioners, time integrators, etc.

For examples of using MFEM, see the examples/ and miniapps/ directories, as well as the OpenGL visualization tool GLVis which is available at https://glvis.org.

License

MFEM is distributed under the terms of the BSD-3 license. All new contributions must be made under this license. See LICENSE and NOTICE for details.

SPDX-License-Identifier: BSD-3-Clause
LLNL Release Number: LLNL-CODE-806117
DOI: 10.11578/dc.20171025.1248

mfem's Issues

Nonlinear term calculation

Hi,

I am testing a convection diffusion problem.

The diffusion term is treated implicitly by solving a Helmholtz equation. The nonlinear convection term 0.5*(\partial u_i u_j/\partial x_j+ u_j \partial u_i/\partial x_j) is treated explicitly, The geometry I am testing is a fully-developed plane channel.

Without the convection term, the result is very good. However, when the convection term is turned on, it diverged. It appears that some small error increases uncontrollably due to the nonlinearity. I suspect that a filter/dealiasing approach is necessary.

I use these lines to calculate u_0 \partial u_0/\partial x_0
wrks2[0]->ProjectGridFunction(u[0]); //Project u[0] to 1 order lower FEC wrks2[0]
wrks[0]->ProjectGridFunction(wrks2[0]); //Project back to the original FEC wrks[0]
wrks[0]->GetDerivative(1,0,wrks[1]); // Find derivative \partial u_0/\partial x_0
wrks[1]->Times(wrks[0]); // u_0 \partial u_0/\partial x_0

where I added function Times in GridFunction as

void GridFunction::Times(const GridFunction &v)
{
for (int i = 0; i < size; i++)
data[i] *= v.data[i];
}

The diverge problem persists. Do you think I did anything wrong here?

Thank you,
Frank

MatrixFunctionCoefficient

Hi,
While I use MatrixFunctionCoefficient for k in example 5, I got the following error:

error: call of overloaded ‘MatrixFunctionCoefficient(int&, void (&)(mfem::Matrix&))’ is ambiguous
../fem/../mesh/../fem/coefficient.hpp:386: note: candidates are: mfem::MatrixFunctionCoefficient::MatrixFunctionCoefficient(int, void ()(const mfem::Vector&, double, mfem::DenseMatrix&))
../fem/../mesh/../fem/coefficient.hpp:378: note: mfem::MatrixFunctionCoefficient::MatrixFunctionCoefficient(int, void ()(const mfem::Vector&, mfem::DenseMatrix&))

How do I fix it?
mp

Example 1 help

For EX1, actually it is Delta u = const which is Poisson solver and not Laplace 
solver as in head comment. Now suppose, if I consider that I have a solid rod 
and I want to set outer boundary temperature is say 1000.0 degrees. I want to 
see distribution radially. 

To do this, shall I change:

ConstantCoefficient Temp(1000.0);

and change b->AddDomainIntegrator in LinearForm and a->AddDomainIntegrator in 
Bilinear to incorporate Temp instead of one?

Thanks and regs

Original issue reported on code.google.com by [email protected] on 30 Jul 2013 at 11:37

Wrapper to the hypre eigensolver

Adding a wrapper to the hypre eigensolver and a corresponding modification of 
ex3/ex3p, please? 

Andrew Knyazev

Original issue reported on code.google.com by [email protected] on 9 Aug 2012 at 9:56

DGElasticityIntegrator

Dear All,

There is a DGDiffusionIntegrator that complements a DiffusionIntegrator for implementation of a DG method (SIPG, NIPG, etc). I wonder how difficult it would be to create a DGElasticityIntegrator that being added to ElasticityIntegrator could be used for implementation of a DG method (same SIPG, NIPG) for problems of linear elasticity (I would use it to solve an elastic wave equation in a fractured domain, because fractures can be naturally taken into account in DG).

Thanks,
Mikhail

Dofs vs True dofs in Parallel classes

Dear all,

I don't quite understand what is the difference between dofs and true dofs when a problem is being solved in parallel. What makes true dofs true? Why do we need them?

Thank you,
Mikhail

error message

Hi,
I'm testing example 5 by changing the coefficient k to a function.

The error message is as follows:
ex25.cpp: In function ‘int main(int, char*)’:
ex25.cpp:159: error: no matching function for call to ‘mfem::VectorFEDomainLFIntegrator::VectorFEDomainLFIntegrator(mfem::FunctionCoefficient&)’
../fem/lininteg.hpp:179: note: candidates are: mfem::VectorFEDomainLFIntegrator::VectorFEDomainLFIntegrator(mfem::VectorCoefficient&)
../fem/lininteg.hpp:172: note: mfem::VectorFEDomainLFIntegrator::VectorFEDomainLFIntegrator(const mfem::VectorFEDomainLFIntegrator&)
make: ** [ex25] Error 1

Can you explain how to fix it?
mp

Report for Dev release

In OpenSuse 12.3, I tried to run this dev series while I got following error 
which may help you to fix the development issues.
Best,
Abhijit

make
cd lib; g++  -O3 -I/usr/X11R6/include  -DGLVIS_USE_LIBPNG -DGLVIS_USE_FREETYPE 
-I/usr/X11R6/include/freetype2 -I/usr/include/freetype2 -DGLVIS_MULTISAMPLE=4 
-DGLVIS_MS_LINEWIDTH=1.4 -I../../mfem -c aux_gl.cpp
cd lib; g++  -O3 -I/usr/X11R6/include  -DGLVIS_USE_LIBPNG -DGLVIS_USE_FREETYPE 
-I/usr/X11R6/include/freetype2 -I/usr/include/freetype2 -DGLVIS_MULTISAMPLE=4 
-DGLVIS_MS_LINEWIDTH=1.4 -I../../mfem -c aux_vis.cpp
cd lib; gcc  -O3 -I/usr/X11R6/include -I../../mfem -c gl2ps.c
cd lib; g++  -O3 -I/usr/X11R6/include  -DGLVIS_USE_LIBPNG -DGLVIS_USE_FREETYPE 
-I/usr/X11R6/include/freetype2 -I/usr/include/freetype2 -DGLVIS_MULTISAMPLE=4 
-DGLVIS_MS_LINEWIDTH=1.4 -I../../mfem -c material.cpp
cd lib; g++  -O3 -I/usr/X11R6/include  -DGLVIS_USE_LIBPNG -DGLVIS_USE_FREETYPE 
-I/usr/X11R6/include/freetype2 -I/usr/include/freetype2 -DGLVIS_MULTISAMPLE=4 
-DGLVIS_MS_LINEWIDTH=1.4 -I../../mfem -c openglvis.cpp
cd lib; g++  -O3 -I/usr/X11R6/include  -DGLVIS_USE_LIBPNG -DGLVIS_USE_FREETYPE 
-I/usr/X11R6/include/freetype2 -I/usr/include/freetype2 -DGLVIS_MULTISAMPLE=4 
-DGLVIS_MS_LINEWIDTH=1.4 -I../../mfem -c threads.cpp
threads.cpp: In static member function ‘static void* 
communication_thread::execute(void*)’:
threads.cpp:560:73: error: invalid user-defined conversion from 
‘std::basic_istream<char>’ to ‘int’ [-fpermissive]
In file included from /usr/include/c++/4.7/ios:45:0,
                 from /usr/include/c++/4.7/ostream:40,
                 from /usr/include/c++/4.7/iostream:40,
                 from ../../mfem/fem/../general/array.hpp:15,
                 from ../../mfem/fem/fem.hpp:15,
                 from ../../mfem/mfem.hpp:13,
                 from vsdata.hpp:16,
                 from visual.hpp:20,
                 from threads.cpp:16:
/usr/include/c++/4.7/bits/basic_ios.h:113:7: note: candidate is: 
std::basic_ios<_CharT, _Traits>::operator void*() const [with _CharT = char; 
_Traits = std::char_traits<char>] <near match>
/usr/include/c++/4.7/bits/basic_ios.h:113:7: note:   no known conversion for 
implicit ‘this’ parameter from ‘void*’ to ‘int’
threads.cpp:580:72: error: invalid user-defined conversion from ‘std::istream 
{aka std::basic_istream<char>}’ to ‘int’ [-fpermissive]
In file included from /usr/include/c++/4.7/ios:45:0,
                 from /usr/include/c++/4.7/ostream:40,
                 from /usr/include/c++/4.7/iostream:40,
                 from ../../mfem/fem/../general/array.hpp:15,
                 from ../../mfem/fem/fem.hpp:15,
                 from ../../mfem/mfem.hpp:13,
                 from vsdata.hpp:16,
                 from visual.hpp:20,
                 from threads.cpp:16:
/usr/include/c++/4.7/bits/basic_ios.h:113:7: note: candidate is: 
std::basic_ios<_CharT, _Traits>::operator void*() const [with _CharT = char; 
_Traits = std::char_traits<char>] <near match>
/usr/include/c++/4.7/bits/basic_ios.h:113:7: note:   no known conversion for 
implicit ‘this’ parameter from ‘void*’ to ‘int’
make: *** [lib/threads.o] Error 1

Original issue reported on code.google.com by [email protected] on 24 Jul 2013 at 4:22

ex5.cpp

Propose moving mfem repo to llnl organization

Hi @tzanio -- I'd like to propose moving the MFEM organization into the LLNL organization. I can create an mfem team for you there to duplicate the mfem organization access.

Additionally, GitHub will handle page redirects from the old organization to the new if we do it as a transfer of the repositories.

Thoughts? Concerns?

Testing mfem-dev hookup

Test

Original issue reported on code.google.com by tzanio on 2 Oct 2013 at 8:52

Any plan to add two operators into mfem?

Dear Tzanio,

I am glad to hear that new version 3.1 is released and the non-conforming mesh is supported completely. I wish the mfem can be applied in electromagnetic calculation. In this field, two operators are often used:

1. flux
   \phi _m = \int_S {N{\bf{B}} \cdot {\bf{\hat n}}dA = } \int_S {NB_n dA}
2. Electromotive Force
   \xi = \oint_C {E \cdot d\ell }

Is there any plan to add these two operators into mfem?

Thanks,
Tang Laoya

Matrix bandwidth optimization

Hi, all,

Will a matrix bandwidth optimization, such as Cuthill–McKee algorithm, make the AMG solver faster?

Thank you,
Frank

Fix the value at one dof

Hi,

I am trying to solve \Delta u=1, with u=0 imposed at one point. This is commonly used in CFD for the Poisson equation of pressure.

I modified "periodic-cube.mesh" to make it only periodic in x and y directions. I want to fix the first dof in process 0 to 0.0. Here is how I tried to do this:

Array pdofs(1); Vector pfix(1);
pdofs=0; pfix=0.0;
if(myid==0) a->EliminateVDofs( pdofs, pfix,*b);

Here is the result.

Could you please explain to me where I did wrong?

Thank you very much,
Frank

How to formulate shell elements?

Hello,

We are currently evaluating mfem with regard to shell elements, i.e. 2 dimensional elements that live in 3 dimensional space. What is the canonical way to describe that problem within mfem?

We've made experience that using 2d elements on a 3d mesh essentially makes the mesh 2d by ignoring the z coordinate.

In our mesh file we have

dimension
3
elements
N
1 2 0 1 2
[...]

In this case mfem just crashes.

Another case we tried is to set dimension to 2 and use element type triangle with 3d vertices, but mfem ignores the third coordinate.

Is there an example for this type of problem?
Is this the right place to ask or is there any mailing list for discussion?

Thanks in advance!

Coordinates of points associated with DOFs

Dear All,

Is there a way to get coordinates of points associated with degrees of freedom? Can I get them for one finite element only, or can I get all of them? For example, I'd like to print the DOF points and the corresponding value of the solution:

DOF (point) | Solution (value)
x | y  | z  | U

Thanks,
Mikhail

Mesh file format for MFEM

I see that NETGEN with OCC can generate mesh which could be exported in 
different file formats. Can you kindly tell me which file format would be 
useful or do I need to add any extra line like MFEM mesh v1.0.... etc?

Original issue reported on code.google.com by [email protected] on 29 Jul 2013 at 6:04

Legendre polynomials as basis for H1 space

Dear all,

I was recently trying to switch from Lagrange basis that is used for H1 to Legendre, and I seem to understand how to do that, but it looks like I need to change MFEM code for that and rebuild the library - some changes in Poly_1D would be needed, and H1_(Segment|Quadrilateral|Hexahedron)Element classes are heavily based on it, so I can't, for example, inherit from Poly_1D and then feed H1_*Element object with an instance of my Poly_1D derived class. One solution would be to pass a pointer to a Poly_1D (or derived class) object to the H1_*Element constructor - the same way it's done for many other functions in MFEM. It doesn't seem hard, but before I start doing that for my project, I am wondering if there is a better solution, and, maybe, I am just wasting time. So the question is - how to get Legendre basis functions in H1 spaces on Segment, Quadrilateral and Hexahedron elements?

Thank you,
Mikhail

Interpolate a grid function at an arbitrary point

Dear all,

I need to compute the value of my solution (a grid function) at a point which may or may not coincide with a dof point, so I need to interpolate it. I didn't find any existing functions for that, so I'm wondering what is best way to do that within MFEM.

Thank you,
Mikhail

How to build mfem in windows 7 64bit + VS2010/VS2012?

Dear all,
Can I use the MFEM in windows 7 64bit + VS2010/VS2012? How to build it?

Thanks,
Tang Laoya

Original issue reported on code.google.com by [email protected] on 22 Dec 2013 at 11:43

Building library for using LAPACK

Hi,
I need to use methods for DenseMatrix Class, then I have to build (serial) 
library by compiling with LAPACK.

I defined the variable USE_LAPACK       = YES in the makefile, but this setting 
was ignored by the compiler, indeed I see the following compilation line when I 
run make:
cd linalg; g++ -O3  -DMFEM_USE_MEMALLOC      -c densemat.cpp


Finally, my code using DenseMatrixEigensystem doesn't work, and the following 
message is obtained:

DenseMatrixEigensystem::Eval(): Compiled without LAPACK

Could you help me?
Thanks a lot, Pasqua D'Ambra

Original issue reported on code.google.com by [email protected] on 27 Jan 2014 at 9:55

Attachments:

• makefile

Setting HYPRE_DIR

Shouldn't HYPRE_DIR be set before MPICC ?

# Parallel compiler
MPICC      = mpicxx
MPIOPTS    = $(CCOPTS) -I$(HYPRE_DIR)/include
MPIDEBUG   = $(DEBUG_OPTS) -I$(HYPRE_DIR)/include

# The HYPRE library (needed to build the parallel version)
HYPRE_DIR  = ../../hypre-2.8.0b/src/hypre

Regards,
Jian

Original issue reported on code.google.com by [email protected] on 16 Feb 2012 at 8:14

ex11

MassIntegrator on boundary of domain in DG space

Dear All,

I'd like to compute a matrix of integrals over boundary of a domain K :
M = \int_{\partial K} u v,
where u and v are from the finite element space spanned over discontinuous functions

Here is what I'm trying to do with MFEM:

int order = 1;
Mesh mesh(10, 10, Element::QUADRILATERAL, 1, 100, 100);
int dim = mesh.Dimension();
DG_FECollection fec(order, dim);
FiniteElementSpace fespace(&mesh, &fec);
BilinearForm bound_mass(&fespace);
// I tried this:
bound_mass.AddBoundaryIntegrator(new MassIntegrator());
// But it gives a segfault, because it dereferences a NULL pointer in assembling:
// <bilinearform.cpp>
// void BilinearForm::Assemble : 
// if (bbfi.Size())
//   {
//      for (i = 0; i < fes -> GetNBE(); i++)
//      {
//         const FiniteElement &be = *fes->GetBE(i); // Here fes->GetBE(i) return NULL pointer
// Apparently it shouldn't work that way, so I tried:
bound_mass.AddBdrFaceIntegrator(new MassIntegrator());
// but it's not implemented.

I've tried to implement void MassIntegrator:AssembleFaceMatrix to make the second approach work, but I'm sure I did it wrong:

void MassIntegrator::AssembleFaceMatrix(
    const FiniteElement &el1, const FiniteElement &el2,
    FaceElementTransformations &Trans, DenseMatrix &elmat)
{
   int nd1 = el1.GetDof();
   double w;

   elmat.SetSize(nd1);
   shape.SetSize(nd1);

   const IntegrationRule *ir = IntRule;
   if (ir == NULL)
   {
      int order = 2 * el1.GetOrder();

      ir = &IntRules.Get(Trans.FaceGeom, order);
   }

   elmat = 0.0;
   for (int i = 0; i < ir->GetNPoints(); i++)
   {
      const IntegrationPoint &ip = ir->IntPoint(i);
      el1.CalcShape(ip, shape);

      Trans.Face->SetIntPoint(&ip);
      w = Trans.Face->Weight() * ip.weight;
      if (Q)
      {
         w *= Q -> Eval(*Trans.Face, ip);
      }

      AddMult_a_VVt(w, shape, elmat);
   }
}

Could you please advise something on that matter?

Thanks,
Mikhail

Right hand side of PDE

Hi,
In MFEM, I saw examples in which right hand sides of Partial Differential Equations (PDE)
are added by imposing exact solutions.
When you compute numerical solutions of coupled PDE in which the solution of first PDE is
the right hand side (as input) of second one, how do you write the "main" function?
For example, Poisson-Nernst-Planck equation is one example of this case.

Thank you.

mp

Map between indices of elements of Mesh and ParMesh

Dear All,

As far as I understand, after I create a ParMesh based on a Mesh, the indices of the elements are rewritten, so that they are local with respect to the process. Am I right? If so, I'm wondering if I can retrieve an element index in the global Mesh based on its local index in ParMesh. The reason why I need this is because for many of my problems I have cell-wise constant coefficients, i.e. they are constant within one grid cell. Therefore, I've created a simple class to work with such coefficients that computes the value of the coefficient based on the ElementNo of a transformation. Now, I faced a situation when I can't use this approach as ElementNo is local w.r.t. each process and the values of the coefficients that go in the assembling procedures are incorrect.

Any hints would be highly appreciated!

Thanks,
Mikhail

File dated before 1/1/1970 can't be built by VS 2015

Dear administrator,

I found that mfem can't be built by VS 2015 for the following errors:
fatal error C1073: Internal error involving incremental compilation

I googled and found the errors:
https://connect.microsoft.com/VisualStudio/feedback/details/1576822/file-from-year-1601-causes-fatal-error-c1073-internal-error-involving-incremental-compilation

Although it is a bug of VS 2015, but if you can update the files' date it could be better.

Thanks,
Tang Laoya

Infinite loop in GeometryRefiner::Refine

If GeometryRefiner::Refine() is called on a mesh with quadrilaterals and with Times=0, ETimes=1, then lines 625 and 634 in geom.cpp lead to loops that do not increment and so go on infinitely.

The call to GeometryRefiner::Refine() with Times=0, ETimes=1 happens if you call SaveVTK() with ref=0, see gridfunc.cpp line 2148. This seems like a reasonable way to call SaveVTK().

suitesparse

The newest version of SuiteSparse did not work for me. I have both SuiteSparse 5.7.1 installed from source and downloaded from the suite sparse website and also SuiteSparse 4.2.1 installed with macports on macbookpro. The mfem configures and compiles (gcc-mp-4.8) with both versions with no problem. The mfem examples as well. Executing "ex1", however works with SuiteSp4.2.1 and does not with SuiteSp5.7.1. The output is below, have a look and see if you can fix.

===================OUTPUT of mfem (SuiteSparse 5.7.1)
[machine~]./ex1 -m ../data/fichera.mesh
Options used:
--mesh ../data/fichera.mesh
--order 1
--visualization
Mesh Characteristics:
Number of vertices : 31841
Number of edges : 92256
Number of faces : 89088
Number of elements : 28672
Number of bdr elem : 6144
Euler Number : 1
h_min : 0.0625
h_max : 0.0625
kappa_min : 1
kappa_max : 1

NumOfBdrVertices=6146 NumberOfBdrElements=6144
Number of unknowns: 31841

UMFPACK V5.7.1 (Oct 10, 2014): ERROR: ordering failed

UMFPackSolver::SetOperator : umfpack_di_symbolic() failed!
Abort trap: 6
==============================end output mfem(suitesparse 5.7.1)
====================OUTPUT mfem(Suitesparse 4.2.1)
[machine ~]./ex1 -m ../data/fichera.mesh
Options used:
--mesh ../data/fichera.mesh
--order 1
--visualization
Mesh Characteristics:
Number of vertices : 31841
Number of edges : 92256
Number of faces : 89088
Number of elements : 28672
Number of bdr elem : 6144
Euler Number : 1
h_min : 0.0625
h_max : 0.0625
kappa_min : 1
kappa_max : 1

NumOfBdrVertices=6146 NumberOfBdrElements=6144
Number of unknowns: 31841
===============================end of output mfem(suitesparse 4.2.1)

./ex1p -m ../data/periodic-cube.mesh

Hi,

I tried to run "./ex1p -m ../data/periodic-cube.mesh" and visualize with GLVis. It is equivalent to the fully-developed laminar solution in a square duct. However, I only see a perfectly blue cube. It seems that the periodic boundary condition is not applied?

Could you please explain to me why this happened? And, if I just want to simulate a fully-developed laminar solution in a square duct, do I need to modify "periodic-cube.mesh"?

Thank you,
Frank

Initialize Vector data

I've been debugging an application that uses mfem-3.1.
I was tracking down the source of an unused variable using valgrind and found two mfem::Vector functions that allocate, but do not initialize memory:
Vector::Vector (int s)
Vector::SetSize(int s)

Please consider initializing vector data to zero.
It could be as simple as replacing
data = new double[s];
with
data = new doubles;

Is mfem suitable for CFD?

Hi,

Thank you for sharing this nice library.

I have three questions here:

1. Is it possible for mfem to do CFD?
2. I know that spectral-element/hp method commonly uses GLL collection nodes. Is it possible add GLL collection nodes in mfem?
3. Is Neumann boundary condition implemented in mfem?

Thank you,
Frank

New integration rule

Dear all,

I'd like to implement a new integration rule on a segment so that it would be used in square and cube integration rules (the same way it's done right now) and use it along with the existing integration schemes. For example, there is IntegrationRule::UniformRule(), which (as far as I understand) is not used anywhere in the code. Only IntegrationRule::GaussianRule() is. So, basically, my question can be reformulated in the following way: how can I switch to Uniform integration rule working with a high level code - for example, with Example 1. Once I can do that, I'll be able to incorporate a new integration scheme that I need.

Thank you,
Mikhail

Open results in VisIt

Dear all,

I've got some issues trying to open results produced by MFEM in VisIt. For example, after running Example5 and trying to open the file Example5_000000.mfem_root, I've got the following:

VisIt could not read from the file "/u/artemyev/software/mfem-3.0.1/examples/Example5_000000.mfem_root".
The generated error message was:
There was an error opening /u/artemyev/software/mfem-3.0.1/examples/Example5_000000.mfem_root. It may be an invalid file.
VisIt tried using the following file format readers to open the file: Silo
If you know the specific format reader VisIt should use to read this data, you can use Open As... (GUI) or '-o , (CL arg.) and identify that specific reader for VisIt to try. This will possibly give more information on the exact error.

MFEM results are represented not in Silo format, obviously. But I didn't find any JSON support (plugin) in VisIt (I tried version 2.9.2 on Linux Mint 64bit).

So, how should I do that?

Thank you,
Mikhail

Repeating EliminateEssentialBC

Hi,

Sorry for so many questions:) Please take your time to answer my questions.

For my problem, I need to solve Ax=b repeatedly. Fortunately, my A does not change, so A can be assembled only once. However, vector b needs to be updated repeatedly.

My question is how to do "a->EliminateEssentialBC(ess_bdr, x, *b)" repeatedly without destroying ParBilinearForm a? Because constructing ParBilinearForm in every loop is only a waste of time.

The same problem for GetDerivate function: why calculating dshape and inv_jac for every call? why not save these matrix for future use?

Thank you very much for your patience.
Frank

Incompressible NS equation

Hi,

Here is my code solving the incompressible NS equation using MFEM. A mesh for channel flow is also attached here. I added TIMES, ADD, SUBTRACT in gridfunction class.
NS_.zip

This code is NOT fully debugged, but can be a good start if anyone wants to solve NS equation using MFEM.

Hit: For low Reynolds number (say, Re<100), laminar solution can be obtained from this code. However, whenever the convection term is turned on, the code diverges. I have tried both 2/3 and 3/2 dealiasing methods. Neither of them worked. In Nektar++, spectral vanishing viscosity (SVV) or some smooth methods are implemented.

Please let me know if anyone is interested to cooperate. Have a good CFD day.

Thank you,
Frank

Bug in GetDerivative?

Hi,

I tried to test the resolution capability of MFEM by calculating derivatives of sin(k*x). 7 elements are used to divide 2\pi in the x direction and 5th order H1_FECollection is employed. Here are the error when k=1,2,3 and 4. The accuracy is questionable. Is there any way to improve this? I used wrote some spectral element code, which gives much better results than this.

Thank you,
Frank

How to integrate and interpolate?

Hi,

Suppose I have a ParGridFunction u, how to:

1. get the volumetric average
2. interpolate to different order of collection points?

Thank you very much,
Frank

How to copy a SparseMatrix

Dear all,

I can't figure out how I can create a copy of an existing SparseMatrix.
I am solving the following equation:
(M+D)*u_0 = M*(2u_1 - u_2) - <other stuff> + D*u_2,
where M and D are sparse matrices obtained with different bilinear forms corresponding to the same FE space.
I denote the matrix (M+D) as A, and I need it to be different from both M and D, because I need those on the right hand side of the equation.
Doing this:
SparseMatrix A = M; A += D;
copies all internals of the M matrix to the A including pointers (because the SparseMatrix::operator=(const SparseMatrix&) doesn't exist, and the compiler creates its own - a primitive one), so changing the matrix A with A += D; changes the M matrix as well (even if const SparseMatrix &M = m->SpMat();, i.e. the M is supposed to stay const).

Could you please give me a hint how to create a new matrix A=M+D, or is there any other solution for my problem?

Thank you,
Mikhail

How to compute gradient of displacement and other functions in ex2.cpp

I am new to MFEM. In example ex2.cpp I would like to compute gradient of 
displacement field u and then stress sigma(u), strain eps(u), strain energy U 
and elementwise strain energy Ue. The functional form of these quantities as a 
function of u is given by

eps(u) = (grad(u) + grad(u)^T)/2

sigma(u) = 2*mu*eps + lambda*trace(eps) I

U = lmbda/2.0*(trace(eps))^2 + mu*trace(eps^2) 

Ue= \int_(Omega)(U/vol *v dOmega) , where vol is cell volume and v are trial 
functions 

Can you provide an example or illustration for this?
Also, how is it possible to apply a point load at specific location on free end 
of the cantilever beam (ex2.cpp)?

Original issue reported on code.google.com by [email protected] on 17 Mar 2014 at 10:36

GlobalVector() applied to an eigenvector

Dear All,

I'm trying to compute eigenvectors of a problem -Delta u = lambda u with homogeneous Dirichlet boundary conditions following the example 11 (ex11p.cpp). There are only 2 differences:

1. I compute everything with 1 thread, so instead of MPI_COMM_WORLD I use MPI_COMM_SELF
2. I also need to put the eigenvectors all together in a dense matrix, so that one column of the matrix is the eigenvector.

There are no problems with MPI_COMM_SELF, but when I compute the eigenvector, something goes wrong. Here is my addition to the ex11p.cpp right after the eigenvalue problem is solved:

218    DenseMatrix modes(fespace->GetVSize(), nev);
219    for (int i = 0; i < nev; ++i) {
220      const Vector *y = lobpcg->GetEigenvector(i).GlobalVector();
221      Vector x;
222      modes.GetColumnReference(i, x);
223      x = *y;
224    }

This doesn't work for me, however. I suspect there is something in the GlobalVector() therefore I titled the issue that way. But, of course, it's maybe I'm doing something wrong.

Any advice is highly appreciated.

Thanks,
Mikhail

Tag & publish a release of version 3.1.

I'd like to write a package for MFEM in homebrew-science based on the one I wrote in spack. Most package managers (including Homebrew) require a tarball, and won't accept packages that don't have one. This tarball can be created via a tagged release of MFEM version 3.1 on GitHub (see instructions here). Once a tarball is created, I can add an MFEM package to Homebrew, and update my pull request with an MFEM package in spack.

ex3.cpp crash when disable/change the mesh refinement

Dear all,

I am trying to use mfem for my problem. I focus on the electromagnetic 
calculation and use ex3.cpp for my test. For debugging and learning the code of 
mfem, I disabled the refinement from the example:

...
   // 2. Refine the mesh to increase the resolution. In this example we do
   //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
   //    largest number that gives a final mesh with no more than 50,000
   //    elements.
   if(0){   /// disable refinement for test
      int ref_levels =
         (int)floor(log(50000./mesh->GetNE())/log(2.)/mesh->Dimension());
      for (int l = 0; l < ref_levels; l++)
         mesh->UniformRefinement();
   }
   mesh->ReorientTetMesh();
...


The program crashed. Furthermore, the program still crash even I set a smaller 
number refinement number:

      int ref_levels =
         (int)floor(log(30000./mesh->GetNE())/log(2.)/mesh->Dimension());


Could anyone help me to take a look at it?


Thanks,
Tang Laoya 

Original issue reported on code.google.com by [email protected] on 8 Dec 2014 at 2:18

Compilation error with SuiteSparse

Hello,

I am trying to compile MFEM with the flag MFEM_USE_SUITESPARSE ?= YES.
The SUITESPARSE_DIR variable is set to the right directory. I am using the last SuiteSparse sources (v4.4.5), available at http://faculty.cse.tamu.edu/davis/suitesparse.html

The compilation error is (using make serial) :

In file included from linalg.hpp:33:0,
from solvers.cpp:12:
solvers.hpp:23:21: fatal error: umfpack.h: No such file or directory

include < umfpack.h >

MFEM is looking for headers in SuiteSparce_DIR/include but there is no such folder in my SuiteSparse directory. The headers for UMFPACK are in SuiteSparse/UMFPACK, the headers for AMD are in SuiteSparce/AMD etc...

One "working" solution I have found is to add many symlinks in my SuiteSparse directory but this is not really clean.

Am I doing something wrong or the current MFEM makefile was designed to work with other or older versions of the SuiteSparse library ?

Thank you for your help, and thank you for providing MFEM, this is an awsome FEM library !

mfem data directory

How do you generate .mesh file (e.g. beam-tet.mesh) in the data directory of 
the MFEM code? Thanks.

Original issue reported on code.google.com by [email protected] on 13 May 2013 at 9:27

Save results in a specified directory

Hi All,

I'm wondering if there is a way to specify an output directory where I want all the visualization files to be saved? Currently, if I save results in VisIt files, MFEM creates a bunch of directories and corresponding files in a directory from where I call a program. I didn't find any options to change that. That would be especially helpful when one works on a computer cluster where it's typical to have a small home directory where the compiled binary and input files are, and another much larger directory where the big output files are saved.

Thanks,
Mikhail

DerivativeIntegrator does not check the order of integration rule used

Dear MFEM developers,

I have had a problem at the line 949 of bilininteg.cpp in the function
"void DerivativeIntegrator::AssembleElementMatrix2()".
This was due to the order of trial_fe and test_fe was zero for both of them in my particular case,
and consequently, the construction of integration rule at the line 962 was called with order = -1.
I know that calling the DerivativeIntegrator on constant trial&test functions does not make much sense, but this happened during some tests, where I expected the derivative term to simply disappear, i.e. returns zero. Would you mind to add an easy check that treats this issue?

Thanks,
Milan

Asking about MFEM

mind you if you can help me to add MFEM(finite element method library)to integrate with visual studio 2013?

Memory leak in HypreLOBPCG

Hi All,

I've recently come across a memory leak while using HypreLOBPCG, so I'm wondering if you'd like to look at it and check if it's fixable. Here is a minimal example showing the case, which is basically example 11 without comments and extra lines of code for mesh refinement and mesh initialization.

// file lobpcg.cpp
#include "mfem.hpp"
#include <iostream>
using namespace std;
using namespace mfem;

int main(int argc, char *argv[]) {
   int num_procs, myid;
   MPI_Init(&argc, &argv);
   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
   MPI_Comm_rank(MPI_COMM_WORLD, &myid);

   int order = 1;
   int nev = 5;

   Mesh *mesh = new Mesh(10, 10, Element::QUADRILATERAL, 1);
   int dim = mesh->Dimension();

   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
   delete mesh;

   FiniteElementCollection *fec = new H1_FECollection(order, dim);
   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
   HYPRE_Int size = fespace->GlobalTrueVSize();
   if (myid == 0)
      cout << "Number of unknowns: " << size << endl;

   ConstantCoefficient one(1.0);
   Array<int> ess_bdr;
   if (pmesh->bdr_attributes.Size()) {
      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
      ess_bdr = 1;
   }

   ParBilinearForm *a = new ParBilinearForm(fespace);
   a->AddDomainIntegrator(new DiffusionIntegrator(one));
   if (pmesh->bdr_attributes.Size() == 0)
      a->AddDomainIntegrator(new MassIntegrator(one));
   a->Assemble();
   a->EliminateEssentialBCDiag(ess_bdr, 1.0);
   a->Finalize();

   ParBilinearForm *m = new ParBilinearForm(fespace);
   m->AddDomainIntegrator(new MassIntegrator(one));
   m->Assemble();
   // shift the eigenvalue corresponding to eliminated dofs to a large value
   m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
   m->Finalize();

   HypreParMatrix *A = a->ParallelAssemble();
   HypreParMatrix *M = m->ParallelAssemble();

   delete a;
   delete m;

   HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
   amg->SetPrintLevel(0);

   HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
   lobpcg->SetNumModes(nev);
   lobpcg->SetPreconditioner(*amg);
   lobpcg->SetMaxIter(100);
   lobpcg->SetTol(1e-8);
   lobpcg->SetPrecondUsageMode(1);
   lobpcg->SetPrintLevel(1);
   lobpcg->SetMassMatrix(*M);
   lobpcg->SetOperator(*A);

   Array<double> eigenvalues;
   lobpcg->Solve();
   lobpcg->GetEigenvalues(eigenvalues);
   ParGridFunction x(fespace);

   delete lobpcg;
   delete amg;
   delete M;
   delete A;

   delete fespace;
   delete fec;
   delete pmesh;

   MPI_Finalize();
   return 0;
}

When I run it with 1 process and check memory leaks with valgrind there are lots of leaks caused by MPI implementation, which I ignore, and also there is the one from HypreLOBPCG:

==21789== 968 bytes in 1 blocks are definitely lost in loss record 246 of 266
==21789==    at 0x4C29DB4: calloc (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==21789==    by 0x677F80: hypre_CAlloc (in /u/artemyev/projects/mfem/examples/lobpcg)
==21789==    by 0x680AA3: hypre_SeqVectorInitialize (in /u/artemyev/projects/mfem/examples/lobpcg)
==21789==    by 0x68FC11: hypre_ParVectorInitialize (in /u/artemyev/projects/mfem/examples/lobpcg)
==21789==    by 0x4DDDEE: mfem::HypreParVector::HypreParVector(ompi_communicator_t*, int, double*, int*) (hypre.cpp:73)
==21789==    by 0x4E97A9: mfem::HypreLOBPCG::SetOperator(mfem::Operator&) (hypre.cpp:3070)
==21789==    by 0x4C15C1: main (lobpcg.cpp:66)

So I'm wondering if you want to check whether there is indeed a leak.

Thanks,
Mikhail

mat + mat_e of bilinear form

Dear All,

After eliminating bc there is a mat_e sparse matrix such as mat+mat_e=full_matrix of the bilinear form. There are functions like FullMult and FullAddMult that use both mat and mat_e in the matrix-vector multiplication. I am wondering if I can get a full sparse matrix. I didn't find anything existing in the library, so I'm thinking about the following options:

const SparseMatrix &S = billiear_form.FullSpMat();

where

SparseMatrix BilinearForm::FullSpMat() const { return mat + mat_e; } // create a full matrix

Or maybe, something like this:

bilinear_form.RestoreFullMatrix(); // that would involve mat += mat_e
const SparseMatrix &S = bilinear_form.SpMat(); // existing way

Thanks,
Mikhail