diff --git a/INSTALL.md b/INSTALL.md index 6e8a65e0..3871e7c2 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -56,6 +56,7 @@ pip install . -C"with-parallel=Yes" -C"with-gslib=Yes" |------|-------------| | `--with-parallel` | Install both serial and parallel versions of `MFEM` and the wrapper
(note: this option turns on building `metis` and `hypre`) | | `--mfem-branch=` | Download/install MFEM using a specific reference (`git` `branch`, `hash`, or `tag`) | +| `--mfem-miniapps` | Install MFEM with MFEM C++ miniapps | | `--user` | Install in user's site-package | ## Advanced options diff --git a/_build_system/build_config.py b/_build_system/build_config.py index a82bf1f8..71d07c10 100644 --- a/_build_system/build_config.py +++ b/_build_system/build_config.py @@ -36,6 +36,7 @@ def print_config(): print(" when needed, the dependency (mfem/hypre/metis) will be installed under " + bglb.ext_prefix) print(" build mfem : " + ("Yes" if bglb.build_mfem else "No")) + print(" build miniapps: " + ("Yes" if bglb.mfem_miniapps else "No")) print(" build metis : " + ("Yes" if bglb.build_metis else "No")) print(" build hypre : " + ("Yes" if bglb.build_hypre else "No")) print(" build libceed : " + ("Yes" if bglb.build_libceed else "No")) @@ -89,7 +90,7 @@ def initialize_cmd_options(command_obj): command_obj.mfem_source = bglb.mfem_source command_obj.mfem_branch = '' command_obj.mfem_debug = False - command_obj.mfem_build_miniapps = False + command_obj.mfem_miniapps = True command_obj.metis_prefix = '' command_obj.hypre_prefix = '' @@ -146,7 +147,7 @@ def initialize_cmd_options(command_obj): ('mfem-source=', None, 'Specify mfem source location' + 'MFEM source directory. Required to run-swig '), ('mfem-debug', None, 'Build MFME with MFEM_DEBUG enabled'), - ('mfem-build-miniapps', None, 'build MFME Miniapps'), + ('mfem-miniapps', None, 'build MFME Miniapps'), ('hypre-prefix=', None, 'Specify locaiton of hypre' + 'libHYPRE.so must exits under /lib'), ('metis-prefix=', None, 'Specify locaiton of metis' + @@ -221,15 +222,22 @@ def process_cmd_options(command_obj, cfs): assert False, str(command_obj) + " does not have " + attr setattr(command_obj, attr, value) else: - value = cfs.pop(param, "No") if not hasattr(command_obj, attr): assert False, str(command_obj) + " does not have " + attr + if getattr(command_obj, attr): + value = cfs.pop(param, "Yes") + else: + value = cfs.pop(param, "No") + if value.upper() in ("YES", "TRUE", "1"): setattr(command_obj, attr, True) else: setattr(command_obj, attr, False) + if len(cfs) != 0: + assert False, "unknonw input is given " + str(cfs) + def process_setup_options(command_obj, args): for item in args: @@ -298,7 +306,7 @@ def configure_install(self): bglb.run_swig_parallel = bool(self.with_parallel) bglb.mfem_debug = bool(self.mfem_debug) - bglb.mfem_build_miniapps = bool(self.mfem_build_miniapps) + bglb.mfem_miniapps = bool(self.mfem_miniapps) if bglb.build_serial: bglb.build_serial = (not bglb.swig_only and not bglb.ext_only) @@ -439,7 +447,6 @@ def configure_install(self): bglb.build_parallel = False bglb.keep_temp = True - if bglb.libceed_only: bglb.clean_swig = False bglb.run_swig = False @@ -453,7 +460,6 @@ def configure_install(self): bglb.build_libceed = True bglb.keep_temp = True - if bglb.gslib_only: bglb.clean_swig = False bglb.run_swig = False @@ -466,7 +472,7 @@ def configure_install(self): bglb.build_gslib = True bglb.keep_temp = True - bglb.is_configured = True + configure_build = configure_install diff --git a/_build_system/build_mfem.py b/_build_system/build_mfem.py index 2e55853e..67781d1c 100644 --- a/_build_system/build_mfem.py +++ b/_build_system/build_mfem.py @@ -43,7 +43,7 @@ def add_rpath(p, dest): cmake_opts = {'DBUILD_SHARED_LIBS': '1', 'DMFEM_ENABLE_EXAMPLES': '1', - 'DMFEM_ENABLE_MINIAPPS': '0', + 'DMFEM_ENABLE_MINIAPPS': '1', 'DCMAKE_SHARED_LINKER_FLAGS': ldflags, 'DMFEM_USE_ZLIB': '1', 'DCMAKE_CXX_FLAGS': bglb.cxxstd_flag, @@ -56,8 +56,8 @@ def add_rpath(p, dest): if bglb.mfem_debug: cmake_opts['DMFEM_DEBUG'] = 'YES' - if bglb.mfem_build_miniapps: - cmake_opts['DMFEM_ENABLE_MINIAPPS'] = '1' + if not bglb.mfem_miniapps: + cmake_opts['DMFEM_ENABLE_MINIAPPS'] = '0' if bglb.verbose: cmake_opts['DCMAKE_VERBOSE_MAKEFILE'] = '1' diff --git a/_build_system/build_pymfem.py b/_build_system/build_pymfem.py index fb572e86..585ca41b 100644 --- a/_build_system/build_pymfem.py +++ b/_build_system/build_pymfem.py @@ -71,6 +71,7 @@ def write_setup_local(): 'gslibpinc': os.path.join(bglb.gslibp_prefix, 'include'), 'cxxstdflag': bglb.cxxstd_flag, 'build_mfem': '1' if bglb.build_mfem else '0', + 'build_miniapps': '1' if bglb.mfem_miniapps else '0', 'bdist_wheel_dir': bglb.bdist_wheel_dir, } diff --git a/docs/changelog.txt b/docs/changelog.txt index ba4501b8..4374a5f2 100644 --- a/docs/changelog.txt +++ b/docs/changelog.txt @@ -2,6 +2,14 @@ 2025 12 * PyMFEM4.9 - ordering.i is added and made correponding adjustment in *.i files. + - update ex22, ex22p to save ComplexGridFunction as is + - fixed ParComplexGridFunction::Save wrapping + - fixed GetDataArray of Vector/DenseMatrix/DenseTensor/SparseMatrix so that + destructore is not pre-maturely called when return value is not stored to + a variable. + - particleset.hpp, particlevector.hpp wrapper + - miniapp/dpg/utils/*.hpp is wrapped. + - Python version of pmaxwell.py in miniapp/dpg directory 2025 09 * PyMFEM 4.8 diff --git a/mfem/_par/blockoperator.i b/mfem/_par/blockoperator.i index 9ec894c9..165c416b 100644 --- a/mfem/_par/blockoperator.i +++ b/mfem/_par/blockoperator.i @@ -76,9 +76,15 @@ if len(args) == 2: mfem::SparseMatrix *Opr2SparseMat(mfem::Operator *op) { return dynamic_cast(op); } + mfem::SparseMatrix *Opr2SparseMatrix(mfem::Operator *op) { + return dynamic_cast(op); + } mfem::HypreParMatrix *Opr2HypreParMat(mfem::Operator *op) { return dynamic_cast(op); } + mfem::HypreParMatrix *Opr2HypreParMatrix(mfem::Operator *op) { + return dynamic_cast(op); + } %} %include "linalg/blockoperator.hpp" diff --git a/mfem/_par/blockstaticcond.i b/mfem/_par/blockstaticcond.i new file mode 100644 index 00000000..07f34c32 --- /dev/null +++ b/mfem/_par/blockstaticcond.i @@ -0,0 +1,46 @@ +%module(package="mfem._par") blockstaticcond +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "miniapps/dpg/util/blockstaticcond.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_BLOCKSTATICCOND + +%init %{ +import_array(); +%} + +%inline %{ +#include "miniapps/dpg/util/blockstaticcond.cpp" +%} + + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "array.i" +%import "vector.i" +%import "densemat.i" +%import "operators.i" +%import "blockmatrix.i" +%import "blockoperator.i" +%import "fespace.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +OSTREAM_TYPEMAP(std::ostream&) + +%ignore mfem::BlockStaticCondensation::ConvertListToReducedTrueDofs; +%include "miniapps/dpg/util/blockstaticcond.hpp" + +#endif diff --git a/mfem/_par/complex_densemat.i b/mfem/_par/complex_densemat.i new file mode 100644 index 00000000..479ca99f --- /dev/null +++ b/mfem/_par/complex_densemat.i @@ -0,0 +1,30 @@ +%module(package="mfem._par") complex_densemat +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + +%init %{ +import_array(); +%} + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "array.i" +%import "vector.i" +%import "densemat.i" +%import "complex_operator.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +%ignore mfem::ComplexLUFactors::Mult(int m, int n, std::complex *X) const; +%include "linalg/complex_densemat.hpp" diff --git a/mfem/_par/complexstaticcond.i b/mfem/_par/complexstaticcond.i new file mode 100644 index 00000000..641c6f63 --- /dev/null +++ b/mfem/_par/complexstaticcond.i @@ -0,0 +1,48 @@ +%module(package="mfem._par") complexstaticcond +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "miniapps/dpg/util/complexstaticcond.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_COMPLEXSTATICCOND + +%init %{ +import_array(); +%} + +%inline %{ +#include "miniapps/dpg/util/complexstaticcond.cpp" +%} + + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "array.i" +%import "vector.i" +%import "gridfunc.i" +%import "mesh.i" +%import "solvers.i" +%import "operators.i" +%import "blockmatrix.i" +%import "blockoperator.i" +%import "complex_densemat.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +OSTREAM_TYPEMAP(std::ostream&) + +%ignore mfem::ComplexBlockStaticCondensation::ConvertListToReducedTrueDofs; +%include "miniapps/dpg/util/complexstaticcond.hpp" + +#endif diff --git a/mfem/_par/complexweakform.i b/mfem/_par/complexweakform.i new file mode 100644 index 00000000..d133b7b3 --- /dev/null +++ b/mfem/_par/complexweakform.i @@ -0,0 +1,65 @@ +%module(package="mfem._par") complexweakform +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "miniapps/dpg/util/complexweakform.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pylininteg.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_COMPLEXWEAKFORM + +%init %{ +import_array(); +%} + +%inline %{ +#include "miniapps/dpg/util/complexweakform.cpp" +%} + + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "coefficient.i" +%import "gridfunc.i" +%import "mesh.i" +%import "solvers.i" +%import "operators.i" +%import "blockmatrix.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +OSTREAM_TYPEMAP(std::ostream&) + +%import "../common/object_array_typemap.i" +LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array&, + mfem::FiniteElementSpace*) +LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array&, + mfem::FiniteElementCollection*) + + +%pythonprepend mfem::ComplexDPGWeakForm::ComplexDPGWeakForm %{ + if len(args) > 0: + fes_, fecol_ = args + self._fes = fes_ + self._fecol = fecol_ +%} + +%pythonprepend mfem::ComplexDPGWeakForm::SetSpaces %{ + self._fes = fes_ + self._fecol = fecol_ +%} + + +%include "miniapps/dpg/util/complexweakform.hpp" + +#endif diff --git a/mfem/_par/densemat.i b/mfem/_par/densemat.i index 6042d228..aacf9f84 100644 --- a/mfem/_par/densemat.i +++ b/mfem/_par/densemat.i @@ -124,6 +124,23 @@ def __getitem__(self, *args): return _densemat.DenseTensor___getitem__(self, check) %} +%typemap(out) PyObject* mfem::DenseMatrix::GetDataArray() { + // assign self to base object + Py_INCREF($self); + PyObject * base = PyArray_BASE((PyArrayObject *) $1); + PyArray_SetBaseObject((PyArrayObject *) base , $self); + + $result = $1; +} +%typemap(out) PyObject* mfem::DenseTensor::GetDataArray() { + // assign self to base object + Py_INCREF($self); + PyObject * base = PyArray_BASE((PyArrayObject *) $1); + PyArray_SetBaseObject((PyArrayObject *) base, $self); + + $result = $1; +} + %include "linalg/densemat.hpp" %extend mfem::DenseMatrix { diff --git a/mfem/_par/dpg.i b/mfem/_par/dpg.i new file mode 100644 index 00000000..c7d1c9a2 --- /dev/null +++ b/mfem/_par/dpg.i @@ -0,0 +1,51 @@ +%module(package="mfem._par") dpg + +%feature("autodoc", "1"); + +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pylininteg.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +using namespace mfem; +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_WEAKFORM + + +%init %{ +import_array(); +%} + +%include "exception.i" +%include "../common/exception.i" + +%inline %{ +#include "miniapps/dpg/util/pml.cpp" +#include "miniapps/dpg/util/blockstaticcond.cpp" +#include "miniapps/dpg/util/complexstaticcond.cpp" +#include "miniapps/dpg/util/weakform.cpp" +#include "miniapps/dpg/util/complexweakform.cpp" +#include "miniapps/dpg/util/pweakform.cpp" +#include "miniapps/dpg/util/pcomplexweakform.cpp" +%} + +%include "../common/dpg_common.i" + +%include "miniapps/dpg/util/pml.hpp" +%include "miniapps/dpg/util/blockstaticcond.hpp" +%include "miniapps/dpg/util/complexstaticcond.hpp" +%include "miniapps/dpg/util/weakform.hpp" +%include "miniapps/dpg/util/complexweakform.hpp" +%include "miniapps/dpg/util/pweakform.hpp" +%include "miniapps/dpg/util/pcomplexweakform.hpp" + + +#endif diff --git a/mfem/_par/fe_coll.i b/mfem/_par/fe_coll.i index b9c1099a..dbfc27bb 100644 --- a/mfem/_par/fe_coll.i +++ b/mfem/_par/fe_coll.i @@ -33,6 +33,13 @@ import_array1(-1); %import "eltrans.i" %import "lininteg.i" +/* define FiniteElementCollectionPtrArray */ +%import "../common/array_listtuple_typemap.i" +ARRAY_LISTTUPLE_INPUT_SWIGOBJ(mfem::FiniteElementCollection *, 1) +%import "../common/array_instantiation_macro.i" +IGNORE_ARRAY_METHODS(mfem::FiniteElementCollection *) +INSTANTIATE_ARRAY0(FiniteElementCollection *, FiniteElementCollection, 1) + %include "fem/fe_coll.hpp" %pythoncode %{ DG_FECollection = L2_FECollection diff --git a/mfem/_par/handle.i b/mfem/_par/handle.i index 84229488..2d8af480 100644 --- a/mfem/_par/handle.i +++ b/mfem/_par/handle.i @@ -32,6 +32,7 @@ import_array1(-1); %import "../common/exception.i" %import "operators.i" +%import "complex_operator.i" #ifdef MFEM_USE_MPI %import "hypre.i" @@ -60,6 +61,13 @@ GET_RENAME(HypreParMatrix) RESET_RENAME(HypreParMatrix) CONVERT_FROM_RENAME(HypreParMatrix) +//CONSTRUCTOR_WRAP(ComplexOperator) +AS_RENAME(ComplexOperator) +IS_RENAME(ComplexOperator) +GET_RENAME(ComplexOperator) +RESET_RENAME(ComplexOperator) +CONVERT_FROM_RENAME(ComplexOperator) + //CONSTRUCTOR_RENAME(ComplexHypreParMatrix) AS_RENAME(ComplexHypreParMatrix) IS_RENAME(ComplexHypreParMatrix) @@ -97,6 +105,13 @@ GET_WRAP(HypreParMatrix) RESET_WRAP(HypreParMatrix) CONVERT_FROM_WRAP(HypreParMatrix) +//CONSTRUCTOR_WRAP(ComplexOperator) +AS_WRAP(ComplexOperator) +IS_WRAP(ComplexOperator) +GET_WRAP(ComplexOperator) +RESET_WRAP(ComplexOperator) +CONVERT_FROM_WRAP(ComplexOperator) + //CONSTRUCTOR_WRAP(ComplexHypreParMatrix) AS_WRAP(ComplexHypreParMatrix) IS_WRAP(ComplexHypreParMatrix) @@ -113,3 +128,5 @@ RESET_WRAP(PetscParMatrix) CONVERT_FROM_WRAP(PetscParMatrix) #endif + + diff --git a/mfem/_par/hypre.i b/mfem/_par/hypre.i index f4841540..8915b484 100644 --- a/mfem/_par/hypre.i +++ b/mfem/_par/hypre.i @@ -104,11 +104,11 @@ bool is_HYPRE_USING_CUDA(){ is 4, cols is set to be rows. */ -%typemap(in) (int *I, - HYPRE_BigInt *J, - mfem::real_t *data, - HYPRE_BigInt *rows, - HYPRE_BigInt *cols) +%typemap(in) (const int *I, + const HYPRE_BigInt *J, + const mfem::real_t *data, + const HYPRE_BigInt *rows, + const HYPRE_BigInt *cols) (PyArrayObject *tmp_arr1_ = NULL, PyArrayObject *tmp_arr2_ = NULL, PyArrayObject *tmp_arr3_ = NULL, @@ -133,8 +133,11 @@ bool is_HYPRE_USING_CUDA(){ $5 = (HYPRE_BigInt *) PyArray_DATA(tmp_arr5_); } } -%typemap(freearg) (int *I, HYPRE_BigInt *J, - mfem::real_t *data, HYPRE_BigInt *rows, HYPRE_BigInt *cols){ +%typemap(freearg) (const int *I, + const HYPRE_BigInt *J, + const mfem::real_t *data, + const HYPRE_BigInt *rows, + const HYPRE_BigInt *cols){ Py_XDECREF(tmp_arr1_$argnum); Py_XDECREF(tmp_arr2_$argnum); Py_XDECREF(tmp_arr3_$argnum); @@ -144,9 +147,11 @@ bool is_HYPRE_USING_CUDA(){ } } -%typemap(typecheck ) (int *I, HYPRE_BigInt *J, - mfem::real_t *data, HYPRE_BigInt *rows, - HYPRE_BigInt *cols){ +%typemap(typecheck ) (const int *I, + const HYPRE_BigInt *J, + const mfem::real_t *data, + const HYPRE_BigInt *rows, + const HYPRE_BigInt *cols){ /* check if list of 5 numpy array or not */ if (!PyList_Check($input)) $1 = 0; else { diff --git a/mfem/_par/ordering.i b/mfem/_par/ordering.i index cd86978f..0267785b 100644 --- a/mfem/_par/ordering.i +++ b/mfem/_par/ordering.i @@ -13,10 +13,12 @@ import_array1(-1); %} %include "exception.i" -%import "element.i" + %include "../common/typemap_macros.i" %include "../common/exception.i" +%import "array.i" +%import "vector.i" %include "linalg/ordering.hpp" diff --git a/mfem/_par/particleset.i b/mfem/_par/particleset.i new file mode 100644 index 00000000..57fc0a46 --- /dev/null +++ b/mfem/_par/particleset.i @@ -0,0 +1,30 @@ +// +// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved. +// +%module(package="mfem._par") particleset + +%feature("autodoc", "1"); + +%{ +#include "mfem.hpp" +#include "../common/io_stream.hpp" +#include "numpy/arrayobject.h" +#include "../common/pyoperator.hpp" +#include "../common/pyintrules.hpp" +%} + +%init %{ +import_array1(-1); +%} + +%include "exception.i" +%include "../common/typemap_macros.i" +%include "../common/exception.i" +%include "../common/kernel_dispatch.i" + +%import "array.i" +%import "ordering.i" +%import "particlevector.i" + +%include "fem/particleset.hpp" + diff --git a/mfem/_par/particlevector.i b/mfem/_par/particlevector.i new file mode 100644 index 00000000..4d47770a --- /dev/null +++ b/mfem/_par/particlevector.i @@ -0,0 +1,28 @@ +// +// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved. +// +%module(package="mfem._par") particlevector + +%feature("autodoc", "1"); + +%{ +#include "mfem.hpp" +#include "../common/io_stream.hpp" +#include "numpy/arrayobject.h" +#include "../common/pyoperator.hpp" +#include "../common/pyintrules.hpp" +%} + +%init %{ +import_array1(-1); +%} + +%include "exception.i" +%include "../common/typemap_macros.i" +%include "../common/exception.i" + +%import "vector.i" +%import "ordering.i" + +%include "linalg/particlevector.hpp" + diff --git a/mfem/_par/pcomplexweakform.i b/mfem/_par/pcomplexweakform.i new file mode 100644 index 00000000..d6ef0a88 --- /dev/null +++ b/mfem/_par/pcomplexweakform.i @@ -0,0 +1,66 @@ +%module(package="mfem._par") pcomplexweakform +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "miniapps/dpg/util/pcomplexweakform.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_PCOMPLEXWEAKFORM + +%init %{ +import_array(); +%} + +%inline %{ +#include "miniapps/dpg/util/pcomplexweakform.cpp" +%} + + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "array.i" +%import "vector.i" +%import "gridfunc.i" +%import "mesh.i" +%import "solvers.i" +%import "operators.i" +%import "blockmatrix.i" +%import "complex_densemat.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +OSTREAM_TYPEMAP(std::ostream&) + +%import "../common/object_array_typemap.i" +LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array&, + mfem::ParFiniteElementSpace*) +LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array&, + mfem::FiniteElementCollection*) + + +%pythonprepend mfem::ParComplexDPGWeakForm::ParComplexDPGWeakForm %{ + if len(args) > 0: + trial_pfes_, fecol_ = args + self._fes = trial_pfes_ + self._fecol = fecol_ + +%} + +%pythonprepend mfem::ParComplexDPGWeakForm::SetParSpaces %{ + self._fes = trial_pfes_ + self._fecol = fecol_ +%} + +%include "miniapps/dpg/util/pcomplexweakform.hpp" + +#endif diff --git a/mfem/_par/pfespace.i b/mfem/_par/pfespace.i index 569c45ed..fc6dd3dd 100644 --- a/mfem/_par/pfespace.i +++ b/mfem/_par/pfespace.i @@ -79,6 +79,13 @@ def GetFaceNbrElementVDofs(self, i): return vdofs.ToList() %} +/* define FiniteElementSpaceArray */ +%import "../common/array_listtuple_typemap.i" +ARRAY_LISTTUPLE_INPUT_SWIGOBJ(mfem::ParFiniteElementSpace *, 1) +%import "../common/array_instantiation_macro.i" +IGNORE_ARRAY_METHODS(mfem::ParFiniteElementSpace *) +INSTANTIATE_ARRAY0(ParFiniteElementSpace *, ParFiniteElementSpace, 1) + %include "fem/pfespace.hpp" %extend mfem::ParFiniteElementSpace{ diff --git a/mfem/_par/pml.i b/mfem/_par/pml.i new file mode 100644 index 00000000..ba004a04 --- /dev/null +++ b/mfem/_par/pml.i @@ -0,0 +1,44 @@ +%module(package="mfem._par") pml +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "miniapps/dpg/util/pcomplexweakform.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_PML + +%init %{ +import_array(); +%} + +%inline %{ +#include "miniapps/dpg/util/pml.cpp" +%} + + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "array.i" +%import "vector.i" +%import "gridfunc.i" +%import "mesh.i" +%import "solvers.i" +%import "operators.i" +%import "blockmatrix.i" +%import "complex_densemat.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +%include "miniapps/dpg/util/pml.hpp" + +#endif diff --git a/mfem/_par/pweakform.i b/mfem/_par/pweakform.i new file mode 100644 index 00000000..124af89d --- /dev/null +++ b/mfem/_par/pweakform.i @@ -0,0 +1,66 @@ +%module(package="mfem._par") pweakform +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "miniapps/dpg/util/pweakform.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_PWEAKFORM + +%init %{ +import_array(); +%} + +%inline %{ +#include "miniapps/dpg/util/pweakform.cpp" +%} + + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "array.i" +%import "vector.i" +%import "gridfunc.i" +%import "mesh.i" +%import "solvers.i" +%import "operators.i" +%import "blockmatrix.i" +%import "complex_densemat.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +OSTREAM_TYPEMAP(std::ostream&) + +%import "../common/object_array_typemap.i" +LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array&, + mfem::ParFiniteElementSpace*) +LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array&, + mfem::FiniteElementCollection*) + + +%pythonprepend mfem::ParDPGWeakForm::ParDPGWeakForm %{ + if len(args) > 0: + trial_fes_, fecol_ = args + self._fes = trial_pfes_ + self._fecol = fecol_ + +%} + +%pythonprepend mfem::ParDPGWeakForm::SetParSpaces %{ + self._fes = trial_pfes_ + self._fecol = fecol_ +%} + +%include "miniapps/dpg/util/pweakform.hpp" + +#endif diff --git a/mfem/_par/setup.py b/mfem/_par/setup.py index d373cc61..5e47104f 100644 --- a/mfem/_par/setup.py +++ b/mfem/_par/setup.py @@ -40,7 +40,7 @@ def get_extensions(): mfemptpl, numpyinc, mpi4pyinc, mpiinc, hypreinc, metisinc, hyprelib, metis5lib, cc_par, cxx_par, cc_ser, cxx_ser, - cxxstdflag, mfem_outside, + cxxstdflag, mfem_outside, build_miniapps, add_pumi, add_cuda, add_libceed, add_strumpack, add_suitesparse, add_gslibp, bdist_wheel_dir) @@ -67,9 +67,12 @@ def get_extensions(): add_gslibp = '' cxxstdflag = '-std=c++17' mfem_outside = '0' + build_miniapps = '0' mpiinc = '' libraries = ['mfem',] + if build_miniapps != '0': + libraries.append("mfem-common") # remove current directory from path print("__file__", os.path.abspath(__file__)) @@ -126,7 +129,10 @@ def get_extensions(): "submesh", "transfermap", "staticcond", "sidredatacollection", "psubmesh", "ptransfermap", "enzyme", "attribute_sets", "arrays_by_name", - "hyperbolic", "bounds", "integrator", "ordering"] + "hyperbolic", "complex_densemat", + "bounds", "integrator", "ordering", + "dpg", "particleset", "particlevector"] + if mpiinc != '': include_dirs.append(mpiinc) diff --git a/mfem/_par/sparsemat.i b/mfem/_par/sparsemat.i index 3e80ad17..893fea03 100644 --- a/mfem/_par/sparsemat.i +++ b/mfem/_par/sparsemat.i @@ -147,6 +147,28 @@ if len(args) == 1 and isinstance(args[0], csr_matrix): } } +%typemap(out) PyObject* mfem::SparseMatrix::GetDataArray() { + // assign self to base object + Py_INCREF($self); + PyArray_SetBaseObject((PyArrayObject *) $1, $self); + + $result = $1; +} +%typemap(out) PyObject* mfem::SparseMatrix::GetJArray() { + // assign self to base object + Py_INCREF($self); + PyArray_SetBaseObject((PyArrayObject *) $1, $self); + + $result = $1; +} +%typemap(out) PyObject* mfem::SparseMatrix::GetIArray() { + // assign self to base object + Py_INCREF($self); + PyArray_SetBaseObject((PyArrayObject *) $1, $self); + + $result = $1; +} + //%rename(add_sparse) mfem::Add; %include "linalg/sparsemat.hpp" diff --git a/mfem/_par/vector.i b/mfem/_par/vector.i index 8e2d409d..2bff502c 100644 --- a/mfem/_par/vector.i +++ b/mfem/_par/vector.i @@ -168,6 +168,14 @@ XXXPTR_SIZE_IN(mfem::Vector **data_, int asize, mfem::Vector *) IGNORE_ARRAY_METHODS(mfem::Vector *) INSTANTIATE_ARRAY0(Vector *, Vector, 1) +%typemap(out) PyObject* mfem::Vector::GetDataArray() { + // assign self to base object + Py_INCREF($self); + PyArray_SetBaseObject((PyArrayObject *) $1, $self); + + $result = $1; +} + %include "linalg/vector.hpp" %extend mfem::Vector { diff --git a/mfem/_par/weakform.i b/mfem/_par/weakform.i new file mode 100644 index 00000000..251ff80c --- /dev/null +++ b/mfem/_par/weakform.i @@ -0,0 +1,67 @@ +%module(package="mfem._par") weakform +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "miniapps/dpg/util/weakform.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pylininteg.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_WEAKFORM + + +%init %{ +import_array(); +%} + +%inline %{ +#include "miniapps/dpg/util/weakform.cpp" +%} + + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "array.i" +%import "vector.i" +%import "blockvector.i" +%import "bilininteg.i" +%import "lininteg.i" +%import "fespace.i" +%import "blockmatrix.i" +%import "operators.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +OSTREAM_TYPEMAP(std::ostream&) + +%import "../common/object_array_typemap.i" +LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array&, + mfem::FiniteElementSpace*) +LIST_TO_MFEMOBJ_ARRAY_IN(const mfem::Array&, + mfem::FiniteElementCollection*) + + +%pythonprepend mfem::DPGWeakForm::DPGWeakForm %{ + if len(args) > 0: + fes_, fecol_ = args + self._fes = fes_ + self._fecol = fecol_ +%} + +%pythonprepend mfem::DPGWeakForm::SetSpaces %{ + self._fes = fes_ + self._fecol = fecol_ +%} + +%include "miniapps/dpg/util/weakform.hpp" + +#endif diff --git a/mfem/_ser/blockoperator.i b/mfem/_ser/blockoperator.i index 7b98164f..e1a7b20b 100644 --- a/mfem/_ser/blockoperator.i +++ b/mfem/_ser/blockoperator.i @@ -77,6 +77,9 @@ if len(args) == 2: mfem::SparseMatrix *Opr2SparseMat(mfem::Operator *op) { return dynamic_cast(op); } + mfem::SparseMatrix *Opr2SparseMatrix(mfem::Operator *op) { + return dynamic_cast(op); + } %} %include "linalg/blockoperator.hpp" diff --git a/mfem/_ser/blockstaticcond.i b/mfem/_ser/blockstaticcond.i new file mode 100644 index 00000000..71655353 --- /dev/null +++ b/mfem/_ser/blockstaticcond.i @@ -0,0 +1,46 @@ +%module(package="mfem._ser") blockstaticcond +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "miniapps/dpg/util/blockstaticcond.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_BLOCKSTATICCOND + +%init %{ +import_array(); +%} + +%inline %{ +#include "miniapps/dpg/util/blockstaticcond.cpp" +%} + + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "array.i" +%import "vector.i" +%import "densemat.i" +%import "operators.i" +%import "blockmatrix.i" +%import "blockoperator.i" +%import "fespace.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +OSTREAM_TYPEMAP(std::ostream&) + +%ignore mfem::BlockStaticCondensation::ConvertListToReducedTrueDofs; +%include "miniapps/dpg/util/blockstaticcond.hpp" + +#endif diff --git a/mfem/_ser/complex_densemat.i b/mfem/_ser/complex_densemat.i index df906e60..0b7c816b 100644 --- a/mfem/_ser/complex_densemat.i +++ b/mfem/_ser/complex_densemat.i @@ -29,4 +29,5 @@ import_array1(-1); %import "../common/exception.i" %import "../common/io_stream_typemap.i" +%ignore mfem::ComplexLUFactors::Mult(int m, int n, std::complex *X) const; %include "linalg/complex_densemat.hpp" diff --git a/mfem/_ser/complexstaticcond.i b/mfem/_ser/complexstaticcond.i index 6fbf47d0..4294c6cc 100644 --- a/mfem/_ser/complexstaticcond.i +++ b/mfem/_ser/complexstaticcond.i @@ -15,14 +15,13 @@ #include "../common/io_stream.hpp" %} +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_COMPLEXSTATICCOND + %init %{ import_array1(-1); %} -%include "../common/existing_mfem_headers.i" - -#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_COMPLEXSTATICCOND - %inline %{ #include "miniapps/dpg/util/complexstaticcond.cpp" %} @@ -39,13 +38,14 @@ import_array1(-1); %import "solvers.i" %import "operators.i" %import "blockmatrix.i" +%import "blockoperator.i" %import "complex_densemat.i" %import "../common/exception.i" %import "../common/io_stream_typemap.i" OSTREAM_TYPEMAP(std::ostream&) - +%ignore mfem::ComplexBlockStaticCondensation::ConvertListToReducedTrueDofs; %include "miniapps/dpg/util/complexstaticcond.hpp" #endif diff --git a/mfem/_ser/complexweakform.i b/mfem/_ser/complexweakform.i index 1e831a44..4da4e4fd 100644 --- a/mfem/_ser/complexweakform.i +++ b/mfem/_ser/complexweakform.i @@ -16,14 +16,13 @@ #include "../common/io_stream.hpp" %} +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_COMPLEXWEAKFORM + %init %{ import_array1(-1); %} -%include "../common/existing_mfem_headers.i" - -#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_COMPLEXWEAKFORM - %inline %{ #include "miniapps/dpg/util/complexweakform.cpp" %} @@ -34,16 +33,38 @@ import_array1(-1); %import "../common/exception.i" %import "coefficient.i" +%import "fe_coll.i" %import "gridfunc.i" %import "mesh.i" %import "solvers.i" %import "operators.i" +%import "fespace.i" %import "blockmatrix.i" %import "../common/exception.i" %import "../common/io_stream_typemap.i" OSTREAM_TYPEMAP(std::ostream&) +%include "../common/typemap_macros.i" +LIST_TO_MFEMOBJ_ARRAY_IN(mfem::Array&, + mfem::FiniteElementSpace*) +LIST_TO_MFEMOBJ_ARRAY_IN(mfem::Array&, + mfem::FiniteElementCollection*) + + +%pythonprepend mfem::ComplexDPGWeakForm::ComplexDPGWeakForm %{ + if len(args) > 0: + fes_, fecol_ = args + self._fes = fes_ + self._fecol = fecol_ +%} + +%pythonprepend mfem::ComplexDPGWeakForm::SetSpaces %{ + self._fes = fes_ + self._fecol = fecol_ +%} + + %include "miniapps/dpg/util/complexweakform.hpp" #endif diff --git a/mfem/_ser/densemat.i b/mfem/_ser/densemat.i index e176caec..dc80969c 100644 --- a/mfem/_ser/densemat.i +++ b/mfem/_ser/densemat.i @@ -124,6 +124,23 @@ def __getitem__(self, *args): return _densemat.DenseTensor___getitem__(self, check) %} +%typemap(out) PyObject* mfem::DenseMatrix::GetDataArray() { + // assign self to base object + Py_INCREF($self); + PyObject * base = PyArray_BASE((PyArrayObject *) $1); + PyArray_SetBaseObject((PyArrayObject *) base, $self); + + $result = $1; +} +%typemap(out) PyObject* mfem::DenseTensor::GetDataArray() { + // assign self to base object + Py_INCREF($self); + PyObject * base = PyArray_BASE((PyArrayObject *) $1); + PyArray_SetBaseObject((PyArrayObject *) base, $self); + + $result = $1; +} + %include "linalg/densemat.hpp" %extend mfem::DenseMatrix { diff --git a/mfem/_ser/dpg.i b/mfem/_ser/dpg.i new file mode 100644 index 00000000..fe2dcdde --- /dev/null +++ b/mfem/_ser/dpg.i @@ -0,0 +1,43 @@ +%module(package="mfem._ser") dpg +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pylininteg.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +using namespace mfem; +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_WEAKFORM + + +%init %{ +import_array(); +%} + +%include "exception.i" +%include "../common/exception.i" + +%inline %{ +#include "miniapps/dpg/util/pml.cpp" +#include "miniapps/dpg/util/blockstaticcond.cpp" +#include "miniapps/dpg/util/complexstaticcond.cpp" +#include "miniapps/dpg/util/weakform.cpp" +#include "miniapps/dpg/util/complexweakform.cpp" +%} + +%include "../common/dpg_common.i" + +%include "miniapps/dpg/util/pml.hpp" +%include "miniapps/dpg/util/blockstaticcond.hpp" +%include "miniapps/dpg/util/complexstaticcond.hpp" +%include "miniapps/dpg/util/weakform.hpp" +%include "miniapps/dpg/util/complexweakform.hpp" + +#endif diff --git a/mfem/_ser/fe_coll.i b/mfem/_ser/fe_coll.i index 38a8fad4..69037c80 100644 --- a/mfem/_ser/fe_coll.i +++ b/mfem/_ser/fe_coll.i @@ -35,9 +35,12 @@ import_array1(-1); %import "lininteg.i" %import "../common/exception.i" - //%inline %{ - // typedef mfem::L2_FECollection mfem::DG_FECollection; - // %} +/* define FiniteElementCollectionPtrArray */ +%import "../common/array_listtuple_typemap.i" +ARRAY_LISTTUPLE_INPUT_SWIGOBJ(mfem::FiniteElementCollection *, 1) +%import "../common/array_instantiation_macro.i" +IGNORE_ARRAY_METHODS(mfem::FiniteElementCollection *) +INSTANTIATE_ARRAY0(FiniteElementCollection *, FiniteElementCollection, 1) %include "fem/fe_coll.hpp" %pythoncode %{ diff --git a/mfem/_ser/handle.i b/mfem/_ser/handle.i index 2c6d086a..0b6531b8 100644 --- a/mfem/_ser/handle.i +++ b/mfem/_ser/handle.i @@ -56,6 +56,13 @@ GET_RENAME(SparseMatrix) RESET_RENAME(SparseMatrix) CONVERT_FROM_RENAME(SparseMatrix) +//CONSTRUCTOR_WRAP(ComplexOperator) +AS_RENAME(ComplexOperator) +IS_RENAME(ComplexOperator) +GET_RENAME(ComplexOperator) +RESET_RENAME(ComplexOperator) +CONVERT_FROM_RENAME(ComplexOperator) + %include "linalg/handle.hpp" %pythoncode %{ @@ -70,4 +77,11 @@ GET_WRAP(SparseMatrix) RESET_WRAP(SparseMatrix) CONVERT_FROM_WRAP(SparseMatrix) +//CONSTRUCTOR_WRAP(ComplexOperator) +AS_WRAP(ComplexOperator) +IS_WRAP(ComplexOperator) +GET_WRAP(ComplexOperator) +RESET_WRAP(ComplexOperator) +CONVERT_FROM_WRAP(ComplexOperator) + diff --git a/mfem/_ser/ordering.i b/mfem/_ser/ordering.i index 2f1c8363..0b793361 100644 --- a/mfem/_ser/ordering.i +++ b/mfem/_ser/ordering.i @@ -13,9 +13,10 @@ import_array1(-1); %} %include "exception.i" -%import "element.i" %include "../common/typemap_macros.i" %include "../common/exception.i" +%import "array.i" +%import "vector.i" %include "linalg/ordering.hpp" diff --git a/mfem/_ser/particleset.i b/mfem/_ser/particleset.i new file mode 100644 index 00000000..16dfa59f --- /dev/null +++ b/mfem/_ser/particleset.i @@ -0,0 +1,31 @@ +// +// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved. +// +%module(package="mfem._ser") particleset + +%feature("autodoc", "1"); + +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "../common/io_stream.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pyintrules.hpp" +%} + +%init %{ +import_array1(-1); +%} + +%include "exception.i" +%include "../common/typemap_macros.i" +%include "../common/exception.i" + +%import "array.i" +%import "ordering.i" +%import "particlevector.i" + +%include "fem/particleset.hpp" + + + diff --git a/mfem/_ser/particlevector.i b/mfem/_ser/particlevector.i new file mode 100644 index 00000000..4216306e --- /dev/null +++ b/mfem/_ser/particlevector.i @@ -0,0 +1,28 @@ +// +// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved. +// +%module(package="mfem._ser") particlevector + +%feature("autodoc", "1"); + +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "../common/io_stream.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pyintrules.hpp" +%} + +%init %{ +import_array1(-1); +%} + +%include "exception.i" +%include "../common/typemap_macros.i" +%include "../common/exception.i" + +%import "vector.i" +%import "ordering.i" + +%include "linalg/particlevector.hpp" + diff --git a/mfem/_ser/pml.i b/mfem/_ser/pml.i new file mode 100644 index 00000000..5b794aec --- /dev/null +++ b/mfem/_ser/pml.i @@ -0,0 +1,44 @@ +%module(package="mfem._ser") pml +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "miniapps/dpg/util/pcomplexweakform.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_PML + +%init %{ +import_array(); +%} + +%inline %{ +#include "miniapps/dpg/util/pml.cpp" +%} + + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "array.i" +%import "vector.i" +%import "gridfunc.i" +%import "mesh.i" +%import "solvers.i" +%import "operators.i" +%import "blockmatrix.i" +%import "complex_densemat.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +%include "miniapps/dpg/util/pml.hpp" + +#endif diff --git a/mfem/_ser/setup.py b/mfem/_ser/setup.py index 65d55150..fd51bba3 100644 --- a/mfem/_ser/setup.py +++ b/mfem/_ser/setup.py @@ -35,10 +35,10 @@ def get_extensions(): # first load variables from PyMFEM_ROOT/setup_local.py sys.path.insert(0, root) try: - from setup_local import (mfemserbuilddir, mfemserincdir, mfemsrcdir, mfemserlnkdir, - mfemstpl, numpyinc, + from setup_local import (mfemserbuilddir, mfemserincdir, mfemsrcdir, + mfemserlnkdir, mfemstpl, numpyinc, cc_ser, cxx_ser, - cxxstdflag, mfem_outside, + cxxstdflag, mfem_outside, build_miniapps, add_cuda, add_libceed, add_suitesparse, add_gslibs, bdist_wheel_dir) @@ -58,11 +58,14 @@ def get_extensions(): add_libceed = '' add_suitesparse = '' add_gslibs = '' - cxxstdflag = '-std=c++17' + cxxstdflag = '-std=c++17' mfem_outside = '0' + build_miniapps = '0' libraries = ['mfem'] + if build_miniapps != '0': + libraries.append("mfem-common") # remove current directory from path # print("__file__", os.path.abspath(__file__)) @@ -108,9 +111,9 @@ def get_extensions(): "submesh", "transfermap", "staticcond", "sidredatacollection", "enzyme", "attribute_sets", "arrays_by_name", - "hyperbolic", - "complex_densemat", "complexstaticcond", "complexweakform", - "bounds", "integrator", "ordering"] + "hyperbolic", "complex_densemat", + "bounds", "integrator", "ordering", + "dpg", "particleset", "particlevector"] if add_cuda == '1': from setup_local import cudainc diff --git a/mfem/_ser/sparsemat.i b/mfem/_ser/sparsemat.i index d3d0de35..6363767b 100644 --- a/mfem/_ser/sparsemat.i +++ b/mfem/_ser/sparsemat.i @@ -144,6 +144,28 @@ if len(args) == 1 and isinstance(args[0], csr_matrix): } } +%typemap(out) PyObject* mfem::SparseMatrix::GetDataArray() { + // assign self to base object + Py_INCREF($self); + PyArray_SetBaseObject((PyArrayObject *) $1, $self); + + $result = $1; +} +%typemap(out) PyObject* mfem::SparseMatrix::GetJArray() { + // assign self to base object + Py_INCREF($self); + PyArray_SetBaseObject((PyArrayObject *) $1, $self); + + $result = $1; +} +%typemap(out) PyObject* mfem::SparseMatrix::GetIArray() { + // assign self to base object + Py_INCREF($self); + PyArray_SetBaseObject((PyArrayObject *) $1, $self); + + $result = $1; +} + %include "linalg/sparsemat.hpp" diff --git a/mfem/_ser/vector.i b/mfem/_ser/vector.i index d72e8a05..4d9822dd 100644 --- a/mfem/_ser/vector.i +++ b/mfem/_ser/vector.i @@ -174,6 +174,14 @@ XXXPTR_SIZE_IN(mfem::Vector **data_, int asize, mfem::Vector *) IGNORE_ARRAY_METHODS(mfem::Vector *) INSTANTIATE_ARRAY0(Vector *, Vector, 1) +%typemap(out) PyObject* mfem::Vector::GetDataArray() { + // assign self to base object + Py_INCREF($self); + PyArray_SetBaseObject((PyArrayObject *) $1, $self); + + $result = $1; +} + %include "linalg/vector.hpp" %extend mfem::Vector { diff --git a/mfem/_ser/weakform.i b/mfem/_ser/weakform.i new file mode 100644 index 00000000..dddd70c6 --- /dev/null +++ b/mfem/_ser/weakform.i @@ -0,0 +1,68 @@ +%module(package="mfem._ser") weakform +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "miniapps/dpg/util/weakform.hpp" +#include "../common/pyoperator.hpp" +#include "../common/pysolvers.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pyintrules.hpp" +#include "../common/pylininteg.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" +#include "../common/io_stream.hpp" +%} + + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_MINIAPPS_DPG_UTIL_WEAKFORM + +%init %{ +import_array(); +%} + +%inline %{ +#include "miniapps/dpg/util/weakform.cpp" +%} + + +%include "exception.i" +%import "element.i" +%import "../common/exception.i" + +%import "array.i" +%import "vector.i" +%import "blockvector.i" +%import "bilininteg.i" +%import "lininteg.i" +%import "fe_coll.i" +%import "blockmatrix.i" +%import "operators.i" +%import "../common/exception.i" +%import "../common/io_stream_typemap.i" + +OSTREAM_TYPEMAP(std::ostream&) + +%import "../common/object_array_typemap.i" +LIST_TO_MFEMOBJ_ARRAY_IN(mfem::Array&, + mfem::FiniteElementSpace*) +LIST_TO_MFEMOBJ_ARRAY_IN(mfem::Array&, + mfem::FiniteElementCollection*) + + +%pythonprepend mfem::DPGWeakForm::DPGWeakForm %{ + if len(args) > 0: + fes_, fecol_ = args + self._fes = fes_ + self._fecol = fecol_ +%} + +%pythonprepend mfem::DPGWeakForm::SetSpaces %{ + self._fes = fes_ + self._fecol = fecol_ +%} + + +%include "miniapps/dpg/util/weakform.hpp" + +#endif diff --git a/mfem/common/dpg_common.i b/mfem/common/dpg_common.i new file mode 100644 index 00000000..bbdbab7f --- /dev/null +++ b/mfem/common/dpg_common.i @@ -0,0 +1,166 @@ +%import "element.i" +%import "array.i" +%import "vector.i" +%import "blockvector.i" +%import "bilininteg.i" +%import "lininteg.i" +%import "fespace.i" +%import "blockmatrix.i" +%import "operators.i" +%import "blockmatrix.i" +%import "blockoperator.i" +%import "complex_densemat.i" +%import "../common/io_stream_typemap.i" + +OSTREAM_TYPEMAP(std::ostream&) + +%ignore mfem::BlockStaticCondensation::ConvertListToReducedTrueDofs; +%ignore mfem::ComplexBlockStaticCondensation::ConvertListToReducedTrueDofs; + +namespace mfem { +%pythonprepend DPGWeakForm::AddTrialIntegrator %{ + if not hasattr(self, "_integrators"): self._integrators = [] + di = bfi + self._integrators.append(di) + di.thisown=0 ## ownership is C++ side + %} +%pythonprepend DPGWeakForm::AddTestIntegrator %{ + if not hasattr(self, "_integrators"): self._integrators = [] + di = bfi + self._integrators.append(di) + di.thisown=0 ## ownership is C++ side + %} +%pythonprepend DPGWeakForm::AddDomainLFIntegrator %{ + if not hasattr(self, "_integrators"): self._integrators = [] + di = lfi + self._integrators.append(di) + di.thisown=0 ## ownership is C++ side + %} +%pythonprepend ComplexDPGWeakForm::AddTrialIntegrator %{ + if not hasattr(self, "_integrators"): self._integrators = [] + di = bfi_r + if di is not None: + self._integrators.append(di) + di.thisown=0 ## ownership is C++ side + di = bfi_i + if di is not None: + self._integrators.append(di) + di.thisown=0 ## ownership is C++ side + %} +%pythonprepend ComplexDPGWeakForm::AddTestIntegrator %{ + if not hasattr(self, "_integrators"): self._integrators = [] + di = bfi_r + if di is not None: + self._integrators.append(di) + di.thisown=0 ## ownership is C++ side + di = bfi_i + if di is not None: + self._integrators.append(di) + di.thisown=0 ## ownership is C++ side + %} +%pythonprepend ComplexDPGWeakForm::AddDomainLFIntegrator %{ + if not hasattr(self, "_integrators"): self._integrators = [] + di = lfi_r + if di is not None: + self._integrators.append(di) + di.thisown=0 ## ownership is C++ side + di = lfi_i + if di is not None: + self._integrators.append(di) + di.thisown=0 ## ownership is C++ side + %} +} + +//%constant mfem::real_t detJ_r_function(const mfem::Vector &, mfem::CartesianPML *); + +%constant mfem::real_t detJ_r_function(const mfem::Vector &, mfem::CartesianPML *); +%constant mfem::real_t detJ_i_function(const mfem::Vector &, mfem::CartesianPML *); +%constant mfem::real_t abs_detJ_2_function(const mfem::Vector &, mfem::CartesianPML *); +%constant void Jt_J_detJinv_r_function(const mfem::Vector & , mfem::CartesianPML *, + mfem::DenseMatrix & ); +%constant void Jt_J_detJinv_i_function(const mfem::Vector & , mfem::CartesianPML *, + mfem::DenseMatrix &); +%constant void abs_Jt_J_detJinv_2_function(const mfem::Vector &, mfem::CartesianPML *, + mfem::DenseMatrix &); +%constant void detJ_Jt_J_inv_r_function(const mfem::Vector &, mfem::CartesianPML *, + mfem::DenseMatrix &); +%constant void detJ_Jt_J_inv_i_function(const mfem::Vector &, mfem::CartesianPML *, + mfem::DenseMatrix &); +%constant void abs_detJ_Jt_J_inv_2_function(const mfem::Vector &, mfem::CartesianPML *, + mfem::DenseMatrix &); + + +%rename(detJ_r_f) mfem::detJ_r_function; +%rename(detJ_i_f) mfem::detJ_i_function; +%rename(abs_detJ_2_f) mfem::abs_detJ_2_function; +%rename(Jt_J_detJinv_r_f) mfem::Jt_J_detJinv_r_function; +%rename(Jt_J_detJinv_i_f) mfem::Jt_J_detJinv_i_function; +%rename(abs_Jt_J_detJinv_2_f) mfem::abs_Jt_J_detJinv_2_function; +%rename(detJ_Jt_J_inv_r_f) mfem::detJ_Jt_J_inv_r_function; +%rename(detJ_Jt_J_inv_i_f) mfem::detJ_Jt_J_inv_i_function; +%rename(abs_detJ_Jt_J_inv_2_f) mfem::abs_detJ_Jt_J_inv_2_function; + + +// % handle Array input +%include "../common/typemap_macros.i" +LIST_TO_MFEMOBJ_ARRAY_IN(mfem::Array&, + mfem::FiniteElementSpace*) +LIST_TO_MFEMOBJ_ARRAY_IN(mfem::Array&, + mfem::ParFiniteElementSpace*) +LIST_TO_MFEMOBJ_ARRAY_IN(mfem::Array&, + mfem::FiniteElementCollection*) + + +//DPGWeakForm +%pythonprepend mfem::DPGWeakForm::DPGWeakForm %{ + if len(args) > 0: + fes_, fecol_ = args + self._fes = fes_ + self._fecol = fecol_ +%} + +%pythonprepend mfem::DPGWeakForm::SetSpaces %{ + self._fes = fes_ + self._fecol = fecol_ +%} + +//ComplexDPGWeakForm +%pythonprepend mfem::ComplexDPGWeakForm::ComplexDPGWeakForm %{ + if len(args) > 0: + fes_, fecol_ = args + self._fes = fes_ + self._fecol = fecol_ +%} + +%pythonprepend mfem::ComplexDPGWeakForm::SetSpaces %{ + self._fes = fes_ + self._fecol = fecol_ +%} + +//ParDPGWeakForm +%pythonprepend mfem::ParDPGWeakForm::ParDPGWeakForm %{ + if len(args) > 0: + trial_pfes_, fecol_ = args + self._fes = trial_pfes_ + self._fecol = fecol_ + +%} + +%pythonprepend mfem::ParDPGWeakForm::SetParSpaces %{ + self._pfes = trial_pfes_ + self._fecol = fecol_ +%} + +//ParComplexDPGWeakForm +%pythonprepend mfem::ParComplexDPGWeakForm::ParComplexDPGWeakForm %{ + if len(args) > 0: + trial_pfes_, fecol_ = args + self._fes = trial_pfes_ + self._fecol = fecol_ + +%} + +%pythonprepend mfem::ParComplexDPGWeakForm::SetParSpaces %{ + self._pfes = trial_pfes_ + self._fecol = fecol_ +%} diff --git a/mfem/common/mfem_config.i b/mfem/common/mfem_config.i index 2e9ccab3..2047d774 100644 --- a/mfem/common/mfem_config.i +++ b/mfem/common/mfem_config.i @@ -2,7 +2,6 @@ // Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved. // %ignore MFEM_GIT_STRING; -%ignore MFEM_HYPRE_VERSION; %ignore MFEM_SOURCE_DIR; %ignore MFEM_INSTALL_DIR; %ignore MFEM_TIMER_TYPE; diff --git a/mfem/common/numba_coefficient.i b/mfem/common/numba_coefficient.i index c262a0a7..6b031274 100644 --- a/mfem/common/numba_coefficient.i +++ b/mfem/common/numba_coefficient.i @@ -1422,12 +1422,12 @@ def _matrix(func, height=None, width=None, shape=None, td=False, params=None, if height is None and width is None and shape is not None: width = shape[0] height = shape[1] - assert height == shape[0], "height and shape[0] are not consistent" - assert width == shape[1], "width and shape[1] are not consistent" + assert height == shape[1], "height and shape[0] are not consistent" + assert width == shape[0], "width and shape[1] are not consistent" - if shape[0] != shape[1]: - import warnings - warning.warn("Rectangular matrix coefficient is experimental", UserWarning) + #if shape[0] != shape[1]: + # import warnings + # warnings.warn("Rectangular matrix coefficient is experimental", UserWarning) if dependency is None: dependency = [] diff --git a/mfem/common/numba_coefficient_utils.py b/mfem/common/numba_coefficient_utils.py index ec9024fe..e4f3867c 100644 --- a/mfem/common/numba_coefficient_utils.py +++ b/mfem/common/numba_coefficient_utils.py @@ -445,7 +445,7 @@ def _process_dependencies(dependencies, setting): if isinstance(xx, MatrixCoefficient): h1, w1 = x[0].GetHeight(), x[0].GetWidth() h2, w2 = x[1].GetHeight(), x[1].GetWidth() - assert h1 == w1, "matrix must be square" + #assert h1 == w1, "matrix must be square" assert h1 == h2, "real and imaginary has to have the same vdim" assert w1 == w2, "real and imaginary has to have the same vdim" else: @@ -461,8 +461,8 @@ def _process_dependencies(dependencies, setting): else: iscomplex.append(0) - if isinstance(xx, MatrixCoefficient): - assert xx.GetHeight() == xx.GetWidth(), "matrix must be square" + #if isinstance(xx, MatrixCoefficient): + # assert xx.GetHeight() == xx.GetWidth(), "matrix must be square" if isinstance(xx, Coefficient): kinds.append(0) diff --git a/mfem/par.py b/mfem/par.py index 95523511..4841966c 100644 --- a/mfem/par.py +++ b/mfem/par.py @@ -34,6 +34,7 @@ from mfem._par.fe_coll import * from mfem._par.vector import * from mfem._par.complex_operator import * +from mfem._par.complex_densemat import * from mfem._par.complex_fem import * from mfem._par.fespace import * from mfem._par.linearform import * @@ -79,6 +80,8 @@ from mfem._par.quadinterpolator_face import * from mfem._par.attribute_sets import * from mfem._par.ordering import * +from mfem._par.particleset import * +from mfem._par.particlevector import * from mfem._par.fe_base import * from mfem._par.fe_h1 import * @@ -129,6 +132,17 @@ # # modules not a part of standard build # +import importlib.util +import sys +import types +def load_module(module_name, module_code): + spec = importlib.util.spec_from_loader(module_name, loader=None, origin="dynamic") + module = types.ModuleType(module_name) + spec.loader.exec_module(module) if spec.loader else exec(module_code, module.__dict__) + sys.modules[module_name] = module + return module + + try: import mfem._par.pumi as pumi from mfem._par.pumi import * @@ -155,7 +169,14 @@ pass +try: + import mfem._par.dpg as dpg +except: + pass + + # # initialize hypre # Hypre.Init() + diff --git a/mfem/ser.py b/mfem/ser.py index 4c5a7895..47e2be48 100644 --- a/mfem/ser.py +++ b/mfem/ser.py @@ -25,6 +25,7 @@ from mfem._ser.fe_coll import * from mfem._ser.vector import * from mfem._ser.complex_operator import * +from mfem._ser.complex_densemat import * from mfem._ser.complex_fem import * from mfem._ser.fespace import * from mfem._ser.linearform import * @@ -62,6 +63,8 @@ from mfem._ser.quadinterpolator_face import * from mfem._ser.attribute_sets import * from mfem._ser.ordering import * +from mfem._ser.particleset import * +from mfem._ser.particlevector import * from mfem._ser.fe_base import * from mfem._ser.fe_h1 import * @@ -74,15 +77,9 @@ from mfem._ser.fe_nurbs import * from mfem._ser.doftrans import * from mfem._ser.std_vectors import * - from mfem._ser.bounds import * from mfem._ser.integrator import * -try: - from mfem._ser.complex_densemat import * -except ImportError: - pass - from mfem._ser.submesh import * from mfem._ser.transfermap import * from mfem._ser.hyperbolic import * @@ -115,9 +112,11 @@ # # import moduleds built from hpp/cpp under miniapp # + try: - import mfem._ser.complexweakform as complexweakform - import mfem._ser.commlexstaticcond as complexstaticcond + import mfem._ser.dpg as dpg except: pass + + diff --git a/miniapps/dpg/meshes/fichera-waveguide.mesh b/miniapps/dpg/meshes/fichera-waveguide.mesh new file mode 100644 index 00000000..4faba20d --- /dev/null +++ b/miniapps/dpg/meshes/fichera-waveguide.mesh @@ -0,0 +1,93 @@ +MFEM mesh v1.0 + +# +# MFEM Geometry Types (see fem/geom.hpp): +# +# POINT = 0 +# SEGMENT = 1 +# TRIANGLE = 2 +# SQUARE = 3 +# TETRAHEDRON = 4 +# CUBE = 5 +# PRISM = 6 +# PYRAMID = 7 +# + +dimension +3 + +elements +8 +1 5 0 1 4 3 9 10 13 12 +2 5 9 10 13 12 18 19 22 21 +3 5 10 11 14 13 19 20 23 22 +4 5 1 2 5 4 10 11 14 13 +5 5 4 5 8 7 13 14 17 16 +6 5 12 13 16 15 21 22 25 24 +7 5 3 4 7 6 12 13 16 15 +8 5 18 19 22 21 26 27 29 28 + +boundary +28 +1 3 0 3 4 1 +1 3 1 4 5 2 +1 3 3 6 7 4 +1 3 4 7 8 5 +4 3 26 27 29 28 +3 3 0 9 12 3 +3 3 3 12 15 6 +3 3 9 18 21 12 +3 3 12 21 24 15 +3 3 18 26 28 21 +2 3 0 1 10 9 +2 3 9 10 19 18 +2 3 18 19 27 26 +2 3 1 2 11 10 +2 3 10 11 20 19 +5 3 11 14 23 20 +5 3 14 13 22 23 +5 3 19 20 23 22 +5 3 2 5 14 11 +5 3 5 8 17 14 +5 3 8 7 16 17 +5 3 13 14 17 16 +5 3 22 13 16 25 +5 3 16 15 24 25 +5 3 21 22 25 24 +5 3 7 6 15 16 +5 3 29 22 21 28 +5 3 19 22 29 27 + +vertices +30 +3 +0 0 0 +1 0 0 +2 0 0 +0 1 0 +1 1 0 +2 1 0 +0 2 0 +1 2 0 +2 2 0 +0 0 1 +1 0 1 +2 0 1 +0 1 1 +1 1 1 +2 1 1 +0 2 1 +1 2 1 +2 2 1 +0 0 2 +1 0 2 +2 0 2 +0 1 2 +1 1 2 +2 1 2 +0 2 2 +1 2 2 +0 0 3 +1 0 3 +0 1 3 +1 1 3 diff --git a/miniapps/dpg/meshes/scatter.mesh b/miniapps/dpg/meshes/scatter.mesh new file mode 100644 index 00000000..dcc21d89 --- /dev/null +++ b/miniapps/dpg/meshes/scatter.mesh @@ -0,0 +1,247 @@ +MFEM mesh v1.0 + +# +# MFEM Geometry Types (see fem/geom.hpp): +# +# POINT = 0 +# SEGMENT = 1 +# TRIANGLE = 2 +# SQUARE = 3 +# TETRAHEDRON = 4 +# CUBE = 5 +# PRISM = 6 +# PYRAMID = 7 +# + +dimension +2 + +elements +80 +1 3 0 1 11 10 +1 3 10 11 21 20 +1 3 11 12 22 21 +1 3 1 2 12 11 +1 3 2 3 13 12 +1 3 3 4 14 13 +1 3 13 14 24 23 +1 3 12 13 23 22 +1 3 22 23 33 32 +1 3 23 24 34 33 +1 3 33 34 44 43 +1 3 32 33 43 42 +1 3 31 32 42 41 +1 3 21 22 32 31 +1 3 20 21 31 30 +1 3 30 31 41 40 +1 3 40 41 51 50 +1 3 41 42 52 51 +1 3 51 52 62 61 +1 3 50 51 61 60 +1 3 60 61 71 70 +1 3 61 62 72 71 +1 3 71 72 82 81 +1 3 70 71 81 80 +1 3 80 81 91 90 +1 3 81 82 92 91 +1 3 82 83 93 92 +1 3 83 84 94 93 +1 3 73 74 84 83 +1 3 72 73 83 82 +1 3 62 63 73 72 +1 3 63 64 74 73 +1 3 53 54 64 63 +1 3 52 53 63 62 +1 3 42 43 53 52 +1 3 43 44 54 53 +1 3 45 46 56 55 +1 3 55 56 66 65 +1 3 54 55 65 64 +1 3 64 65 75 74 +1 3 74 75 85 84 +1 3 84 85 95 94 +1 3 85 86 96 95 +1 3 75 76 86 85 +1 3 65 66 76 75 +1 3 66 67 77 76 +1 3 76 77 87 86 +1 3 86 87 97 96 +1 3 87 88 98 97 +1 3 88 89 99 98 +1 3 78 79 89 88 +1 3 77 78 88 87 +1 3 67 68 78 77 +1 3 68 69 79 78 +1 3 58 59 69 68 +1 3 57 58 68 67 +1 3 56 57 67 66 +1 3 46 47 57 56 +1 3 47 48 58 57 +1 3 48 49 59 58 +1 3 38 39 49 48 +1 3 28 29 39 38 +1 3 27 28 38 37 +1 3 37 38 48 47 +1 3 36 37 47 46 +1 3 26 27 37 36 +1 3 25 26 36 35 +1 3 35 36 46 45 +1 3 34 35 45 44 +1 3 24 25 35 34 +1 3 14 15 25 24 +1 3 4 5 15 14 +1 3 5 6 16 15 +1 3 15 16 26 25 +1 3 16 17 27 26 +1 3 6 7 17 16 +1 3 7 8 18 17 +1 3 17 18 28 27 +1 3 18 19 29 28 +1 3 8 9 19 18 + +boundary +40 +1 1 0 1 +1 1 1 2 +1 1 2 3 +1 1 3 4 +1 1 4 5 +1 1 5 6 +1 1 6 7 +1 1 7 8 +1 1 8 9 +1 1 91 90 +1 1 92 91 +1 1 93 92 +1 1 94 93 +1 1 95 94 +1 1 96 95 +1 1 97 96 +1 1 98 97 +1 1 99 98 +1 1 10 0 +1 1 20 10 +1 1 30 20 +1 1 40 30 +1 1 50 40 +1 1 60 50 +1 1 70 60 +1 1 80 70 +1 1 90 80 +1 1 9 19 +1 1 19 29 +1 1 29 39 +1 1 39 49 +1 1 49 59 +1 1 59 69 +1 1 69 79 +1 1 79 89 +1 1 89 99 +2 1 44 54 +2 1 45 44 +2 1 55 45 +2 1 54 55 + +vertices +100 +2 +0 0 +0.11111111 0 +0.22222222 0 +0.33333333 0 +0.44444444 0 +0.55555556 0 +0.66666667 0 +0.77777778 0 +0.88888889 0 +1 0 +0 0.11111111 +0.11111111 0.11111111 +0.22222222 0.11111111 +0.33333333 0.11111111 +0.44444444 0.11111111 +0.55555556 0.11111111 +0.66666667 0.11111111 +0.77777778 0.11111111 +0.88888889 0.11111111 +1 0.11111111 +0 0.22222222 +0.11111111 0.22222222 +0.22222222 0.22222222 +0.33333333 0.22222222 +0.44444444 0.22222222 +0.55555556 0.22222222 +0.66666667 0.22222222 +0.77777778 0.22222222 +0.88888889 0.22222222 +1 0.22222222 +0 0.33333333 +0.11111111 0.33333333 +0.22222222 0.33333333 +0.33333333 0.33333333 +0.44444444 0.33333333 +0.55555556 0.33333333 +0.66666667 0.33333333 +0.77777778 0.33333333 +0.88888889 0.33333333 +1 0.33333333 +0 0.44444444 +0.11111111 0.44444444 +0.22222222 0.44444444 +0.33333333 0.44444444 +0.44444444 0.44444444 +0.55555556 0.44444444 +0.66666667 0.44444444 +0.77777778 0.44444444 +0.88888889 0.44444444 +1 0.44444444 +0 0.55555556 +0.11111111 0.55555556 +0.22222222 0.55555556 +0.33333333 0.55555556 +0.44444444 0.55555556 +0.55555556 0.55555556 +0.66666667 0.55555556 +0.77777778 0.55555556 +0.88888889 0.55555556 +1 0.55555556 +0 0.66666667 +0.11111111 0.66666667 +0.22222222 0.66666667 +0.33333333 0.66666667 +0.44444444 0.66666667 +0.55555556 0.66666667 +0.66666667 0.66666667 +0.77777778 0.66666667 +0.88888889 0.66666667 +1 0.66666667 +0 0.77777778 +0.11111111 0.77777778 +0.22222222 0.77777778 +0.33333333 0.77777778 +0.44444444 0.77777778 +0.55555556 0.77777778 +0.66666667 0.77777778 +0.77777778 0.77777778 +0.88888889 0.77777778 +1 0.77777778 +0 0.88888889 +0.11111111 0.88888889 +0.22222222 0.88888889 +0.33333333 0.88888889 +0.44444444 0.88888889 +0.55555556 0.88888889 +0.66666667 0.88888889 +0.77777778 0.88888889 +0.88888889 0.88888889 +1 0.88888889 +0 1 +0.11111111 1 +0.22222222 1 +0.33333333 1 +0.44444444 1 +0.55555556 1 +0.66666667 1 +0.77777778 1 +0.88888889 1 +1 1 diff --git a/miniapps/dpg/pmaxwell.py b/miniapps/dpg/pmaxwell.py new file mode 100644 index 00000000..124d553d --- /dev/null +++ b/miniapps/dpg/pmaxwell.py @@ -0,0 +1,1161 @@ +## +# MFEM Ultraweak DPG Maxwell parallel example +## +# Compile with: make pmaxwell +## +# sample run + +# mpirun -np 1 python pmaxwell.py -o 2 -pref 2 -prob 0 +# mpirun -np 4 python pmaxwell.py -o 3 -sref 0 -pref 3 -rnum 4.8 -sc -prob 0 +# mpirun -np 1 python pmaxwell.py -o 2 -pref 2 -prob 0 -m inline-hex.mesh + +# mpirun -np 4 python pmaxwell.py -m ../../data/star.mesh -o 2 -sref 0 -pref 3 -rnum 0.5 -prob 0 +# mpirun -np 4 python pmaxwell.py -m ../../data/inline-quad.mesh -o 3 -sref 0 -pref 3 -rnum 4.8 -sc -prob 0 +# mpirun -np 4 python pmaxwell.py -m ../../data/inline-hex.mesh -o 2 -sref 0 -pref 1 -rnum 0.8 -sc -prob 0 +# mpirun -np 4 python pmaxwell.py -m ../../data/inline-quad.mesh -o 3 -sref 1 -pref 3 -rnum 4.8 -sc -prob 2 +# mpirun -np 4 python pmaxwell.py -o 3 -sref 1 -pref 2 -rnum 11.8 -sc -prob 3 +# mpirun -np 4 python pmaxwell.py -o 3 -sref 1 -pref 2 -rnum 9.8 -sc -prob 4 + +# AMR run. Note that this is a computationally intensive sample run. +# We recommend trying it on a large machine with more mpi ranks +# mpirun -np 4 pmaxwell -o 3 -sref 0 -pref 15 -prob 1 -theta 0.7 -sc + +# Description: +# This example code demonstrates the use of MFEM to define and solve +# the "ultraweak" (UW) DPG formulation for the Maxwell problem + +# ∇×(1/μ ∇×E) - ω² ϵ E = Ĵ , in Ω +# E×n = E₀ , on ∂Ω + +# It solves the following kinds of problems +# 1) Known exact solutions with error convergence rates +# a) A manufactured solution problem where E is a plane beam +# 2) Fichera "microwave" problem +# 3) PML problems +# a) Generic PML problem with point source given by the load +# b) Plane wave scattering from a square +# c) PML problem with a point source prescribed on the boundary + +# The DPG UW deals with the First Order System +# i ω μ H + ∇ × E = 0, in Ω +# -i ω ϵ E + ∇ × H = J, in Ω +# E × n = E_0, on ∂Ω +# Note: Ĵ = -iωJ + +# The ultraweak-DPG formulation is obtained by integration by parts of both +# equations and the introduction of trace unknowns on the mesh skeleton + +# in 2D +# E is vector valued and H is scalar. +# (∇ × E, F) = (E, ∇ × F) + < n × E , F> +# or (∇ ⋅ AE , F) = (AE, ∇ F) + < AE ⋅ n, F> +# where A = [0 1; -1 0]; + +# E ∈ (L²(Ω))² , H ∈ L²(Ω) +# Ê ∈ H^-1/2(Γₕ), Ĥ ∈ H^1/2(Γₕ) +# i ω μ (H,F) + (E, ∇ × F) + < AÊ, F > = 0, ∀ F ∈ H¹ +# -i ω ϵ (E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω) +# Ê = E₀ on ∂Ω +# ------------------------------------------------------------------------- +# | | E | H | Ê | Ĥ | RHS | +# ------------------------------------------------------------------------- +# | F | (E,∇ × F) | i ω μ (H,F) | < Ê, F > | | | +# | | | | | | | +# | G | -i ω ϵ (E,G) | (H,∇ × G) | | < Ĥ, G × n > | (J,G) | +# where (F,G) ∈ H¹ × H(curl,Ω) + +# in 3D +# E,H ∈ (L^2(Ω))³ +# Ê ∈ H_0^1/2(Ω)(curl, Γₕ), Ĥ ∈ H^-1/2(curl, Γₕ) +# i ω μ (H,F) + (E,∇ × F) + < Ê, F × n > = 0, ∀ F ∈ H(curl,Ω) +# -i ω ϵ (E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω) +# Ê × n = E₀ on ∂Ω +# ------------------------------------------------------------------------- +# | | E | H | Ê | Ĥ | RHS | +# ------------------------------------------------------------------------- +# | F | (E,∇ × F) | i ω μ (H,F) | < n × Ê, F > | | | +# | | | | | | | +# | G | -i ω ϵ (E,G) | (H,∇ × G) | | < n × Ĥ, G > | (J,G) | +# where (F,G) ∈ H(curl,Ω) × H(curl,Ω) + +# Here we use the "Adjoint Graph" norm on the test space i.e., +# ||(F,G)||²ᵥ = ||A^*(F,G)||² + ||(F,G)||² where A is the +# maxwell operator defined by (1) + +# The PML formulation is + +# ∇×(1/μ α ∇×E) - ω² ϵ β E = Ĵ , in Ω +# E×n = E₀ , on ∂Ω + +# where α = |J|⁻¹ Jᵀ J (in 2D it's the scalar |J|⁻¹), +# β = |J| J⁻¹ J⁻ᵀ, J is the Jacobian of the stretching map +# and |J| its determinant. + +# The first order system reads +# i ω μ α⁻¹ H + ∇ × E = 0, in Ω +# -i ω ϵ β E + ∇ × H = J, in Ω +# E × n = E₀, on ∂Ω + +# and the ultraweak formulation is + +# in 2D +# E ∈ (L²(Ω))² , H ∈ L²(Ω) +# Ê ∈ H^-1/2(Ω)(Γₕ), Ĥ ∈ H^1/2(Γₕ) +# i ω μ (α⁻¹ H,F) + (E, ∇ × F) + < AÊ, F > = 0, ∀ F ∈ H¹ +# -i ω ϵ (β E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω) +# Ê = E₀ on ∂Ω +# --------------------------------------------------------------------------------- +# | | E | H | Ê | Ĥ | RHS | +# --------------------------------------------------------------------------------- +# | F | (E,∇ × F) | i ω μ (α⁻¹ H,F) | < Ê, F > | | | +# | | | | | | | +# | G | -i ω ϵ (β E,G) | (H,∇ × G) | | < Ĥ, G × n > | (J,G) | + +# where (F,G) ∈ H¹ × H(curl,Ω) + +## +# in 3D +# E,H ∈ (L^2(Ω))³ +# Ê ∈ H_0^1/2(Ω)(curl, Γ_h), Ĥ ∈ H^-1/2(curl, Γₕ) +# i ω μ (α⁻¹ H,F) + (E,∇ × F) + < Ê, F × n > = 0, ∀ F ∈ H(curl,Ω) +# -i ω ϵ (β E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω) +# Ê × n = E_0 on ∂Ω +# ------------------------------------------------------------------------------- +# | | E | H | Ê | Ĥ | RHS | +# ------------------------------------------------------------------------------- +# | F | ( E,∇ × F) | i ω μ (α⁻¹ H,F) | < n × Ê, F > | | | +# | | | | | | | +# | G | -iωϵ (β E,G) | (H,∇ × G) | | < n × Ĥ, G > | (J,G) | +# where (F,G) ∈ H(curl,Ω) × H(curl,Ω) + +# For more information see https:##doi.org/10.1016/j.camwa.2021.01.017 +from mpi4py import MPI +from numba import njit, void, int32, int64, float64, complex128, types +from mfem.common.bessel import yv as yn +from mfem.common.bessel import jv as jn +import os +from os.path import expanduser, join, dirname + +import numpy as np +from numpy import pi, exp + +import mfem.par as mfem +mdpg = mfem.dpg + + +num_procs = MPI.COMM_WORLD.size +myid = MPI.COMM_WORLD.rank +smyid = '.'+'{:0>6d}'.format(myid) + +visport = 19916 + + +def VisualizeField(sock, vishost, visport, + gf, title, + x=0, y=0, w=400, h=400, keys='', vec=False): + + pmesh = gf.ParFESpace().GetParMesh() + + newly_opened = False + + while True: + if sock is None: + sock = mfem.socketstream("localhost", visport) + sock.precision(8) + newly_opened = True + + sock << "parallel " << str(num_procs) << " " << str(myid) << "\n" + sock << "solution\n" + sock << pmesh << gf + + if newly_opened: + sock << "window_title '" << title << "'\n" + sock << "window_geometry " + sock << x << " " << y << " " << w << " " << h << "\n" + + if keys != '': + sock << "keys " << keys << "\n" + else: + sock << "keys maaAc\n" + if vec: + sock << "vvv" + sock.flush() + if sock.is_open(): + break + + return sock + + +prob_type = ("plane_wave", + "fichera_oven", + "pml_general", + "pml_plane_wave_scatter", + "pml_pointsource") + + +class Functions: + def __init__(self, dim, dimc, omega, epsilon, mu, iprob): + + prob = prob_type[iprob] + + @njit(complex128[:](float64[:])) + def maxwell_solution(x): + E = np.zeros(dim, dtype=np.complex128) + + if prob == "plane_wave": + E[0] = exp(1j * omega * np.sum(x)) + + elif prob == "pml_plane_wave_scatter": + E[1] = exp(1j * omega * x[0]) + + elif prob == "fichera_oven": + if abs(x[2] - 3.0) < 1e-10: + E[0] = np.sin(pi*x[1]) + elif prob == "pml_pointsource": + k = omega * np.sqrt(epsilon * mu) + shift = np.array([-0.5]*dim) + + if dim == 2: + x0 = x[0] + shift[0] + x1 = x[1] + shift[1] + r = np.sqrt(x0 * x0 + x1 * x1) + beta = k * r + + # Bessel functions + Ho = jn(0, beta).real + 1j * yn(0, beta).real + Ho_r = -k * jn(1, beta).real + 1j * yn(1, beta).real + Ho_rr = -k * k * (1. / beta * + (jn(1, beta).real + 1j * yn(1, beta).real) - + (jn(2, beta).real + 1j * yn(2, beta).real)) + + # First derivatives + r_x = x0 / r + r_y = x1 / r + r_xy = -(r_x / r) * r_y + r_xx = (1.0 / r) * (1.0 - r_x * r_x) + + val = 0.25 * 1j * Ho + val_xx = 0.25 * 1j * (r_xx * Ho_r + r_x * r_x * Ho_rr) + val_xy = 0.25 * 1j * (r_xy * Ho_r + r_x * r_y * Ho_rr) + E[0] = 1j / k * (k * k * val + val_xx) + E[1] = 1j / k * val_xy + else: + x0 = x[0] + shift[0] + x1 = x[1] + shift[1] + x2 = x[2] + shift[2] + r = sqrt(x0 * x0 + x1 * x1 + x2 * x2) + + r_x = x0 / r + r_y = x1 / r + r_z = x2 / r + r_xx = (1.0 / r) * (1.0 - r_x * r_x) + r_yx = -(r_y / r) * r_x + r_zx = -(r_z / r) * r_x + + val = exp(1j * k * r) / r + val_r = val / r * (1j * k * r - 1.) + val_rr = val / (r * r) * (-k * k * r * r + - 2. * 1j * k * r + 2.) + + val_xx = val_rr * r_x * r_x + val_r * r_xx + val_yx = val_rr * r_x * r_y + val_r * r_yx + val_zx = val_rr * r_x * r_z + val_r * r_zx + alpha = 1j * k / 4. / pi / k / k + E[0] = alpha * (k * k * val + val_xx) + E[1] = alpha * val_yx + E[2] = alpha * val_zx + else: + assert False, "should not come here" + return E + + @njit(complex128[:](float64[:])) + def maxwell_solution_curl(x): + curlE = np.zeros(dimc, dtype=np.complex128) + if prob == "plane_wave": + pw = exp(1j * omega * np.sum(x)) + if dim == 3: + curlE[0] = 0.0 + curlE[1] = 1j * omega * pw + curlE[2] = -1j * omega * pw + else: + curlE[0] = -1j * omega * pw + elif prob == "pml_plane_wave_scatter": + pw = exp(1j * omega * (x[0])) + curlE[0] = 1j * omega * pw + else: + assert False, "should not come here" + return curlE + + @njit(complex128[:](float64[:])) + def maxwell_solution_curlcurl(x): + curlcurlE = np.zeros(dim, dtype=np.complex128) + if prob == "plane_wave": + pw = exp(1j * omega * np.sum(x)) + if dim == 3: + curlcurlE[0] = 2. * omega * omega * pw + curlcurlE[1] = - omega * omega * pw + curlcurlE[2] = - omega * omega * pw + else: + curlcurlE[0] = omega * omega * pw + curlcurlE[1] = -omega * omega * pw + elif prob == "pml_plane_wave_scatter": + pw = np.exp(1j * omega * x[0]) + curlcurlE[1] = omega * omega * pw + else: + assert False, "should not come here" + return curlcurlE + + self.maxwell_solution = maxwell_solution + self.maxwell_solution_curl = maxwell_solution_curl + self.maxwell_solution_curlcurl = maxwell_solution_curlcurl + + @njit(complex128[:](float64[:])) + def curlH_exact(x): + # ∇ × H = ∇ × ∇ × E / ω μ + curlcurlE = maxwell_solution_curlcurl(x) + return 1j*curlcurlE / (omega * mu) + + @njit(complex128[:](float64[:])) + def H_exact(x): + H = maxwell_solution_curl(x)*1j/omega/mu + return H + + @njit(complex128[:](float64[:])) + def hatE_exact(x): + if dim == 3: + E = maxwell_solution(x) + return E + else: + ret = np.zeros(2, dtype=np.complex128) + E = maxwell_solution(x) + ret[0] = E[1] + ret[1] = -E[0] + return ret + + @mfem.jit.vector(vdim=dim) + def E_exact_r(x): + E = maxwell_solution(x) + return E.real + + @mfem.jit.vector(vdim=dim) + def E_exact_i(x): + E = maxwell_solution(x) + return E.imag + + @mfem.jit.vector(vdim=dim) + def H_exact_r(x): + return H_exact(x).real + + @mfem.jit.vector(vdim=dim) + def H_exact_i(x): + return H_exact(x).imag + + @mfem.jit.vector(vdim=dim) + def hatE_exact_r(x): + ret = hatE_exact(x) + return ret.real + + @mfem.jit.vector(vdim=dim) + def hatE_exact_i(x): + ret = hatE_exact(x) + return ret.imag + + self.E_exact_r = E_exact_r + self.E_exact_i = E_exact_i + self.H_exact_r = H_exact_r + self.H_exact_i = H_exact_i + self.hatE_exact_r = hatE_exact_r + self.hatE_exact_i = hatE_exact_i + + # J = -i ω ϵ E + ∇ × H + # J_r + iJ_i = -i ω ϵ (E_r + i E_i) + ∇ × (H_r + i H_i) + @mfem.jit.vector(vdim=dim) + def rhs_func_r(x): + # J_r = ω ϵ E_i + ∇ × H_r + E = maxwell_solution(x) + curlH = curlH_exact(x) + return omega * epsilon * E.imag + curlH.real + + @mfem.jit.vector(vdim=dim) + def rhs_func_i(x): + # J_r = ω ϵ E_i + ∇ × H_r + E = maxwell_solution(x) + curlH = curlH_exact(x) + return -omega * epsilon * E.real + curlH.imag + + self.rhs_func_r = rhs_func_r + self.rhs_func_i = rhs_func_i + + @mfem.jit.vector(vdim=dim) + def source_function(x): + center = np.zeros(dim)+0.5 + r = 0.0 + for i in range(dim): + r += (x[i] - center[i])**2 + + n = 5.0 * omega * np.sqrt(epsilon * mu) / pi + coeff = pow(n, 2) / pi + alpha = -pow(n, 2) * r + f = np.zeros(dim) + f[0] = -omega * coeff * exp(alpha)/omega + return f + + self.source_function = source_function + + +def run(meshfile='', + order=1, + delta_order=1, + prob=0, + sr=0, + pr=1, + epsilon=1.0, + mu=1.0, + rnum=1.0, + theta=0.0, + static_cond=False, + visualization=False): + + omega = 2.*pi*rnum + with_pml = False + exact_known = False + + if prob == 0: + exact_known = True + mesh_file = expanduser( + join(dirname(__file__), '..', '..', 'data', meshfile)) + + elif prob == 1: + meshfile = "meshes/fichera-waveguide.mesh" + omega = 5.0 + rnum = omega/(2.*pi) + mesh_file = expanduser( + join(dirname(__file__), meshfile)) + + elif prob == 2: + with_pml = True + mesh_file = expanduser( + join(dirname(__file__), '..', '..', 'data', meshfile)) + + else: + with_pml = True + meshfile = "meshes/scatter.mesh" + mesh_file = expanduser( + join(dirname(__file__), meshfile)) + + mesh = mfem.Mesh(mesh_file, 1, 1) + dim = mesh.Dimension() + assert dim > 1, "Dimension = 1 is not supported in this example" + + dimc = 3 if dim == 3 else 1 + for i in range(sr): + mesh.UniformRefinement() + mesh.EnsureNCMesh(False) + + pml = None + if with_pml: + length = mfem.doubleArray2D(dim, 2) + length.Assign(0.25) + pml = mdpg.CartesianPML(mesh, length) + pml.SetOmega(omega) + pml.SetEpsilonAndMu(epsilon, mu) + + pmesh = mfem.ParMesh(MPI.COMM_WORLD, mesh) + del mesh + + attr = mfem.intArray() + attrPML = mfem.intArray() + if with_pml: + pml.SetAttributes(pmesh, attr, attrPML) + + # Define spaces + TrialSpace = {"E_space": 0, + "H_space": 1, + "hatE_space": 2, + "hatH_space": 3, } + TestSpace = {"F_space": 0, + "G_space": 1, } + + # Vector L2 L2 space for E + E_fec = mfem.L2_FECollection(order-1, dim) + E_fes = mfem.ParFiniteElementSpace(pmesh, E_fec, dim) + + # Vector L2 L2 space for H + H_fec = mfem.L2_FECollection(order-1, dim) + H_fes = mfem.ParFiniteElementSpace(pmesh, H_fec, dimc) + + # H^-1/2 (curl) space for Ê + test_order = order + delta_order + + if dim == 3: + hatE_fec = mfem.ND_Trace_FECollection(order, dim) + hatH_fec = mfem.ND_Trace_FECollection(order, dim) + F_fec = mfem.ND_FECollection(test_order, dim) + else: + hatE_fec = mfem.RT_Trace_FECollection(order-1, dim) + hatH_fec = mfem.H1_Trace_FECollection(order, dim) + F_fec = mfem.H1_FECollection(test_order, dim) + + hatE_fes = mfem.ParFiniteElementSpace(pmesh, hatE_fec) + hatH_fes = mfem.ParFiniteElementSpace(pmesh, hatH_fec) + G_fec = mfem.ND_FECollection(test_order, dim) + + # trial_fes = mfem.ParFiniteElementSpaceArray() + # test_fec = mfem.FiniteElementCollectionArray() + trial_fes = [] + test_fec = [] + + trial_fes.append(E_fes) + trial_fes.append(H_fes) + trial_fes.append(hatE_fes) + trial_fes.append(hatH_fes) + test_fec.append(F_fec) + test_fec.append(G_fec) + + # Bilinear form coefficients + one = mfem.ConstantCoefficient(1.0) + eps2omeg2 = mfem.ConstantCoefficient(epsilon*epsilon*omega*omega) + mu2omeg2 = mfem.ConstantCoefficient(mu*mu*omega*omega) + muomeg = mfem.ConstantCoefficient(mu*omega) + negepsomeg = mfem.ConstantCoefficient(-epsilon*omega) + epsomeg = mfem.ConstantCoefficient(epsilon*omega) + negmuomeg = mfem.ConstantCoefficient(-mu*omega) + + # for the 2D case + rot_mat = mfem.DenseMatrix(np.array([[0, 1.], [-1, 0]])) + rot = mfem.MatrixConstantCoefficient(rot_mat) + epsrot = mfem.ScalarMatrixProductCoefficient(epsomeg, rot) + negepsrot = mfem.ScalarMatrixProductCoefficient(negepsomeg, rot) + + if with_pml: + epsomeg_cf = mfem.RestrictedCoefficient(epsomeg, attr) + negepsomeg_cf = mfem.RestrictedCoefficient(negepsomeg, attr) + eps2omeg2_cf = mfem.RestrictedCoefficient(eps2omeg2, attr) + muomeg_cf = mfem.RestrictedCoefficient(muomeg, attr) + negmuomeg_cf = mfem.RestrictedCoefficient(negmuomeg, attr) + mu2omeg2_cf = mfem.RestrictedCoefficient(mu2omeg2, attr) + epsrot_cf = mfem.MatrixRestrictedCoefficient(epsrot, attr) + negepsrot_cf = mfem.MatrixRestrictedCoefficient(negepsrot, attr) + else: + epsomeg_cf = epsomeg + negepsomeg_cf = negepsomeg + eps2omeg2_cf = eps2omeg2 + muomeg_cf = muomeg + negmuomeg_cf = negmuomeg + mu2omeg2_cf = mu2omeg2 + epsrot_cf = epsrot + negepsrot_cf = negepsrot + + detJ_r = mdpg.PmlCoefficient(mdpg.detJ_r_function, pml) + detJ_i = mdpg.PmlCoefficient(mdpg.detJ_i_function, pml) + abs_detJ_2 = mdpg.PmlCoefficient(mdpg.abs_detJ_2_function, pml) + detJ_Jt_J_inv_r = mdpg.PmlMatrixCoefficient( + dim, mdpg.detJ_Jt_J_inv_r_function, pml) + detJ_Jt_J_inv_i = mdpg.PmlMatrixCoefficient( + dim, mdpg.detJ_Jt_J_inv_i_function, pml) + abs_detJ_Jt_J_inv_2 = mdpg.PmlMatrixCoefficient( + dim, mdpg.abs_detJ_Jt_J_inv_2_function, pml) + negmuomeg_detJ_r = mfem.ProductCoefficient(negmuomeg, detJ_r) + negmuomeg_detJ_i = mfem.ProductCoefficient(negmuomeg, detJ_i) + muomeg_detJ_r = mfem.ProductCoefficient(muomeg, detJ_r) + mu2omeg2_detJ_2 = mfem.ProductCoefficient(mu2omeg2, abs_detJ_2) + epsomeg_detJ_Jt_J_inv_i = mfem.ScalarMatrixProductCoefficient( + epsomeg, detJ_Jt_J_inv_i) + epsomeg_detJ_Jt_J_inv_r = mfem.ScalarMatrixProductCoefficient( + epsomeg, detJ_Jt_J_inv_r) + negepsomeg_detJ_Jt_J_inv_r = mfem.ScalarMatrixProductCoefficient( + negepsomeg, detJ_Jt_J_inv_r) + muomeg_detJ_Jt_J_inv_r = mfem.ScalarMatrixProductCoefficient( + muomeg, detJ_Jt_J_inv_r) + negmuomeg_detJ_Jt_J_inv_i = mfem.ScalarMatrixProductCoefficient( + negmuomeg, detJ_Jt_J_inv_i) + negmuomeg_detJ_Jt_J_inv_r = mfem.ScalarMatrixProductCoefficient( + negmuomeg, detJ_Jt_J_inv_r) + mu2omeg2_detJ_Jt_J_inv_2 = mfem.ScalarMatrixProductCoefficient( + mu2omeg2, abs_detJ_Jt_J_inv_2) + eps2omeg2_detJ_Jt_J_inv_2 = mfem.ScalarMatrixProductCoefficient( + eps2omeg2, abs_detJ_Jt_J_inv_2) + negmuomeg_detJ_r_restr = mfem.RestrictedCoefficient( + negmuomeg_detJ_r, attrPML) + negmuomeg_detJ_i_restr = mfem.RestrictedCoefficient( + negmuomeg_detJ_i, attrPML) + muomeg_detJ_r_restr = mfem.RestrictedCoefficient(muomeg_detJ_r, attrPML) + mu2omeg2_detJ_2_restr = mfem.RestrictedCoefficient( + mu2omeg2_detJ_2, attrPML) + epsomeg_detJ_Jt_J_inv_i_restr = mfem.MatrixRestrictedCoefficient( + epsomeg_detJ_Jt_J_inv_i, attrPML) + epsomeg_detJ_Jt_J_inv_r_restr = mfem.MatrixRestrictedCoefficient( + epsomeg_detJ_Jt_J_inv_r, attrPML) + negepsomeg_detJ_Jt_J_inv_r_restr = mfem.MatrixRestrictedCoefficient( + negepsomeg_detJ_Jt_J_inv_r, attrPML) + muomeg_detJ_Jt_J_inv_r_restr = mfem.MatrixRestrictedCoefficient( + muomeg_detJ_Jt_J_inv_r, attrPML) + negmuomeg_detJ_Jt_J_inv_i_restr = mfem.MatrixRestrictedCoefficient( + negmuomeg_detJ_Jt_J_inv_i, attrPML) + negmuomeg_detJ_Jt_J_inv_r_restr = mfem.MatrixRestrictedCoefficient( + negmuomeg_detJ_Jt_J_inv_r, attrPML) + mu2omeg2_detJ_Jt_J_inv_2_restr = mfem.MatrixRestrictedCoefficient( + mu2omeg2_detJ_Jt_J_inv_2, attrPML) + eps2omeg2_detJ_Jt_J_inv_2_restr = mfem.MatrixRestrictedCoefficient( + eps2omeg2_detJ_Jt_J_inv_2, attrPML) + + if with_pml and dim == 2: + epsomeg_detJ_Jt_J_inv_i_rot = mfem.MatrixProductCoefficient( + epsomeg_detJ_Jt_J_inv_i, rot) + epsomeg_detJ_Jt_J_inv_r_rot = mfem.MatrixProductCoefficient( + epsomeg_detJ_Jt_J_inv_r, rot) + negepsomeg_detJ_Jt_J_inv_r_rot = mfem.MatrixProductCoefficient( + negepsomeg_detJ_Jt_J_inv_r, rot) + epsomeg_detJ_Jt_J_inv_i_rot_restr = mfem.MatrixRestrictedCoefficient(epsomeg_detJ_Jt_J_inv_i_rot, + attrPML) + epsomeg_detJ_Jt_J_inv_r_rot_restr = mfem.MatrixRestrictedCoefficient(epsomeg_detJ_Jt_J_inv_r_rot, + attrPML) + negepsomeg_detJ_Jt_J_inv_r_rot_restr = mfem.MatrixRestrictedCoefficient(negepsomeg_detJ_Jt_J_inv_r_rot, + attrPML) + + a = mdpg.ParComplexDPGWeakForm(trial_fes, test_fec) + a.StoreMatrices() # needed for AMR + + # (E,∇ × F) + a.AddTrialIntegrator(mfem.TransposeIntegrator(mfem.MixedCurlIntegrator(one)), + None, + TrialSpace["E_space"], + TestSpace["F_space"]) + # -i ω ϵ (E , G) = i (- ω ϵ E, G) + a.AddTrialIntegrator(None, + mfem.TransposeIntegrator( + mfem.VectorFEMassIntegrator(negepsomeg_cf)), + TrialSpace["E_space"], + TestSpace["G_space"]) + # (H,∇ × G) + a.AddTrialIntegrator(mfem.TransposeIntegrator(mfem.MixedCurlIntegrator(one)), + None, + TrialSpace["H_space"], + TestSpace["G_space"]) + # < n×Ĥ ,G> + a.AddTrialIntegrator(mfem.TangentTraceIntegrator(), None, + TrialSpace["hatH_space"], + TestSpace["G_space"]) + # test integrators + # (∇×G ,∇× δG) + a.AddTestIntegrator(mfem.CurlCurlIntegrator(one), None, + TestSpace["G_space"], + TestSpace["G_space"]) + # (G,δG) + a.AddTestIntegrator(mfem.VectorFEMassIntegrator(one), None, + TestSpace["G_space"], + TestSpace["G_space"]) + + if dim == 3: + # i ω μ (H, F) + a.AddTrialIntegrator(None, mfem.TransposeIntegrator( + mfem.VectorFEMassIntegrator(muomeg_cf)), + TrialSpace["H_space"], + TestSpace["F_space"]) + # < n×Ê,F> + a.AddTrialIntegrator(mfem.TangentTraceIntegrator(), None, + TrialSpace["hatE_space"], + TestSpace["F_space"]) + + # test integrators + # (∇×F,∇×δF) + a.AddTestIntegrator(mfem.CurlCurlIntegrator(one), None, + TestSpace["F_space"], + TestSpace["F_space"]) + # (F,δF) + a.AddTestIntegrator(mfem.VectorFEMassIntegrator(one), None, + TestSpace["F_space"], + TestSpace["F_space"]) + # μ^2 ω^2 (F,δF) + a.AddTestIntegrator(mfem.VectorFEMassIntegrator(mu2omeg2_cf), None, + TestSpace["F_space"], + TestSpace["F_space"]) + # -i ω μ (F,∇ × δG) = i (F, -ω μ ∇ × δ G) + a.AddTestIntegrator(None, mfem.MixedVectorWeakCurlIntegrator(negmuomeg_cf), + TestSpace["F_space"], + TestSpace["G_space"]) + # -i ω ϵ (∇ × F, δG) + a.AddTestIntegrator(None, mfem.MixedVectorCurlIntegrator(negepsomeg_cf), + TestSpace["F_space"], + TestSpace["G_space"]) + # i ω μ (∇ × G,δF) + a.AddTestIntegrator(None, mfem.MixedVectorCurlIntegrator(muomeg_cf), + TestSpace["G_space"], + TestSpace["F_space"]) + # i ω ϵ (G, ∇ × δF ) + a.AddTestIntegrator(None, mfem.MixedVectorWeakCurlIntegrator(epsomeg_cf), + TestSpace["G_space"], + TestSpace["F_space"]) + # ϵ^2 ω^2 (G,δG) + a.AddTestIntegrator(mfem.VectorFEMassIntegrator(eps2omeg2_cf), None, + TestSpace["G_space"], + TestSpace["G_space"]) + else: + # i ω μ (H, F) + a.AddTrialIntegrator(None, mfem.MixedScalarMassIntegrator(muomeg_cf), + TrialSpace["H_space"], + TestSpace["F_space"]) + # < n×Ê,F> + a.AddTrialIntegrator(mfem.TraceIntegrator(), None, + TrialSpace["hatE_space"], + TestSpace["F_space"]) + # test integrators + # (∇F,∇δF) + a.AddTestIntegrator(mfem.DiffusionIntegrator(one), None, + TestSpace["F_space"], + TestSpace["F_space"]) + # (F,δF) + a.AddTestIntegrator(mfem.MassIntegrator(one), None, + TestSpace["F_space"], + TestSpace["F_space"]) + # μ^2 ω^2 (F,δF) + a.AddTestIntegrator(mfem.MassIntegrator(mu2omeg2_cf), None, + TestSpace["F_space"], + TestSpace["F_space"]) + # -i ω μ (F,∇ × δG) = i (F, -ω μ ∇ × δ G) + a.AddTestIntegrator(None, + mfem.TransposeIntegrator( + mfem.MixedCurlIntegrator(negmuomeg_cf)), + TestSpace["F_space"], + TestSpace["G_space"]) + # -i ω ϵ (∇ × F, δG) = i (- ω ϵ A ∇ F,δG), A = [0 1; -1; 0] + a.AddTestIntegrator(None, mfem.MixedVectorGradientIntegrator(negepsrot_cf), + TestSpace["F_space"], + TestSpace["G_space"]) + # i ω μ (∇ × G,δF) = i (ω μ ∇ × G, δF ) + a.AddTestIntegrator(None, mfem.MixedCurlIntegrator(muomeg_cf), + TestSpace["G_space"], + TestSpace["F_space"]) + # i ω ϵ (G, ∇ × δF ) = i (ω ϵ G, A ∇ δF) = i ( G , ω ϵ A ∇ δF) + a.AddTestIntegrator(None, + mfem.TransposeIntegrator( + mfem.MixedVectorGradientIntegrator(epsrot_cf)), + TestSpace["G_space"], + TestSpace["F_space"]) + # ϵ^2 ω^2 (G, δG) + a.AddTestIntegrator(mfem.VectorFEMassIntegrator(eps2omeg2_cf), None, + TestSpace["G_space"], + TestSpace["G_space"]) + if with_pml: + # trial integrators + # -i ω ϵ (β E , G) = -i ω ϵ ((β_re + i β_im) E, G) + # = (ω ϵ β_im E, G) + i (- ω ϵ β_re E, G) + a.AddTrialIntegrator(mfem.TransposeIntegrator(mfem.VectorFEMassIntegrator( + epsomeg_detJ_Jt_J_inv_i_restr)), + mfem.TransposeIntegrator(mfem.VectorFEMassIntegrator( + negepsomeg_detJ_Jt_J_inv_r_restr)), + TrialSpace["E_space"], TestSpace["G_space"]) + if dim == 3: + # trial integrators + # i ω μ (α^-1 H, F) = i ω μ ( (α^-1_re + i α^-1_im) H, F) + # = (- ω μ α^-1_im, H,F) + i *(ω μ α^-1_re H, F) + a.AddTrialIntegrator( + mfem.TransposeIntegrator(mfem.VectorFEMassIntegrator( + negmuomeg_detJ_Jt_J_inv_i_restr)), + mfem.TransposeIntegrator(mfem.VectorFEMassIntegrator( + muomeg_detJ_Jt_J_inv_r_restr)), + TrialSpace["H_space"], TestSpace["F_space"]) + # test integrators + # μ^2 ω^2 (|α|^-2 F,δF) + a.AddTestIntegrator( + mfem.VectorFEMassIntegrator( + mu2omeg2_detJ_Jt_J_inv_2_restr), None, + TestSpace["F_space"], TestSpace["F_space"]) + # -i ω μ (α^-* F,∇ × δG) = i (F, - ω μ α^-1 ∇ × δ G) + # = i (F, - ω μ (α^-1_re + i α^-1_im) ∇ × δ G) + # = (F, - ω μ α^-1_im ∇ × δ G) + i (F, - ω μ α^-1_re ∇×δG) + a.AddTestIntegrator(mfem.MixedVectorWeakCurlIntegrator( + negmuomeg_detJ_Jt_J_inv_i_restr), + mfem.MixedVectorWeakCurlIntegrator( + negmuomeg_detJ_Jt_J_inv_r_restr), + TestSpace["F_space"], TestSpace["G_space"]) + # -i ω ϵ (β ∇ × F, δG) = -i ω ϵ ((β_re + i β_im) ∇ × F, δG) + # = (ω ϵ β_im ∇ × F, δG) + i (- ω ϵ β_re ∇ × F, δG) + a.AddTestIntegrator(mfem.MixedVectorCurlIntegrator( + epsomeg_detJ_Jt_J_inv_i_restr), + mfem.MixedVectorCurlIntegrator( + negepsomeg_detJ_Jt_J_inv_r_restr), + TestSpace["F_space"], TestSpace["G_space"]) + # i ω μ (α^-1 ∇ × G,δF) = i ω μ ((α^-1_re + i α^-1_im) ∇ × G,δF) + # = (- ω μ α^-1_im ∇ × G,δF) + i (ω μ α^-1_re ∇ × G,δF) + a.AddTestIntegrator(mfem.MixedVectorCurlIntegrator( + negmuomeg_detJ_Jt_J_inv_i_restr), + mfem.MixedVectorCurlIntegrator( + muomeg_detJ_Jt_J_inv_r_restr), + TestSpace["G_space"], TestSpace["F_space"]) + # i ω ϵ (β^* G, ∇×δF) = i ω ϵ ( (β_re - i β_im) G, ∇×δF) + # = (ω ϵ β_im G, ∇×δF) + i ( ω ϵ β_re G, ∇×δF) + a.AddTestIntegrator(mfem.MixedVectorWeakCurlIntegrator( + epsomeg_detJ_Jt_J_inv_i_restr), + mfem.MixedVectorWeakCurlIntegrator( + epsomeg_detJ_Jt_J_inv_r_restr), + TestSpace["G_space"], TestSpace["F_space"]) + # ϵ^2 ω^2 (|β|^2 G,δG) + a.AddTestIntegrator(mfem.VectorFEMassIntegrator( + eps2omeg2_detJ_Jt_J_inv_2_restr), None, + TestSpace["G_space"], TestSpace["G_space"]) + + else: + # trial integrators + # i ω μ (α^-1 H, F) = i ω μ ( (α^-1_re + i α^-1_im) H, F) + # = (- ω μ α^-1_im, H,F) + i *(ω μ α^-1_re H, F) + a.AddTrialIntegrator( + mfem.MixedScalarMassIntegrator(negmuomeg_detJ_i_restr), + mfem.MixedScalarMassIntegrator(muomeg_detJ_r_restr), + TrialSpace["H_space"], TestSpace["F_space"]) + # test integrators + # μ^2 ω^2 (|α|^-2 F,δF) + a.AddTestIntegrator(mfem.MassIntegrator(mu2omeg2_detJ_2_restr), None, + TestSpace["F_space"], TestSpace["F_space"]) + # -i ω μ (α^-* F,∇ × δG) = (F, ω μ α^-1 ∇ × δ G) + # =(F, - ω μ α^-1_im ∇ × δ G) + i (F, - ω μ α^-1_re ∇×δG) + a.AddTestIntegrator( + mfem.TransposeIntegrator( + mfem.MixedCurlIntegrator(negmuomeg_detJ_i_restr)), + mfem.TransposeIntegrator( + mfem.MixedCurlIntegrator(negmuomeg_detJ_r_restr)), + TestSpace["F_space"], TestSpace["G_space"]) + # -i ω ϵ (β ∇ × F, δG) = i (- ω ϵ β A ∇ F,δG), A = [0 1; -1; 0] + # = (ω ϵ β_im A ∇ F, δG) + i (- ω ϵ β_re A ∇ F, δG) + a.AddTestIntegrator(mfem.MixedVectorGradientIntegrator( + epsomeg_detJ_Jt_J_inv_i_rot_restr), + mfem.MixedVectorGradientIntegrator( + negepsomeg_detJ_Jt_J_inv_r_rot_restr), + TestSpace["F_space"], TestSpace["G_space"]) + # i ω μ (α^-1 ∇ × G,δF) = i (ω μ α^-1 ∇ × G, δF ) + # = (- ω μ α^-1_im ∇ × G,δF) + i (ω μ α^-1_re ∇ × G,δF) + a.AddTestIntegrator(mfem.MixedCurlIntegrator(negmuomeg_detJ_i_restr), + mfem.MixedCurlIntegrator(muomeg_detJ_r_restr), + TestSpace["G_space"], TestSpace["F_space"]) + # i ω ϵ (β^* G, ∇ × δF ) = i ( G , ω ϵ β A ∇ δF) + # = ( G , ω ϵ β_im A ∇ δF) + i ( G , ω ϵ β_re A ∇ δF) + a.AddTestIntegrator( + mfem.TransposeIntegrator(mfem.MixedVectorGradientIntegrator( + epsomeg_detJ_Jt_J_inv_i_rot_restr)), + mfem.TransposeIntegrator(mfem.MixedVectorGradientIntegrator( + epsomeg_detJ_Jt_J_inv_r_rot_restr)), + TestSpace["G_space"], TestSpace["F_space"]) + # ϵ^2 ω^2 (|β|^2 G,δG) + a.AddTestIntegrator(mfem.VectorFEMassIntegrator( + eps2omeg2_detJ_Jt_J_inv_2_restr), None, + TestSpace["G_space"], TestSpace["G_space"]) + + # RHS + fncs = Functions(dim, dimc, omega, epsilon, mu, prob) + if prob == 0: + f_rhs_r = fncs.rhs_func_r + f_rhs_i = fncs.rhs_func_i + a.AddDomainLFIntegrator(mfem.VectorFEDomainLFIntegrator(f_rhs_r), + mfem.VectorFEDomainLFIntegrator(f_rhs_i), + TestSpace["G_space"]) + elif prob == 2: + f_source = fncs.source_function + a.AddDomainLFIntegrator(mfem.VectorFEDomainLFIntegrator(f_source), + None, + TestSpace["G_space"]) + + hatEex_r = fncs.hatE_exact_r + hatEex_i = fncs.hatE_exact_i + + if myid == 0: + txt = "\n Ref |" + " Dofs |" + " ω |" + if exact_known: + txt = txt + " L2 Error |" + " Rate |" + + txt = txt + " Residual |" + " Rate |" + " PCG it |" + print(txt) + + if exact_known: + print("-"*82) + else: + print("-"*60) + + res0 = 0. + err0 = 0. + dof0 = 0 + + elements_to_refine = mfem.intArray() + + # visualizaiton socket + E_out_r = None + H_out_r = None + + E_r = mfem.ParGridFunction() + E_i = mfem.ParGridFunction() + H_r = mfem.ParGridFunction() + H_i = mfem.ParGridFunction() + + if static_cond: + a.EnableStaticCondensation() + + for it in range(pr+1): + a.Assemble() + + ess_tdof_list = mfem.intArray() + ess_bdr = mfem.intArray() + + if pmesh.bdr_attributes.Size() != 0: + ess_bdr.SetSize(pmesh.bdr_attributes.Max()) + ess_bdr.Assign(1) + hatE_fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list) + if with_pml: + ess_bdr.Assign(0) + ess_bdr[1] = 1 + + # shift the ess_tdofs + for j in range(ess_tdof_list.Size()): + ess_tdof_list[j] += E_fes.GetTrueVSize() + H_fes.GetTrueVSize() + + offsets = mfem.intArray(5) + offsets[0] = 0 + offsets[1] = E_fes.GetVSize() + offsets[2] = H_fes.GetVSize() + offsets[3] = hatE_fes.GetVSize() + offsets[4] = hatH_fes.GetVSize() + + offsets.PartialSum() + offsetsl = offsets.ToList() + + x = mfem.Vector([0.]*(2*offsetsl[-1])) + if prob != 2: + hatE_gf_r = mfem.ParGridFunction(hatE_fes, x, offsetsl[2]) + hatE_gf_i = mfem.ParGridFunction( + hatE_fes, x, offsetsl[-1] + offsetsl[2]) + if dim == 3: + hatE_gf_r.ProjectBdrCoefficientTangent(hatEex_r, ess_bdr) + hatE_gf_i.ProjectBdrCoefficientTangent(hatEex_i, ess_bdr) + else: + hatE_gf_r.ProjectBdrCoefficientNormal(hatEex_r, ess_bdr) + hatE_gf_i.ProjectBdrCoefficientNormal(hatEex_i, ess_bdr) + Ah = mfem.OperatorPtr() + X = mfem.Vector() + B = mfem.Vector() + + a.FormLinearSystem(ess_tdof_list, x, Ah, X, B) + + Ahc = Ah.AsComplexOperator() + + BlockA_r = mfem.Opr2BlockOpr(Ahc.real()) + BlockA_i = mfem.Opr2BlockOpr(Ahc.imag()) + num_blocks = BlockA_r.NumRowBlocks() + + # this is to debug matrix + # for i in range(num_blocks): + # for j in range(num_blocks): + # (mfem.Opr2HypreParMatrix(BlockA_r.GetBlock(i, j))).Print(str(i)+"_"+str(j)+"r") + # (mfem.Opr2HypreParMatrix(BlockA_i.GetBlock(i, j))).Print(str(i)+"_"+str(j)+"i") + + tdof_offsets = mfem.intArray(2*num_blocks+1) + + tdof_offsets[0] = 0 + skip = 0 if static_cond else 2 + k = 2 if static_cond else 0 + for i in range(num_blocks): + tdof_offsets[i+1] = trial_fes[i+k].GetTrueVSize() + tdof_offsets[num_blocks+i+1] = trial_fes[i+k].GetTrueVSize() + tdof_offsets.PartialSum() + + blockA = mfem.BlockOperator(tdof_offsets) + for i in range(num_blocks): + for j in range(num_blocks): + blockA.SetBlock(i, j, BlockA_r.GetBlock(i, j)) + blockA.SetBlock(i, j+num_blocks, BlockA_i.GetBlock(i, j), -1.0) + blockA.SetBlock(i+num_blocks, j+num_blocks, + BlockA_r.GetBlock(i, j)) + blockA.SetBlock(i+num_blocks, j, BlockA_i.GetBlock(i, j)) + X.Assign(0.0) + M = mfem.BlockDiagonalPreconditioner(tdof_offsets) + + if not static_cond: + E_mat = mfem.Opr2HypreParMatrix(BlockA_r.GetBlock(0, 0)) + H_mat = mfem.Opr2HypreParMatrix(BlockA_r.GetBlock(1, 1)) + + solver_E = mfem.HypreBoomerAMG(E_mat) + solver_E.SetPrintLevel(0) + solver_E.SetSystemsOptions(dim) + solver_H = mfem.HypreBoomerAMG(H_mat) + solver_H.SetPrintLevel(0) + solver_H.SetSystemsOptions(dim) + M.SetDiagonalBlock(0, solver_E) + M.SetDiagonalBlock(1, solver_H) + M.SetDiagonalBlock(num_blocks, solver_E) + M.SetDiagonalBlock(num_blocks+1, solver_H) + + hatE_mat = mfem.Opr2HypreParMatrix(BlockA_r.GetBlock(skip, skip)) + solver_hatE = mfem.HypreAMS(hatE_mat, hatE_fes) + solver_hatE.SetPrintLevel(0) + + hatH_mat = mfem.Opr2HypreParMatrix(BlockA_r.GetBlock(skip+1, skip+1)) + if dim == 2: + solver_hatH = mfem.HypreBoomerAMG(hatH_mat) + solver_hatH.SetPrintLevel(0) + else: + solver_hatH = mfem.HypreAMS(hatH_mat, hatH_fes) + solver_hatH.SetPrintLevel(0) + + M.SetDiagonalBlock(skip, solver_hatE) + M.SetDiagonalBlock(skip+1, solver_hatH) + M.SetDiagonalBlock(skip+num_blocks, solver_hatE) + M.SetDiagonalBlock(skip+num_blocks+1, solver_hatH) + + cg = mfem.CGSolver(MPI.COMM_WORLD) + cg.SetRelTol(1e-6) + cg.SetMaxIter(10000) + cg.SetPrintLevel(0) + cg.SetPreconditioner(M) + cg.SetOperator(blockA) + cg.Mult(B, X) + + # for i in range(num_blocks): + # delete &M.GetDiagonalBlock(i); + + num_iter = cg.GetNumIterations() + + a.RecoverFEMSolution(X, x) + + residuals = a.ComputeResidual(x) + + residual = residuals.Norml2() + maxresidual = residuals.Max() + globalresidual = residual * residual + + maxresidual = MPI.COMM_WORLD.allreduce(maxresidual, op=MPI.MAX) + globalresidual = MPI.COMM_WORLD.allreduce(globalresidual, op=MPI.SUM) + globalresidual = np.sqrt(globalresidual) + + E_r.MakeRef(E_fes, x, 0) + E_i.MakeRef(E_fes, x, offsetsl[-1]) + + H_r.MakeRef(H_fes, x, offsetsl[1]) + H_i.MakeRef(H_fes, x, offsetsl[-1]+offsetsl[1]) + + dofs = 0 + for i in range(len(trial_fes)): + dofs += trial_fes[i].GlobalTrueVSize() + + if exact_known: + E_ex_r = fncs.E_exact_r + E_ex_i = fncs.E_exact_i + H_ex_r = fncs.H_exact_r + H_ex_i = fncs.H_exact_i + E_err_r = E_r.ComputeL2Error(E_ex_r) + E_err_i = E_r.ComputeL2Error(E_ex_i) + H_err_r = H_r.ComputeL2Error(E_ex_r) + H_err_i = H_r.ComputeL2Error(E_ex_i) + L2Error = np.sqrt(E_err_r*E_err_r + E_err_i*E_err_i + + H_err_r*H_err_r + H_err_i*H_err_i) + rate_err = 0 if it == 0 else dim * \ + np.log(err0/L2Error)/np.log(dof0/dofs) + err0 = L2Error + + rate_res = 0.0 if it == 0 else dim * \ + np.log(res0/globalresidual)/np.log(dof0/dofs) + + res0 = globalresidual + dof0 = dofs + + if myid == 0: + txt = ("{:5d}".format(it) + " | " + + "{:10d}".format(dof0) + " | " + + " {:.1f}".format(2.0*rnum) + " π | ") + if exact_known: + txt = txt + (" {:.3e}".format(err0) + " | " + + " {:.2f} ".format(rate_err) + " | ") + + txt = txt + (" {:.3e}".format(res0) + " | " + + " {:.2f} ".format(rate_res) + " | " + + "{:6d}".format(num_iter) + " | ") + print(txt) + + if visualization: + keys = "jRcml\n" if it == 0 and dim == 2 else "" + E_out_r = VisualizeField(E_out_r, "localhost", visport, E_r, + "Numerical Electric field (real part)", 0, 0, 500, 500, keys) + H_out_r = VisualizeField(H_out_r, "localhost", visport, H_r, + "Numerical Magnetic field (real part)", 0, 0, 500, 500, keys) + + if it == pr: + break + + if theta > 0.0: + elements_to_refine.SetSize(0) + for iel in range(pmesh.GetNE()): + if residuals[iel] > theta * maxresidual: + elements_to_refine.Append(iel) + pmesh.GeneralRefinement(elements_to_refine, 1, 1) + else: + pmesh.UniformRefinement() + + if with_pml: + pml.SetAttributes(pmesh) + for i in range(len(trial_fes)): + trial_fes[i].Update(False) + a.Update() + + +if __name__ == "__main__": + from mfem.common.arg_parser import ArgParser + + parser = ArgParser(description='Ex40 (Eikonal queation)') + parser.add_argument('-m', '--mesh', + default="inline-quad.mesh", + action='store', type=str, + help='Mesh file to use.') + parser.add_argument('-o', '--order', + action='store', default=1, type=int, + help="Finite element order (polynomial degree).") + + parser.add_argument("-rnum", "--number-of-wavelengths", + action='store', default=1.0, type=float, + help="Number of wavelengths") + parser.add_argument("-mu", "--permeability", + action='store', default=1.0, type=float, + help="Permeability of free space (or 1/(spring constant)).") + parser.add_argument("-eps", "--permittivity", + action='store', default=1.0, type=float, + help="Permittivity of free space (or mass constant).") + parser.add_argument("-prob", "--problem", + action='store', default=0, type=int, + help="\n".join(("Problem case" + " 0: plane wave, 1: Fichera 'oven', " + " 2: Generic PML problem with point source given as a load " + " 3: Scattering of a plane wave, " + " 4: Point source given on the boundary"))) + parser.add_argument("-do", "--delta-order", + action='store', default=1, type=int, + help="Order enrichment for DPG test space.") + parser.add_argument("-theta", "--theta", + action='store', default=0.0, type=float, + help="Theta parameter for AMR") + parser.add_argument("-sref", "--serial-ref", + action='store', default=0, type=int, + help="Number of parallel refinements.") + parser.add_argument("-pref", "--parallel-ref", + action='store', default=1, type=int, + help="Number of parallel refinements.") + parser.add_argument("-sc", "--static-condensation", + action='store_true', default=False, + help="Enable static condensation.") + parser.add_argument('-no-vis', '--no-visualization', + action='store_true', + default=False, + help='Disable or disable GLVis visualization') + + args = parser.parse_args() + if myid == 0: + parser.print_options(args) + + visualization = not args.no_visualization + + run(meshfile=args.mesh, + order=args.order, + prob=args.problem, + sr=args.serial_ref, + pr=args.parallel_ref, + epsilon=args.permittivity, + mu=args.permeability, + delta_order=args.delta_order, + rnum=args.number_of_wavelengths, + theta=args.theta, + static_cond=args.static_condensation, + visualization=visualization)