ex3.cpp crash when disable/change the mesh refinement

[GoogleCodeExporter]

ex3.cpp crash when disable/change the mesh refinement

[GoogleCodeExporter]

Hi Tang Laoya,

I can't reproduce the crash you are observing. Can you tell us which version of 
the library you are using and what is the command line?

Thanks,
Tzanio

Original comment by tzanio on 8 Dec 2014 at 10:00

• Added labels: Type-Defect, Priority-High

[GoogleCodeExporter]

Dear Prof. Tzanio,

Thank you very much for your kindly reply. I used the latest version 
http://releases.mfem.googlecode.com/hg/mfem-9e90b6710d.tgz
and some modifications are done for running in Windows 7 64bit + VS 2012/VS2013 
+ Intel C++ 15 (till now I still haven't successfully build the project by C++ 
compiler of VS2013).  

When decreasing or removing the refinement, the code crashed at line 810 of 
'bilininteg.cpp':
w *= Q->Eval(Trans, ip);

The command line I used is 
ex3 beam-hex.mesh

I used that file because I noticed that there are multiple domains so I can 
test multiple materials.


BTW: I am trying to study the setting of sources 
   // 4. Set up the linear form b(.) which corresponds to the right-hand side
   //    of the FEM linear system, which in this case is (f,phi_i) where f is
   //    given by the function f_exact and phi_i are the basis functions in the
   //    finite element fespace.
   VectorFunctionCoefficient f(3, f_exact);
   LinearForm *b = new LinearForm(fespace);
   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
   b->Assemble();

and boundary conditions:
   Array<int> ess_bdr(mesh->bdr_attributes.Max());
   ess_bdr = 1;
   a->EliminateEssentialBC(ess_bdr, x, *b);

Is there any more detailed document to describe the code?

In step 6 of this example, the EM diffusion operator is 
curl mu^(-1) curl + sigma I
What's the 'sigma I' stand for?

For my problem, the Maxwell equation (for electric field E) is
curl mu^(-1) curl E - k_0^2 * epsilon = i*omega*mu*J

where k_0 is wave number of vacuum, epsilon is relative permittivity, omega is 
frequency angle and J is current density of source.

How to apply the EM diffusion operator?

Thanks,
Tang Laoya

Original comment by [email protected] on 9 Dec 2014 at 12:28

[GoogleCodeExporter]

Sorry, the equation is
curl mu^(-1) curl E - k_0^2 * epsilon E = i*omega*mu*J 

Thanks,
Tang Laoya

Original comment by [email protected] on 9 Dec 2014 at 12:58

[GoogleCodeExporter]

Dear Prof. Tzanio,

Another question: I noticed that the FEM matrix in ex3.cpp is real number, 
however, for the time-harmonic electromagnetic field, the global FEM matrix 
should be a complex matrix. Is there anything I missed?

Thanks,
Tang Laoya

Original comment by [email protected] on 9 Dec 2014 at 3:28

[GoogleCodeExporter]

Hi Tang Laoya,

I still can't reproduce your bug. Can you verify that you see a crash starting 
from a fresh copy of the 9e90b6710d release?

We do not have support for complex numbers in MFEM, so the problem solved in 
ex3 is the time-domain definite Maxwell equations, where the mass matrix 
coefficient is positive. You can still solve the indefinite problem with 
complex coefficients if you reformulate it as the equivalent 2x2 block system 
for the real and imaginary parts of the electric field and the source current.

Cheers,
Tzanio

Original comment by tzanio on 10 Dec 2014 at 10:10

[GoogleCodeExporter]

Dear Prof. Tzanio,

Thank you very much for your kindly reply. I verified that the crash version is 
a copy of 9e90b6710d, but there are some modification for Windows building, and 
the modification can only work on Intel C++ 15 compiler. I don't know whether 
it is the modification lead to the crash. 

On the other hand, thanks for your suggestion on reformulating the matrix, as 
you suggested, I can assemble the matrix and right hand side two times for real 
part and imaginary part respectively, and then form the complex matrix and then 
use external solver to solve the problem. 

I still have two questions:
1) how to apply the external current source? For my problem, the current source 
is added on a good conductive metal, such as copper pole. When solve the steady 
electric field, I can act the copper pole as equi-potential object so that 
copper pole can be merged into one super-node and the external current source 
are added into that super-node. But how to do for time-harmonic electromagnetic 
field?

2) how to set material when considering anisotropic (the conductivity and 
permittivity are tensor)?

Thanks,
Tang Laoya

Original comment by [email protected] on 11 Dec 2014 at 12:29

[GoogleCodeExporter]

Hi Tang Laoya,

I am not sure I understand your first question, if you have prescribed current 
in portion of the domain you have to  identify that subdomain (e.g. with a 
different attribute in the mesh) and define a piecewise-constant vector field J 
that is given in the metal and zero otherwise.

You can have a matrix coefficient in the VectorFEMassIntegrator (where the 
conductivity is), so anisotropy there shouldn't be a problem. The 
CurlCurlIntegrator does not have a matrix coefficient option, but adding it 
shouldn't be too difficult if you take a careful look at the code in 
fem/bilininteg.cpp.

Hope this helps,
Tzanio

Original comment by tzanio on 14 Dec 2014 at 8:16

[GoogleCodeExporter]

Hi Tang Laoya,

Please let us know if you still see a problem in the new release (3.0). We have 
improved the Windows support (including your changes).

Cheers,
Tzanio

Original comment by tzanio on 28 Jan 2015 at 6:14

[GoogleCodeExporter]

Dear Prof. tzanio,

That's great! Thank you very for your update. I will test the latest
version ASAP.

Original comment by [email protected] on 28 Jan 2015 at 8:03

[GoogleCodeExporter]

Dear Prof. Tzanio,

I test the latest version by ex3.cpp. However, this time the code crashed
at line 72:
   mesh = new Mesh(imesh, 1, 1);

The build environment is Windows 7 64bit + VS2013, or Windows 7 64bit +
VS2013 + Intel C++ 2015.
Could you please help me to take a look at it?

BTW: I noticed that there are much more examples provided. The most
attractive one is ex6.cpp, which solves non-conforming adaptive mesh
refinement for the Laplace problem. I have two questions:
1. could you please also provide an example to solve the problem in ex3.cpp
by non-conforming adaptive mesh refinement?
2. without adaptive mesh refinement, but there are to set of meshes (size
of one set of meshes is meters, but the size of another set of meshes
inside the former is about millimeters), can I solve the problem by
non-conforming meshes?

Original comment by [email protected] on 28 Jan 2015 at 9:40

[GoogleCodeExporter]

Dear Prof. Tzanio,

Further check shows that it crashed at line 1820 of mesh.cpp:
   string mesh_type;

Original comment by [email protected] on 28 Jan 2015 at 1:22

[GoogleCodeExporter]

Hi Tang Laoya,

Regarding your questions about nonconforming meshes. They should work for all 
FE spaces we support, including H(curl) and H(div). Just call 
mesh.GeneralRefinement and a.ConformingAssemble/x.ConformingPronlogate before 
and after the solve, as in ex6.cpp. 

The challenge is of course estimating the error in your solution to determine 
which elements to refine. It may be sufficient to just do predefined 
refinements in a region you care about. In that case I would suggest marking 
that region with an attribute and then have a function to refine elements with 
that attribute (just loop over all elements, collect the right indices and call 
GeneralRefinement). You can then call the function repeatedly to refine the 
region as many times as you like before the computation.

Best regards,
Jakub

Original comment by [email protected] on 28 Jan 2015 at 7:22

[GoogleCodeExporter]

Hi Tang Laoya,

Regarding the crash in line 1820 of mesh.cpp: how did you run ex3, and what was 
the output of the crash/error? We have tested it with VS2013 (on Windows 8), 
and ex3 (e.g. run without arguments) worked fine.

Also, did you modify the mfem-3.0 release at all? Line 1820 in the distribution 
should be "getline(input, mesh_type);"

Thanks,
Tzanio

Original comment by tzanio on 28 Jan 2015 at 8:22

[GoogleCodeExporter]

Dear Prof. Tzanio,

Thank you very much for your kindly reply.
I am sorry I mean line 1818:
   string mesh_type;

Further test I found that the crash problem solved after I added the
preprocessor  'MFEM_THREAD_SAFE' in config.hpp.
Which preprocessors do you suggest to use in Windows system? How do you
modify config.hpp?

Original comment by [email protected] on 29 Jan 2015 at 1:18

[GoogleCodeExporter]

Hi Tang Laoya,

I am glad to hear that you resolved the problem. I am not sure why you needed 
MFEM_THREAD_SAFE, we typically use this only with OpenMP on Linux.

Regarding config.hpp, you just create a copy from config.hpp.in (as described 
at the end of the INSTALL file) and uncomment the #defines you want. All of 
them are optional and some depend on external libraries. I don't have any 
recommendation except that we usually enable MFEM_USE_MEMALLOC.

Hope this helps,
Tzanio

Original comment by tzanio on 29 Jan 2015 at 4:33

• Changed state: Started
• Added labels: OpSys-Windows

[GoogleCodeExporter]

Dear Prof. Tzanio,

Thank you very much for your kindly reply. As you described, I uncommented
#define MFEM_USE_MEMALLOC

in config.hppp, but the code still crash, until I uncommented the line
#define MFEM_THREAD_SAFE

The code works.

But it is strange that you can work on Windows without 'MFEM_THREAD_SAFE'.

Original comment by [email protected] on 29 Jan 2015 at 5:29

[GoogleCodeExporter]

Dear Prof. Tzanio,

I noticed that you have updated MFEM to 3.0.1 and I tested your latest version 
and found that the crash problem was fixed when compiled by Intel C++ compiler 
and Visual C++ compiler. I also noticed that if the 'MFEM_USE_OPENMP' is 
enabled, the 'MFEM_THREAD_SAFE' should also be enabled, else the code would 
crash, you have already commented in config.hpp. Thank you very much.

On the other hand, I noticed the mesh refinement in this example is by 
'UniformRefinement', which means the refinement to every element every time, 
and furthermore, the refinement only support for quadrilateral and hexahedral 
elements, is there any function to do adaptive refinement according to error 
analysis, for example, superconvergent patch recovery and support for 
triangular and tetrahedral?
[1] O. C. Zienkiewicz and J. Z. Zhu, The superconvergent patch recovery and a 
posteriori error estimates. Part 1: The recovery technique [J]. Int. J. Numer. 
Meth. Eng., 1992, 33:1331-1364.
[2] O. C. Zienkiewicz and J. Z. Zhu, The superconvergent patch recovery and a 
posteriori error estimates. Part 2: Error estimates and adaptivity[J]. Int. J. 
Numer. Meth. Eng., 1992, 33:1365-1382.
[3] Elena Zhelezina, Adaptive Finite Element Method for the Numerical 
Simulation of Electric, Magnetic and Acoustic Fields (Ph. D. thesis), 2005.


P.S.: I noticed that MUMPS is a very powerful solver, which can solve both real 
and complex number matrix, do you have any plan to provide interface to it?

Thanks,
Tang Laoya

Original comment by [email protected] on 5 Feb 2015 at 12:37

[GoogleCodeExporter]

HI Tang Laoya,

I am glad to hear that the new version is working better for you. As Jakub 
replied in https://code.google.com/p/mfem/issues/detail?id=11#c12, you can 
(easily) change Example 3 to include non-conforming adaptive refinement 
(similar to Example 6), provided you have an error estimator that determines 
which elements to refine. We currently don't have any sophisticated error 
estimators in MFEM, but may consider adding some in the future. We will also 
consider adding support for MUMPS.

Thanks,
Tzanio

Original comment by tzanio on 6 Feb 2015 at 2:11

[GoogleCodeExporter]

Original comment by tzanio on 6 Feb 2015 at 2:11

• Changed state: Fixed