diff --git a/.gitignore b/.gitignore index eb4c6658c..8165839ed 100644 --- a/.gitignore +++ b/.gitignore @@ -135,7 +135,10 @@ AppDir/ 3rdparty/mmg3d/mmg3d4/ 3rdparty/mmg3d/tag-tar-158547 3rdparty/mshmet/ff-flags -3rdparty/mshmet/mshmet.2012.04.25/ +3rdparty/mshmet/src/FAIRE +3rdparty/mshmet/src/FAIT +3rdparty/mshmet/src/sourceslib/compil.date +3rdparty/mshmet/src/tag-tar 3rdparty/mumps-seq/FAIRE-4.10.0 3rdparty/mumps-seq/Makefile.inc 3rdparty/mumps/FAIRE-4.10.0 @@ -270,6 +273,7 @@ plugin/seq/ff-c++ plugin/seq/ff-get-dep plugin/seq/ff-pkg-download plugin/seq/ffmaster +plugin/seq/ffmaster.dSYM plugin/seq/include/ plugin/seq/load.link *speedtest.out diff --git a/3rdparty/Makefile.am b/3rdparty/Makefile.am index 19207604d..fddea1144 100755 --- a/3rdparty/Makefile.am +++ b/3rdparty/Makefile.am @@ -76,7 +76,6 @@ ipopt/patch-IpBlas \ ff-petsc/Makefile \ ff-petsc/Makefile-PETSc.inc \ ff-petsc/Makefile.complex \ -ff-petsc/petsc-3.11.2.patch \ getall # FFCS: See [[file:../../../configure.ac::tools_problems_all_platforms]] for reasons why some tools may be deactivated @@ -118,7 +117,7 @@ install-other: tag-install-other tag-install-other: if test -n "$(TOOL_COMPILE_hpddm)"; then $(MAKE) install-hpddm; else true; fi if test -n "$(TOOL_COMPILE_htool)"; then $(MAKE) install-htool; else true; fi - if test -n "$(TOOL_COMPILE_bem)"; then $(MAKE) install-boost install-bemtool; else true; fi + if test -n "$(TOOL_COMPILE_bemtool)"; then $(MAKE) install-boost install-bemtool; else true; fi touch tag-install-other WHERE-OTHER: lib/WHERE.hpddm lib/WHERE.htool lib/WHERE.bemtool lib/WHERE.boost @@ -169,12 +168,12 @@ lib/WHERE.htool: ./getall # begin BEMTOOL install-bemtool: ./getall - if test -n "$(TOOL_COMPILE_bem)"; then ./getall -o bemtool -a; $(MAKE) include/BemTool/done.tag lib/WHERE.bemtool; else true;fi + if test -n "$(TOOL_COMPILE_bemtool)"; then ./getall -o bemtool -a; $(MAKE) include/BemTool/done.tag lib/WHERE.bemtool; else true;fi reinstall-bemtool:./getall -rm lib/WHERE.bemtool -rm ./pkg/bemtool.zip -rm -rf include/BemTool - -test -n "$(TOOL_COMPILE_bem)" && $(MAKE) install-bemtool || true + -test -n "$(TOOL_COMPILE_bemtool)" && $(MAKE) install-bemtool || true include/BemTool/done.tag: ./pkg/bemtool.zip -if test -f ./pkg/bemtool.zip; then \ cd include; rm -rf BemTool BemTool-*; unzip -q ../pkg/bemtool.zip; mv BemTool-* BemTool; touch BemTool/done.tag; \ @@ -188,12 +187,12 @@ lib/WHERE.bemtool: ./getall # begin BOOST install-boost: ./getall - if test -n "$(TOOL_COMPILE_bem)"; then ./getall -o Boost -a; $(MAKE) boost/done.tag lib/WHERE.boost; else true; fi + if test -n "$(TOOL_COMPILE_bemtool)"; then ./getall -o Boost -a; $(MAKE) boost/done.tag lib/WHERE.boost; else true; fi reinstall-boost: ./getall -rm lib/WHERE.boost -rm ./pkg/boost_for_bemtool.tar.gz -rm -rf boost - -test -n "$(TOOL_COMPILE_bem)" && $(MAKE) install-boost || true + -test -n "$(TOOL_COMPILE_bemtool)" && $(MAKE) install-boost || true boost/done.tag: ./pkg/boost_for_bemtool.tar.gz if test -f ./pkg/boost_for_bemtool.tar.gz; then \ rm -rf boost boost_*; tar xzf pkg/boost_for_bemtool.tar.gz; mv boost_for_bemtool boost; \ diff --git a/3rdparty/ff-petsc/Makefile b/3rdparty/ff-petsc/Makefile index 0bf0c5a98..3c1d9f326 100644 --- a/3rdparty/ff-petsc/Makefile +++ b/3rdparty/ff-petsc/Makefile @@ -203,7 +203,7 @@ $(SRCDIR)/tag-make-real:$(SRCDIR)/tag-conf-real $(SRCDIR)/tag-install-real :$(SRCDIR)/tag-make-real cd $(SRCDIR) && $(MAKE) PETSC_DIR=$(PETSC_DIR) PETSC_ARCH=fr install-lib - cd $(SRCDIR) && $(MAKE) PETSC_DIR=$(PETSC_DIR) PETSC_ARCH=fr check + cd $(SRCDIR) && export HWLOC_HIDE_ERRORS=2 && $(MAKE) PETSC_DIR=$(PETSC_DIR) PETSC_ARCH=fr check test -e $(DIR_INSTALL_REAL)/include/petsc.h test -e $(DIR_INSTALL_REAL)/lib/petsc/conf/petscvariables touch $@ diff --git a/3rdparty/ff-petsc/petsc-cmake-seq.patch b/3rdparty/ff-petsc/petsc-cmake-seq.patch index 60b790f1f..ed4a878a1 100644 --- a/3rdparty/ff-petsc/petsc-cmake-seq.patch +++ b/3rdparty/ff-petsc/petsc-cmake-seq.patch @@ -66,9 +66,9 @@ diff --color -Naur petsc-3.19.0-p1/config/BuildSystem/config/packages/PTScotch.p @@ -9,6 +9,8 @@ self.download = ['git://https://gitlab.inria.fr/scotch/scotch.git', 'https://gitlab.inria.fr/scotch/scotch/-/archive/'+self.gitcommit+'/scotch-'+self.gitcommit+'.tar.gz', - 'http://ftp.mcs.anl.gov/pub/petsc/externalpackages/scotch-'+self.gitcommit+'.tar.gz'] + 'https://web.cels.anl.gov/projects/petsc/download/externalpackages/scotch-'+self.gitcommit+'.tar.gz'] + self.download_macos = ['https://gitlab.inria.fr/scotch/scotch/-/archive/v6.1.2/scotch-v6.1.2.tar.gz', -+ 'http://ftp.mcs.anl.gov/pub/petsc/externalpackages/scotch-v6.1.2.tar.gz'] ++ 'http://web.cels.anl.gov/projects/petsc/download/externalpackages/scotch-v6.1.2.tar.gz'] self.downloaddirnames = ['scotch','petsc-pkg-scotch'] self.liblist = [['libptesmumps.a','libptscotchparmetisv3.a','libptscotch.a','libptscotcherr.a','libesmumps.a','libscotch.a','libscotcherr.a'],['libptesmumps.a','libptscotchparmetis.a','libptscotch.a','libptscotcherr.a','libesmumps.a','libscotch.a','libscotcherr.a'], ['libptesmumps.a','libptscotchparmetis.a','libptscotch.a','libptscotcherr.a','libesmumps.a','libscotch.a','libscotcherr.a']] diff --git a/3rdparty/getall b/3rdparty/getall index 71ae6d902..85a0ce4ee 100755 --- a/3rdparty/getall +++ b/3rdparty/getall @@ -171,10 +171,10 @@ download('parmmg','https://github.com/prj-/ParMmg/archive/9cb2f7a22ef590d196d028 'https://github.com/prj-', 'parmmg.zip', '32035de70272121137c0869e8c49e4bd'); -download('PETSc','https://www.mcs.anl.gov/petsc/mirror/release-snapshots/petsc-3.19.2.tar.gz', - 'https://www.mcs.anl.gov/petsc/mirror/release-snapshots/', - 'petsc-3.19.2.tar.gz', - '9d53f784ca2ab37bf25b7991ceefa1d9'); +download('PETSc','https://web.cels.anl.gov/projects/petsc/download/release-snapshots/petsc-3.20.2.tar.gz', + 'https://web.cels.anl.gov/projects/petsc/download/release-snapshots/', + 'petsc-3.20.2.tar.gz', + '1e170a5f096433ca21aa643c80c749eb'); download('htool','https://github.com/htool-ddm/htool/archive/946875d79d0036afb4dc2c0c13c165a607d830df.zip', 'https://github.com/htool-ddm/', @@ -187,10 +187,10 @@ download('hpddm','https://codeload.github.com/hpddm/hpddm/zip/7113b9a6b77fceee3f 'hpddm.zip', '6910b7b974f0b60d9c247c666e7f3862'); -download('bemtool','https://github.com/PierreMarchand20/BemTool/archive/435b9d8ea86228b9ec4b526fc915aa35a7a8bec6.zip', +download('bemtool','https://github.com/PierreMarchand20/BemTool/archive/629c44513698405b58c50650cba69419474062ad.zip', 'https://github.com/PierreMarchand20/BemTool', 'bemtool.zip', - '3f9c95580efeb4d2b8f1fa977f7fd395'); + '869832f5cbec4dfb2c16e2d94bad0b7d'); download('Boost','https://www.ljll.math.upmc.fr/~tournier/boost_for_bemtool.tar.gz', 'https://www.boost.org', diff --git a/AutoGeneratedFile.tar.gz b/AutoGeneratedFile.tar.gz index 402c3d8ad..d0462e0d9 100644 Binary files a/AutoGeneratedFile.tar.gz and b/AutoGeneratedFile.tar.gz differ diff --git a/CHANGELOG.md b/CHANGELOG.md index 0654f6bc3..277ec3101 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,43 @@ # Changelog All notable changes to this project will be documented in this file. +## [4.14] +### Added +- Finite element BDM2 and BDM2ortho in test, Bug in BDM2ortho corrected the 4 sept 2014 in version: v4.13-130-g1af52457 +- Conversion of matrix or transpose of matrix in `int[int][int]` array to get the structure of sparse matrix. + see tutorial/sparse-matrix.edp example at end + ``` + matrix A = va(Ph,Vh); + int[int] a = A, at= A'; + ```` +- a meshL finite function can be see as real function with 1, or 2 parameters + ``` + meshL ThL = segment(10); fespace VhL(ThL,P1); VhL u= x; + cout << u(0.5) << endl; + cout << u(0.5,0) << endl; + ``` +- Exemple to code convolution of 2 function with one with a small support + too be not to expanxive + see tutorial/Convolution-Sample.edp example +- Support for dense blocks in PETSc matrices +- GenEO for saddle-point examples with PCHPDDM in PETSc +- Distributed ParaView output on `meshS` +- Interface to `mmg2d` for two-dimensional `mesh` +- Support for Mmg parameters `localParameter`, `-nosizreq`, -hgradreq` + +### Changed +- PETSc 3.20.2 + +### Deprecated +- + +### Removed +- + +### Fixed +- bug in P3pnc3d in vectorial case (thank to loic.balaziatchynillama@cea.fr ) +- in segment(10,region=1,label=ll); region is now used.. + ## [4.13] ### Added @@ -48,7 +85,7 @@ All notable changes to this project will be documented in this file. - add New Finite element 2d on mesh : RT0dc (discontinios RT0 ) in plugin Element_Mixte see example plugin/RT0dc.edp - and P1nc (Crouziex-Raviat) + bulle : name P1bnc in plugin Element_P1ncdc + and P1nc (Crouziex-Raviart) + bulle : name P1bnc in plugin Element_P1ncdc and P1nc totaly discontinous + bulle ; name P1bdcnc in plugin Element_P1ncdc see example plugin/example testp1dcnc.edp for akram.beni-hamad@inria.fr @@ -604,7 +641,7 @@ All notable changes to this project will be documented in this file. ### Changed - The main distribution is now on Github -[Unreleased]: https://github.com/FreeFem/FreeFem-sources/compare/v4.12..develop +[4.14]: https://github.com/FreeFem/FreeFem-sources/compare/v4.13..v4.14 [4.13]: https://github.com/FreeFem/FreeFem-sources/compare/v4.12..v4.13 [4.12]: https://github.com/FreeFem/FreeFem-sources/compare/v4.11..v4.12 [4.11]: https://github.com/FreeFem/FreeFem-sources/compare/v4.10..v4.11 diff --git a/bin/change-dylib-gfortran b/bin/change-dylib-gfortran new file mode 100755 index 000000000..e4670027f --- /dev/null +++ b/bin/change-dylib-gfortran @@ -0,0 +1,49 @@ +#!/bin/bash +## change path gfortran lib in freefem install +## usage : change-dylibgfortran oldlibdir newlibdir libnam1 libnam2 list of dylyb file pathname to change ... +## F. hecht sep. 2023 on MoOS Ventura , Monterey +## remark the file $oldlibdir/$libnam1 , ... must exist ... +olddir=$1 +shift +newdir=$1 +shift +lib1=$1 +shift +lib2=$1 +shift +lib3=$1 +shift +fi + +if test -f "$olddir/$lib1" -a -f "$olddir/$lib2" -a -f "$olddir/$lib3" \ + -a -f "$newdir/$lib1" -a -f "$newdir/$lib2" -a -f "$newdir/$lib3" ; then +for ff in $@; do +ch=0 +for lib in $lib1 $lib2 $lib3; do + old=$olddir/$lib + new=$newdir/$lib + nnl=`otool -L $ff|grep $lib|awk '{print $1}'` + case $nnl in + /* ) old=$nnl + esac +# echo "old .. " $old $nnl +# otool -L $ff|awk "/$lib/ "'{print $1}' + nn=`otool -L $ff|grep $old|wc -l` +# echo "nn=" $nn $nnl $old + otool -L $ff| grep brew + if test $nn -ne 0 ; then +# echo ..change $old to $new in $ff + ch=1 + echo change $old to $new1 in $ff + install_name_tool -change $old $new $ff + fi + done + if test $ch -eq 1 ; then + codesign --remove-signature $ff + codesign -s - $ff + fi + done +else + echo "$olddir/$lib1" "$olddir/$lib2" "$olddir/$lib3" do some not exist! + echo "$newdir/$lib1" "$newdir/$lib2" "$newdir/$lib3" do not exist! +fi \ No newline at end of file diff --git a/bin/compile-version-MacOS b/bin/compile-version-MacOS new file mode 100755 index 000000000..54411a578 --- /dev/null +++ b/bin/compile-version-MacOS @@ -0,0 +1,208 @@ +exec 3>&1 # save stdout stream +echo2() { echo "$@"; echo "$@" 1>&3; } +gitbranch=$1 +dircompilation=$2 +optionconf= +nameopt= +export dirdmgpkg= +case $3 in +*generic*) optionconf+=--enable-generic;; +*debug*) optionconf+=--enable-debug;nameopt+=-g;; +*optim*) optionconf+=--enable-optim;nameopt+=-O;; +esac + +echo " option de configuration" $optionconf +echo ========================================================== + +## stop on error .... + +set -e +#-x +dirscript=`dirname $0` +case $dirscript in +".") dirscript=$PWD;; +[^/]*) dirscript=$PWD/"$dirscript";; +esac +# compute the number of jobs possible ... +njob=$(( $(sysctl -n hw.ncpu)*2 )) +# echo $dirscript +mvdylib=$dirscript/move-dylib +changedylibgfortran=$dirscript/change-dylib-gfortran +#myver=`git describe --tags` +arch=`arch` +case $arch in +*x86*) ARCH=Intel;; +*i386*) ARCH=Intel;; +*arm64*) ARCH=Apple-Silicon;; +*) echo2 "unknow arch, stop" ; exit 1;; +esac; +echo ARCH = $ARCH ... +tmp=/tmp +cd $dircompilation + +test -d $gitbranch-ff++ || git clone -b $gitbranch https://github.com/FreeFem/FreeFem-sources.git $gitbranch-ff++ +cd $gitbranch-ff++ +ffscr=$dircompilation/$gitbranch-ff++ +echo +ver=`git describe --tags` +verb=`git describe --tags| sed s/v// | sed 's/-.*$//'` +dirpkg=$dircompilation/pkg-ff-$gitbranch-$ver +dirdmgpkg=FreeFEM-$ver-$ARCH$nameopt +LOG=$dircompilation/log-ff-$ver-$ARCH$nameopt.log +echo the LOG file $LOG +exec 1>$LOG 2>&1 # save log is $LOG + +echo2 compilation FreeFEM version ${ver} on Mac ${arch} in $(pwd) +echo2 ========================================================== +prefix=/Applications/FreeFem++.app/Contents/ff-$verb +echo2 prefix = $prefix +test -d $prefix || mkdir -p $prefix +## check gfortran +gfortran --version +test -d $prefix/gnu/ || mkdir $prefix/gnu/ +echo2 copy lib fortran in $prefix/gnu/ +## copy lib fortran in $prefix/gnu/ +echo2 ========================================================== +( +fftmp=/tmp/ff++-tmp-$$ +test -d fftmp ||mkdir $fftmp +cd $ffsrc +printf " program sizeofint\n integer i\n end\n" >$fftmp/ffff.f +gfortran $fftmp/ffff.f -o $fftmp/ffff +(export DYLD_PRINT_LIBRARIES=1;$fftmp/ffff ) 2>&1 | awk '{print $NF}'| egrep -v '/usr/lib|$fftmp/ffff|grep .dylib' >$prefix/gnu/list-dylib-gfortran +rm -rf $fftmp +) +brewlibgfortran=`grep gfortran $prefix/gnu/list-dylib-gfortran` +brewlibquadmath=`grep libquadmath $prefix/gnu/list-dylib-gfortran` +( +cd $(dirname $brewlibgfortran) +cp -f libgfortran.5.dylib libquadmath.0.dylib libgcc_s.1.1.dylib $prefix/gnu/ +cd $prefix/gnu/ +ln -sf libgfortran.?.dylib libgfortran.dylib +ln -sf libquadmath.?.dylib libquadmath.dylib +## inutile ???? je ne retouve plus cette lib ???? FH.. +test -d /opt/homebrew && libemutls_w=`find /opt/homebrew -name libemutls_w.a` +if [ -f "$libemutls_w" ] ; then + cp $libemutls_w . +fi +$mvdylib $prefix/gnu $prefix/gnu/*.dylib +) +dirlibgfortranold=$(dirname $brewlibgfortran) +dirlibgfortrannew=$prefix/gnu +libgfortrandylib=$(basename $brewlibgfortran) +libquadmathdylib=$(basename $brewlibquadmath) +cd - +test -d pkg || mkdir -v pkg +( +cd pkg +test -f mpich-4.1.2.tar.gz || wget https://www.mpich.org/static/downloads/4.1.2/mpich-4.1.2.tar.gz +## compile mpi https://www.mpich.org/static/downloads/4.1.2/mpich-4.1.2.tar.gz + +## compile gsl https://ftp.gnu.org/gnu/gsl/gsl-2.7.1.tar.gz +test -f gsl-2.7.1.tar.gz || wget https://ftp.gnu.org/gnu/gsl/gsl-2.7.1.tar.gz +## compile hdf5 https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.14/hdf5-1.14.2/src/hdf5-1.14.2.tar.bz2 +test -f hdf5-1.14.2.tar.bz2 || wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.14/hdf5-1.14.2/src/hdf5-1.14.2.tar.bz2 + +tar zxf gsl-2.7.1.tar.gz +tar zxf hdf5-1.14.2.tar.bz2 +tar zxf mpich-4.1.2.tar.gz +test -f $prefix/lib/pkgconfig/gsl.pc || ( +cd gsl-2.7.1 +./configure --prefix=$prefix +make -j $njob +make install +) +test -f $prefix/bin/h5c++ || ( +cd hdf5-1.14.2 +./configure '--enable-cxx' 'CC=clang' 'CXX=clang++' --prefix=$prefix +make -j $njob +make install +) + +test -f $prefix/lib/pkgconfig/mpich.pc || ( +cd mpich-4.1.2 +./configure --prefix=$prefix "FFLAGS=-fallow-argument-mismatch -L$prefix/gnu -lgfortran -lquadmath" 'CC=clang' 'CXX=clang++' +make -j $njob +make install +) +) +## configuration of Freefem++ et compilation of ff++ +cd $gitbranch-ff++ +test -x src/nw/FreeFem++ || ( +echo2 ================================ +echo2 configuration de freefem++ +echo2 ================================ +echo2 compilation de freefem++ +echo2 ================================ +export PATH=$prefix/bin:$PATH +autoreconf -i +./configure "--prefix=$prefix" '--enable-download' $optionconf '--enable-m64' 'CXXFLAGS=-Wno-undefined-var-template' "FLIBS=-L$prefix/gnu -lgfortran -lquadmath" +echo2 ================================ +echo2 compilation de petc +echo2 ================================ +test -d $prefix/ff-petsc && ( + echo2 $prefix/ff-petsc try to reuse ???? +) +test -d $prefix/ff-petsc || ( +test -f 3rdparty/ff-petsc/done || (cd 3rdparty/ff-petsc && make + mkdir $prefix/ff-petsc + chown -R $USER $prefix/ff-petsc + make petsc-slepc +touch done ) +) + +echo2 ================================ +echo2 compilation de freefem++ +echo2 ================================ + +./reconfigure +make -j $njob +) + +echo2 ================================ +echo2 instalation de freefem++ in $dirpkg +echo2 ================================ +## copy pre installed in DESTDIR ($dirpkg) +mkdir -p $dirpkg/$prefix || echo2 dir $dirpkg/$prefix exist +echo2 rsync -av $prefix/. $dirpkg/$prefix/. +rsync -av $prefix/. $dirpkg/$prefix/. +make install DESTDIR=$dirpkg +# remove ref on fg +echo2 ================================ +echo2 fix gfortran lib distribution $dirpkg/$prefix +echo2 ================================ +tmpdylib=/tmp/dylib-ff-$$.txt +tmplst=/tmp/otool-ff-$$.txt +echo2 find $dirpkg/$prefix -name '*.dylib' +echo2 .. $changedylibgfortran $dirlibgfortranold $dirlibgfortrannew $libgfortrandylib $libquadmathdylib +find $dirpkg/$prefix -name '*.dylib' > $tmpdylib + $changedylibgfortran $dirlibgfortranold $dirlibgfortrannew $libgfortrandylib $libquadmathdylib `cat $tmpdylib` + ## verification + echo2 dylib to brew in installation for debug !!!! + cat $tmpdylib| xargs otool -L | egrep -i 'brew|/opt/' ||echo2 OK + +echo2 ================================ +## build of PACKAGE directory +echo2 build of PACKAGE/DMG directory /tmp/$dirdmgpkg +echo2 ================================ + + +test -d /tmp/$dirdmgpkg || ( +mkdir /tmp/$dirdmgpkg +cd /tmp/$dirdmgpkg +mv $dirpkg/Applications/FreeFem++.app /tmp/$dirdmgpkg +cp $dirscript/script/Install-app.sh /tmp/$dirdmgpkg +fftestedp=`find . -name Laplace.edp` +sed -e "s/@VV@/$ver/" -e "s%@APPFF@%$dirdmgpkg%" -e s%@FFTEST@%$fftestedp% <$dirscript/script/README.md.in > README.md + +) +cp $dirscript/script/Install-app.sh /tmp/$dirdmgpkg + +exec 1>&3 # save log is closo $LOG +cat $LOG >> /tmp/$dirdmgpkg/build-log.log + +# build dmg file $HOME/desktop/$dirdmgpkg.dmg +echo hdiutil create -srcfolder /tmp/$dirdmgpkg -fs HFS+ $HOME/desktop/$dirdmgpkg.dmg +test -f $HOME/desktop/$dirdmgpkg.dmg && rm $HOME/desktop/$dirdmgpkg.dmg +hdiutil create -srcfolder /tmp/$dirdmgpkg -fs HFS+ $HOME/desktop/$dirdmgpkg.dmg +echo /tmp/$dirdmgpkg \ No newline at end of file diff --git a/bin/move-dylib b/bin/move-dylib index c78f86bb4..bbad5748d 100755 --- a/bin/move-dylib +++ b/bin/move-dylib @@ -5,8 +5,9 @@ shift echo move to = $to for ff in $@; do file=`greadlink -f $ff` -echo move-dylib: change dylib in $file :: $ff +echo move-dylib: change dylib in $file :: $ff version 1 chmod a+w $file +#codesign --remove-signature $file install_name_tool -id $file $file for dylib in `otool -L $file | awk '/compatibility/ {print $1}' `; do echo "# $dylib " diff --git a/bin/script/Install-app.sh b/bin/script/Install-app.sh new file mode 100755 index 000000000..f68bbe573 --- /dev/null +++ b/bin/script/Install-app.sh @@ -0,0 +1,99 @@ +#!/bin/bash +# Apple +# hdiutil create -srcfolder /Users/hecht/Desktop/FF-dmg/FreeFEM-4.12-Apple-M1 -fs HFS+ FreeFEM-4.12-Apple-M1.dmg + +scriptpath=$0 +cd `dirname $scriptpath` +ppwd=`pwd` +vffapp=$ppwd/FreeFem++.app +vffver=$ppwd/FreeFem++.app/Contents/ff-* +ffver=`basename $vffver` +ver=`expr "//$ffver" : '//ff-\(.*\)'` +vffexe=$vffver/bin/FreeFem++ +tyexe=`file -b $vffexe` +case $tyexe-`arch` in + *arm64*-*arm64*) ARCH=Apple-Silicon ;echo same arch arm64 arm64 ok;; + *x86*-*x86*) ARCH=Intel; echo same arch x86 x86 ok;; + *x86*-*i386*) ARCH=Intel; echo same arch x86 i386 ok;; + *) echo bad arch : sorry, abort!; exit 1;; + esac +#sudo cp -rf FreeFem++.app /Applications +#sudo xattr -rc /Applications/FreeFem++.app +# echo verif arch +echo Install FreeFEM version $ver in /Applications/FreeFem++.app/Contents/$ffver/bin ARCH $ARCH + +cd +filerc=".profile" +case "$SHELL" in + *zsh) filerc=".zprofile";; + *bash) filerc=".bash_profile";; + *) echo erreur SHELL $SHELL; exit 1;; +esac; +echo filerc : $filerc + +dirff=/Applications/FreeFem++.app/Contents/$ffver/bin +ffexec=/Applications/FreeFem++.app/Contents/$ffver/bin/FreeFem++ +ff=$(which FreeFem++) +dd=$(dirname "$ff") + +if [ ! -x $ffexec ] ; then + echo " No FreeFem++.app $ffver so copy FreeFem++.app in /Application (with sudo do be sure)" + echo ============================= + sudo cp -rf $vffapp /Applications + test -L /usr/local/bin/FreeFem++-CoCoa && sudo rm /usr/local/bin/FreeFem++-CoCoa + sudo cp -f /Applications/FreeFem++.app/Contents/$ffver/bin/FreeFem++-CoCoa /usr/local/bin +fi + +if test ! -d "$dd" ; then + echo " bug no directory $dd !" + exit 1; +fi + +if [ -d $dirff ] ; then + # verification of com.apple.quarantine + quarantine=`xattr -l /Applications/FreeFem++.app/Contents/$ffver/bin/FreeFem++|grep com.apple.quarantine|wc -l` + if [ $quarantine -ne 0 ] ; then + echo " remove: com.apple.quarantine from FreeFem++ app (need sudo )" + echo do sudo xattr -rc /Applications/FreeFem++.app/ + sudo xattr -rc /Applications/FreeFem++.app/ + echo " verif quarantine " + xattr -l /Applications/FreeFem++.app/Contents/$ffver/bin/FreeFem++ + fi +else + echo " OK : No quarantine " +fi +# echo test FreeFem++ with no graphics +/Applications/FreeFem++.app/Contents/$ffver/bin/FreeFem++ /Applications/FreeFem++.app/Contents/$ffver/share/FreeFEM/$ver/examples/tutorial/Laplace.edp -nw + +if [ $dd != $dirff ] ; then + echo Update $filerc + echo ============================= + + test -e $filerc || touch $filerc + echo add export PATH=$dirff:\$PATH " >>$filerc" + grep -v FreeFem++.app/Contents $filerc >$filerc.new + echo export PATH=$dirff:\"\$PATH\" >>$filerc.new + diff -u $filerc.new $filerc + if test $? -eq 1; then + mv -f $filerc.new $filerc + chmod a+x $filerc + else + rm $filerc.new + fi + +fi + +cd ~/Desktop +echo Create Link of example and Doc +echo ============================= +test -L 'FreeFEM doc and Examples' && rm 'FreeFEM doc and Examples' +test -e 'FreeFEM doc and Examples' || ln -s /Applications/FreeFem++.app/Contents/$ffver/share/FreeFEM/ 'FreeFEM doc and Examples' +cd - + + +vvv=`$dirff/FreeFem++ -nw| grep 'version'| wc -l` +if [ $vvv -eq 0 ]; then + echo Error missing install missing lib or quarantine ? +else + echo FreeFem++ work +fi \ No newline at end of file diff --git a/bin/script/README.md.in b/bin/script/README.md.in new file mode 100644 index 000000000..6ed1bc3e5 --- /dev/null +++ b/bin/script/README.md.in @@ -0,0 +1,41 @@ +Version @VV@ of FreeFEM +======================= + +Be carefull the Macos security (com.apple.quarantine) + +To install FreeFEM + +- remove (into dustbin) by hand your old FreeFem++ from the application folder +- mount the dmg file by a double click +- Launch a new terminal windows (in folder Utility) +- type + cd /Volumes/@APPFF@ + +- then type + bash ./Install-app.sh + +- try by typing + FreeFem++-Cocoa + +You should receive a message beginning with a few lines and then +Usage: /Applications/FreeFem++.app/Contents/ff-@VV@/bin/FreeFem++ [FreeFEM arguments] filename [script arguments] +FreeFEM arguments: + +A test to see if FreeFem++ works without graphic +/Applications/FreeFem++.app/Contents/ff-@VV@/bin/FreeFem++ /Applications/FreeFem++.app/@FFTEST@ -nw + + +- If it does not work it is most likely due to the security system of Apple. See https://support.apple.com/en-us/HT202491 +or +https://support.apple.com/fr-fr/HT202491 + +if you have a message saying that FreeFem++.app is corrupted and it should be sent to the dustbin then cancel the message box and type in the terminal + sudo xattr -rc /Applications/FreeFem++.app + +and retry and recheck the Autorisation (see https://support.apple.com/en-us/HT202491 ) + +You can also go to the system settings (under the apple at the top left in the menu) and select the confidentiality and security item and then down at the end you should see something like "open anyway". Sorry! This is because some libraries are in Fortran and there is no Apple Fortran so Freefem can't be signed with Apple. + + + + diff --git a/configure.ac b/configure.ac index c45860d3f..8cbc75124 100644 --- a/configure.ac +++ b/configure.ac @@ -9,7 +9,7 @@ # number in debian/changelog. -AC_INIT([FreeFEM],[4.13],[frederic.hecht@sorbonne-universite.fr],[FreeFEM]) +AC_INIT([FreeFEM],[4.14],[frederic.hecht@sorbonne-universite.fr],[FreeFEM]) dnl : ${CFLAGS=""} dnl : ${CXXFLAGS=""} dnl : ${FCFLAGS=""} @@ -2449,8 +2449,8 @@ fi ##echo @@@@@@@@@@@@ $ff_petsc_ok @@ "hpddm.$DYLIB_SUFFIX $FF_PETCS_DYLIB" TOOL_DISABLE(hpddm,["hpddm.$DYLIB_SUFFIX hpddm_substructuring.$DYLIB_SUFFIX $FF_PETSC_DYLIB"]) -TOOL_DISABLE(bem,["bem.$DYLIB_SUFFIX RadiativeTransfer_htool.$DYLIB_SUFFIX"]) -if test -n "$TOOL_COMPILE_bem" -a -z "$ff_HTOOL_INCLUDE"; then +TOOL_DISABLE(bemtool,["bem.$DYLIB_SUFFIX RadiativeTransfer_htool.$DYLIB_SUFFIX"]) +if test -n "$TOOL_COMPILE_bemtool" -a -z "$ff_HTOOL_INCLUDE"; then TOOL_COMPILE_htool="htool" fi ff_inc_htool=$ff_HTOOL_INCLUDE @@ -2461,7 +2461,7 @@ fi if test -n "$ff_inc_htool"; then ff_inc_htool="$ff_inc_htool -DWITH_bemtool -DWITH_htool" fi -if test -z "$TOOL_COMPILE_bem"; then +if test -z "$TOOL_COMPILE_bemtool"; then enable_bemtool=no ff_inc_htool="" fi diff --git a/examples/3dCurve/tutomesh1d.edp b/examples/3dCurve/tutomesh1d.edp index eed60beae..8561190d7 100644 --- a/examples/3dCurve/tutomesh1d.edp +++ b/examples/3dCurve/tutomesh1d.edp @@ -1,7 +1,7 @@ load "msh3" // bizarre int nn= 4; int[int] ll=[10,11]; -meshL Th= Sline(nn,[x*nn,0,0],label=ll); // +meshL Th= segment(nn,[x*nn,0,0],label=ll,region=1); // // print information of Th // bug label corrige ?? cout << " Element : "< 1.0E-6) + ThVelocity = trunc(ThPressure, abs(gamma) < 1.0E-2, label = -111111, new2old = n2o); + else { + ThVelocity = ThPressure; + n2o.resize(ThVelocity.nt); + n2o = 0:ThVelocity.nt - 1; + } + fespace WhUnneeded(ThPressure, [Vk, Vk]); // again, unneeded velocity space since this has two layers of overlap + int[int] restriction = restrict(Wh, WhUnneeded, n2o); + int[int] I = 0:restriction.n - 1; + real[int] C(I.n); + C = 1; + two2one = [I, restriction, C]; + two2one.resize(Wh.ndof, WhUnneeded.ndof); +} +Mat A(AUnneeded, restriction = two2one); // proper Mat for the velocity with a single layer of overlap + +real f = -100.0; +func nu = int(y * n) % 2 == 0 ? 0.4999 : 0.35; +func E = int(y * n) % 2 == 0 ? 1.0E+7 : 200.0E+9; +func mu = E / (2.0*(1.0 + nu)); +func lambda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu)); +varf vPbA(def2(u), def2(v)) = int2d(ThVelocity)(2.0 * mu * (epsilon(u)' * epsilon(v))) + int2d(ThVelocity)(f * vB) + on(4, u = 0.0, uB = 0.0); +varf vPbC([p], [q]) = int2d(ThPressure)(-1.0/lambda * p * q); +varf vPbB(def2(u), [q]) = int2d(ThPressure)(div(u) * q); + +tgv = -2.0; +matrix NA; +matrix NB = vPbB(Wh, Ph); +NB.thresholding(1.0E-16); +matrix NC = vPbC(Ph, Ph); +NC.thresholding(1.0E-16); + +real[int] rhs(Wh.ndof + Ph.ndof); +rhs(0:Wh.ndof - 1) = vPbA(0, Wh); +rhs(Wh.ndof:Wh.ndof + Ph.ndof - 1) = 0.0; +if(tgv < 0.0) { + assert(abs(tgv + 2.0) < 1.0E-6); + NA = vPbA(Wh, Wh, tgv = -10.0); + NA.thresholding(1.0E-16); + varf vPbRHS(def2(u), def2(v)) = on(4, u = 0.0, uB = 0.0); + real[int] tmp = vPbRHS(0, Wh); + real[int] rhsA = NA * tmp; + real[int] rhsB = NB * tmp; + rhs(0:Wh.ndof - 1) -= rhsA; + rhs(Wh.ndof:rhs.n - 1) -= rhsB; + varf vPbG(def2(u), def2(v)) = on(4, u = 1.0, uB = 1.0); + rhsA = vPbG(0, Wh); + setBC(NA, rhsA, -2.0); + if(sym) + symmetrizeCSR(NA); + for [i, j, v : NB] if(abs(rhsA(j) - 1.0) < 1.0E-2) v = 0.0; +} +else { + NA.thresholding(1.0E-16); + NA = vPbA(Wh, Wh); +} + +A = NA; +Mat B(C, A, NB); +C = NC; +Mat Elas = [[A, B'], + [B, C ]]; +set(Elas, sparams = "-pc_type none"); +if(HasType("PC", "hpddm")) { + string common = " -ksp_converged_reason -ksp_view_final_residual -ksp_monitor_true_residual -ksp_type fgmres"; + set(Elas, sparams = common + " -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self " + + "-prefix_push fieldsplit_1_ " + + " -ksp_type fgmres -ksp_rtol 1.0E-1 -pc_type hpddm -pc_hpddm_has_neumann -ksp_monitor -pc_hpddm_ksp_pc_side right -pc_hpddm_ksp_rtol 1.0E-1 -pc_hpddm_schur_precondition geneo " + + " -prefix_push pc_hpddm_levels_1_ " + + " -eps_nev 10 -eps_threshold 0.5 -st_share_sub_ksp -eps_gen_non_hermitian -sub_pc_factor_mat_solver_type mumps -sub_pc_type cholesky -eps_tol 1.0E-2 " + + " -prefix_pop -prefix_push pc_hpddm_coarse_ " + + " -correction balanced -pc_type cholesky -pc_factor_mat_solver_type mumps " + + " -prefix_pop" + + " -prefix_pop", setup = 1); + set(A, parent = Elas, bs = sym ? 1 : 2, sparams = + "-prefix_push fieldsplit_0_ -pc_type hpddm -ksp_pc_side right -pc_hpddm_has_neumann -pc_hpddm_define_subdomains" + + " -prefix_push pc_hpddm_levels_1_ " + + " -eps_nev 20 -eps_threshold 0.3 -st_share_sub_ksp -sub_pc_factor_mat_solver_type mumps -sub_pc_type cholesky " + + " -prefix_pop -prefix_push pc_hpddm_coarse_ " + + " -pc_type cholesky -pc_factor_mat_solver_type mumps " + + " -prefix_pop" + + " -prefix_pop"); + set(C, parent = Elas); + + real[int] sol(Wh.ndof + Ph.ndof); + sol = Elas^-1 * rhs; + Wh def2(u); + u[] = sol(0:Wh.ndof - 1); + Ph p; + p[] = sol(Wh.ndof:sol.n - 1); + plotMPI(ThVelocity, def2(u), [Vk, Vk], def2, real, cmm = "Global velocity (overlap = 1)"); + real alpha = 1.0E4; + mesh ThMoved = movemesh(ThVelocity, [x + alpha * u, y + alpha * uB]); + u[] = mpirank; + plotMPI(ThMoved, u, P2, unneeded, real, cmm = "Global moved solution"); + plotD(ThPressure, p, cmm = "Global pressure (overlap = 2)"); + p = nu; + plotMPI(ThPressure, p, P0, unneeded, real, cmm = "Material coefficient"); +} diff --git a/examples/hpddm/laplace-RT-2d-PETSc.edp b/examples/hpddm/laplace-RT-2d-PETSc.edp index 451d96899..7bbe39a09 100644 --- a/examples/hpddm/laplace-RT-2d-PETSc.edp +++ b/examples/hpddm/laplace-RT-2d-PETSc.edp @@ -23,7 +23,7 @@ fespace Wh(Th, Pk); // local finite element space Mat A; int[int] n2o; macro ThN2O()n2o// -DmeshBuild(Th); +DmeshCreate(Th); { macro ThPostProcessD(D) { Wh def(u), def(v); diff --git a/examples/hpddm/laplace-RT-3d-PETSc.edp b/examples/hpddm/laplace-RT-3d-PETSc.edp index 2e80a9321..770ac0133 100644 --- a/examples/hpddm/laplace-RT-3d-PETSc.edp +++ b/examples/hpddm/laplace-RT-3d-PETSc.edp @@ -14,7 +14,7 @@ mesh3 ThGlobal = Th; Mat A; int[int] n2o; macro ThN2O()n2o// -DmeshBuild(Th); +DmeshCreate(Th); func Pk = [RT03d, P0]; fespace Wh(Th, Pk); { diff --git a/examples/hpddm/laplace-adapt-3d-PETSc.edp b/examples/hpddm/laplace-adapt-3d-PETSc.edp index 1c0eb37f6..6e3288922 100644 --- a/examples/hpddm/laplace-adapt-3d-PETSc.edp +++ b/examples/hpddm/laplace-adapt-3d-PETSc.edp @@ -22,7 +22,7 @@ mesh3 ThBackup = Th3; Mat A; int[int] n2o; NewMacro Th3N2O() n2o EndMacro -DmeshBuild(Th3); +DmeshCreate(Th3); MatCreate(Th3, A, P1); fespace Vh(Th3, P1); Vh u; diff --git a/examples/hpddm/laplace-adapt-dist-3d-PETSc.edp b/examples/hpddm/laplace-adapt-dist-3d-PETSc.edp index de0ef7c19..3537134ad 100644 --- a/examples/hpddm/laplace-adapt-dist-3d-PETSc.edp +++ b/examples/hpddm/laplace-adapt-dist-3d-PETSc.edp @@ -40,7 +40,7 @@ metParMmg = met(rest); varf vPoisson(u, v) = int3d(Th3)(grad(u)' * grad(v)) + int3d(Th3)(v) + on(1,2,3, u=0); Mat A; -DmeshBuild(Th3); +DmeshCreate(Th3); MatCreate(Th3, A, P1); fespace Vh(Th3, P1); Vh u; diff --git a/examples/hpddm/laplace-lagrange-PETSc.edp b/examples/hpddm/laplace-lagrange-PETSc.edp index 428f7eefe..855bc4230 100644 --- a/examples/hpddm/laplace-lagrange-PETSc.edp +++ b/examples/hpddm/laplace-lagrange-PETSc.edp @@ -58,11 +58,7 @@ if(mpirank == 0) { plotMPI(Th, sol, Pk, def, real, params); } { - real[int] split(pF.n); - split = 1; - if(mpirank == 0) - split(pF.n - 1) = 2; - set(N, sparams = "-ksp_view_pre -pc_type fieldsplit -fieldsplit_0_pc_type lu -fieldsplit_0_mat_mumps_cntl_3 1e-12 -fieldsplit_0_mat_mumps_icntl_24 1 -pc_fieldsplit_detect_saddle_point", fields = split); + set(N, sparams = "-ksp_view_pre -pc_type fieldsplit -fieldsplit_0_pc_type lu -fieldsplit_0_mat_mumps_cntl_3 1e-12 -fieldsplit_0_mat_mumps_icntl_24 1 -pc_fieldsplit_detect_saddle_point"); KSPSolve(N, pF, pSol); if(mpirank == 0) pSol.resize(pSol.n - 1); diff --git a/examples/hpddm/maxwell-3d-PETSc.edp b/examples/hpddm/maxwell-3d-PETSc.edp index fe55afd8c..3ebf8834d 100644 --- a/examples/hpddm/maxwell-3d-PETSc.edp +++ b/examples/hpddm/maxwell-3d-PETSc.edp @@ -24,7 +24,7 @@ fespace Wh(Th, Pk); Mat A; int[int] n2o; macro ThN2O()n2o// EOM -DmeshBuild(Th); +DmeshCreate(Th); { macro ThPkPart()Edge03ds0// EOM special FE for the domain partitioning MatCreate(Th, A, Pk); diff --git a/examples/hpddm/maxwell-3d-surf-PETSc-complex.edp b/examples/hpddm/maxwell-3d-surf-PETSc-complex.edp index 7aa523f74..8cb824ce4 100644 --- a/examples/hpddm/maxwell-3d-surf-PETSc-complex.edp +++ b/examples/hpddm/maxwell-3d-surf-PETSc-complex.edp @@ -34,7 +34,7 @@ Mat A; int[int] n2o; macro ThSN2O()n2o// ThG = ThS; -DmeshBuild(ThS); +DmeshCreate(ThS); { Mat H, T; MatCreate(ThS, H, P1); @@ -85,7 +85,7 @@ if(!NoGraphicWindow) { meshS ThOutS = ThOutG; int[int] n2o;// macro ThOutSN2O()n2o// - DmeshBuild(ThOutS); + DmeshCreate(ThOutS); { Mat B, T; MatCreate(ThOutS, B, P1); diff --git a/examples/hpddm/mf-2d-SLEPc.edp b/examples/hpddm/mf-2d-SLEPc.edp index bfba0740f..a63c9f70d 100644 --- a/examples/hpddm/mf-2d-SLEPc.edp +++ b/examples/hpddm/mf-2d-SLEPc.edp @@ -8,7 +8,7 @@ mesh Th = square(40, 40); func Pk = P1; fespace Vh(Th, Pk); varf vPb(u, v) = int2d(Th)(dx(u) * dx(v) + dy(u) * dy(v)) + on(1, 2, 3, 4, u = 0.0); -DmeshBuild(Th); +DmeshCreate(Th); Mat A; MatCreate(Th, A, Pk); matrix Loc = vPb(Vh, Vh, tgv = -2); diff --git a/examples/hpddm/oseen-2d-PETSc.edp b/examples/hpddm/oseen-2d-PETSc.edp index 94e15ad75..89ace5885 100644 --- a/examples/hpddm/oseen-2d-PETSc.edp +++ b/examples/hpddm/oseen-2d-PETSc.edp @@ -29,7 +29,7 @@ func b2 = -2 * (2*x - 1) * (1 - (2*y - 1)^2); // mesh //-------------------- mesh Th = square(n, n); -DmeshBuild(Th); +DmeshCreate(Th); fespace Ph(Th, P0); Ph part; PartitionCreate(Th, part[], P0); diff --git a/examples/hpddm/restriction-2d-PETSc.edp b/examples/hpddm/restriction-2d-PETSc.edp index 0cf023e2b..2c7cbb2fc 100644 --- a/examples/hpddm/restriction-2d-PETSc.edp +++ b/examples/hpddm/restriction-2d-PETSc.edp @@ -41,7 +41,7 @@ func Pk = P2; } macro ThUserPartitioning()part// EOM - DmeshBuild(Th); + DmeshCreate(Th); MatCreate(Th, A, Pk); } diff --git a/examples/hpddm/stokes-2d-SLEPc.edp b/examples/hpddm/stokes-2d-SLEPc.edp index 029a3b810..1b1386295 100644 --- a/examples/hpddm/stokes-2d-SLEPc.edp +++ b/examples/hpddm/stokes-2d-SLEPc.edp @@ -21,7 +21,7 @@ int s = getARGV("-split", 1); // refinement factor mesh Th = square(getARGV("-global", 40), getARGV("-global", 40), [x, y]); fespace VhV(Th, PkV); // local finite element space for the velocities fespace VhP(Th, PkP); // local finite element space for the pressure -DmeshBuild(Th); +DmeshCreate(Th); Mat A, M; { macro def(u)def2(u)// diff --git a/examples/hpddm/stokes-adapt-3d-PETSc.edp b/examples/hpddm/stokes-adapt-3d-PETSc.edp new file mode 100644 index 000000000..ae4d6747e --- /dev/null +++ b/examples/hpddm/stokes-adapt-3d-PETSc.edp @@ -0,0 +1,114 @@ +// run with MPI: ff-mpirun -np 4 script.edp +// NBPROC 4 + +load "msh3" +load "mshmet" +load "mmg" +load "aniso" +load "PETSc" +macro dimension()3// EOM // 2D or 3D +include "macro_ddm.idp" // additional DDM functions + +macro def(i)[i, i#B, i#C, i#D]// EOM +macro init(i)[i, i, i, i]// EOM +func Pk = [P2, P2, P2, P1]; +macro grad(u)[dx(u), dy(u), dz(u)]// EOM +macro div(u)(dx(u) + dy(u#B) + dz(u#C))// EOM + +real anisomax = getARGV("-anisomax", 4.0); +bool aniso = bool(anisomax > 1.0); +mesh3 Th; +{ + mesh ThGlobal2d = square(getARGV("-global", 12), getARGV("-global", 12), [x, y]); // global mesh + ThGlobal2d = trunc(ThGlobal2d, (x <= 0.5) || (y <= 0.5), label = 5); + ThGlobal2d = trunc(ThGlobal2d, (y >= 0.25) || (x >= 0.25), label = 5); + //mesh Th2d = movemesh(ThGlobal2d, [-x, y]); + //ThGlobal2d = ThGlobal2d + Th2d; + Th = buildlayers(ThGlobal2d, getARGV("-global", 12) / 2, zbound = [0, 0.4]); +} +mesh3 ThG = Th; // need to backup the original global mesh + +Mat A; +int[int] n2o; // need to know how to go from the local to the global mesh +if(mpirank==0) cout << "it = 0, Th.nv = " << ThG.nv << endl; +macro ThN2O()n2o// +MatCreate(Th, A, Pk); + +fespace Wh(Th, Pk); // complete space [u, v, w, p] +fespace Zh(Th, P2); // velocity space +fespace Qh(Th, P1); // pressure space +Qh onlyP; +Zh onlyU; + +varf vPb(def(u), def(v)) = int3d(Th)( + grad(u)' * grad(v) + + grad(uB)' * grad(vB) + + grad(uC)' * grad(vC) + - div(u) * vD - div(v) * uD + ) + + on(0, 1, 3, 5, u = 0, uB = 0, uC = 0) + + on(2, u = 1000*y*(0.5-y)*z*(0.4-z), uB = 0, uC = 0); + +Wh def(u); + +for(int i = 0; i < 5; ++i) { + real[int] rhs = vPb(0, Wh, tgv = -1); + A = vPb(Wh, Wh, tgv = -1); + set(A, sparams = "-pc_type lu"); // add "-pc_factor_mat_solver_type mumps" to run on 1 proc + u[] = A^-1*rhs; + if(!NoGraphicWindow) { + onlyP = uD; + macro def1(i)i// + plotMPI(Th, onlyP, P1, def1, real, cmm = "Pressure for it = " + (i+1)); + onlyU = u; + plotMPI(Th, onlyU, P2, def1, real, cmm = "x-velocity for it = " + (i+1)); + } + { + fespace WhG(ThG, Pk); + fespace ZhG(ThG, P2); + fespace QhG(ThG, P1); + WhG def(uG), def(uReduce); + u[] .*= A.D; + int[int] rest = restrict(Wh, WhG, n2o); + uReduce[](rest) = u[]; + mpiAllReduce(uReduce[], uG[], mpiCommWorld, mpiSUM); + if(mpirank == 0) { + // Since mshmet cannot regularize FE functions like adaptmesh, + // first regularize solutions on [0,1]; + ZhG uMag = sqrt(uG^2+uGB^2+uGC^2); + real unorm = uMag[].max-uMag[].min; + QhG p = uGD; + real pnorm = p[].max-p[].min; + // Define normalized sol'n + def(uReduce) = [uG/unorm, uGB/unorm, uGC/unorm, uGD/pnorm]; + // create metric + real[int] met = mshmet(ThG, [uReduce, uReduceB, uReduceC], uReduceD, normalization = bool(getARGV("-normalization",0)), aniso = aniso, hmin = 1.0e-3, hmax = 0.2, err = 1.0e-2, verbosity = verbosity); + if(aniso) boundaniso(6, met, anisomax); // restrict anisotropy + fespace METhG(ThG,[P1,P1,P1,P1,P1,P1]); + METhG [m11,m21,m22,m31,m32,m33]; + QhG mm; + if(!aniso) mm[] = met; + else { + m11[] = met; + mm = m11+m22+m33; // trace of anisotropic metric tensor + } + plot(mm); + ThG = mmg3d(ThG, metric = met, hmin = 1.0e-3, hmax = 0.2, hgrad = -1, verbose=verbosity-(verbosity==0)); + cout << "it = " << (i+1) << ", Th.nv = " << ThG.nv << endl; + plot(ThG); + } + broadcast(processor(0), ThG); + def(uG) = def(uG); + Th = ThG; + Mat Adapt; + MatCreate(Th, Adapt, Pk) // decompose the adapted mesh + A = Adapt; // replace the old Jacobian + def(u) = init(0.0); // just to trigger a resize + rest.resize(u[].n); + rest = restrict(Wh, WhG, n2o); + u[] = uG[](rest); + } +} +DmeshSave(Th, "stokes-adapt-3d"); +ofstream sol("stokes-adapt-3d_" + mpirank + "_" + mpisize + ".sol"); +sol << u[]; diff --git a/examples/hpddm/stokes-block-2d-PETSc.edp b/examples/hpddm/stokes-block-2d-PETSc.edp index cfb002015..7ce8274e7 100644 --- a/examples/hpddm/stokes-block-2d-PETSc.edp +++ b/examples/hpddm/stokes-block-2d-PETSc.edp @@ -27,7 +27,7 @@ Mat dA, dC; ThGlobal = trunc(ThGlobal, (x < 0.5) || (y < 0.5), label = 5); Th = movemesh(ThGlobal, [-x, y]); Th = ThGlobal + Th; - DmeshBuild(Th); + DmeshCreate(Th); { macro def(u)def2(u)// macro init(u)init2(u)// diff --git a/examples/hpddm/stokes-block-hpddm-2d-PETSc.edp b/examples/hpddm/stokes-block-hpddm-2d-PETSc.edp new file mode 100644 index 000000000..310f7a516 --- /dev/null +++ b/examples/hpddm/stokes-block-hpddm-2d-PETSc.edp @@ -0,0 +1,98 @@ +load "PETSc" +include "macro_ddm.idp" + +macro def2(i)[i, i#B]// +macro init2(i)[i, i]// +macro grad(u)[dx(u), dy(u)]// +macro div(u)(dx(u) + dy(u#B))// + +int global = getARGV("-global", 50); +bool sym = usedARGV("-sym") != -1; +int[int] l = [2, 2, 1, 2]; +mesh Th = square(global, global, [x, y], label = l, flags = 1); // global mesh +func Vk = P2; +func Pk = P1; + +macro ThRefinementFactor()getARGV("-split", 1)// +DmeshCreate(Th); +Mat A, C; +{ + macro def(i)def2(i)// + macro init(i)init2(i)// + MatCreate(Th, A, [Vk, Vk]); +} +{ + MatCreate(Th, C, Pk); +} +fespace Wh(Th, [Vk, Vk]); +fespace Ph(Th, Pk); + +varf vPbA(def2(u), def2(v)) = int2d(Th)(grad(u)' * grad(v) + grad(uB)' * grad(vB)) + on(1, u = 1.0, uB = 0.0) + on(2, u = 0.0, uB = 0.0); +varf vPbB(def2(u), [q]) = int2d(Th)(-div(u) * q); + +tgv = -2.0; +matrix NA; +matrix NB = vPbB(Wh, Ph); +NB.thresholding(1.0E-16); + +real[int] rhs(Wh.ndof + Ph.ndof); +rhs(0:Wh.ndof - 1) = vPbA(0, Wh); +rhs(Wh.ndof:Wh.ndof + Ph.ndof - 1) = 0.0; +if(tgv < 0.0) { + assert(abs(tgv + 2.0) < 1.0E-6); + NA = vPbA(Wh, Wh, tgv = -10.0); + if (!sym) + NA.thresholding(1.0E-16); + varf vPbRHS(def2(u), def2(v)) = on(1, u = 1.0, uB = 0.0) + on(2, u = 0.0, uB = 0.0); + real[int] tmp = vPbRHS(0, Wh); + real[int] rhsA = NA * tmp; + real[int] rhsB = NB * tmp; + rhs(0:Wh.ndof - 1) -= rhsA; + rhs(Wh.ndof:rhs.n - 1) -= rhsB; + varf vPbG(def2(u), def2(v)) = on(1, 2, u = 1.0, uB = 1.0); + rhsA = vPbG(0, Wh); + setBC(NA, rhsA, -2.0); + if(sym) + symmetrizeCSR(NA); + for [i, j, v : NB] if(abs(rhsA(j) - 1.0) < 1.0E-2) v = 0.0; +} +else { + NA.thresholding(1.0E-16); + NA = vPbA(Wh, Wh); +} + +A = NA; +Mat B(C, A, NB); +Mat Stokes = [[A, B'], + [B, 0 ]]; +set(Stokes, sparams = "-pc_type none"); +if(HasType("PC", "hpddm")) { + string common = " -ksp_converged_reason -ksp_view_final_residual -ksp_monitor_true_residual -ksp_type fgmres"; + set(Stokes, sparams = common + " -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self " + + "-prefix_push fieldsplit_1_ " + + " -ksp_type fgmres -pc_type hpddm -pc_hpddm_has_neumann -ksp_monitor -pc_hpddm_ksp_rtol 1.0E-2 -pc_hpddm_schur_precondition geneo " + + " -prefix_push pc_hpddm_levels_1_ " + + " -eps_nev 20 -eps_threshold 0.5 -st_share_sub_ksp -eps_gen_non_hermitian -sub_pc_factor_mat_solver_type mumps -sub_pc_type cholesky -eps_tol 1.0E-2 " + + " -prefix_pop -prefix_push pc_hpddm_coarse_ " + + " -correction balanced -pc_type cholesky -pc_factor_mat_solver_type mumps " + + " -prefix_pop" + + " -prefix_pop", setup = 1); + set(A, parent = Stokes, bs = sym ? 1 : 2, sparams = + "-prefix_push fieldsplit_0_ -pc_type hpddm -ksp_pc_side right -pc_hpddm_has_neumann -pc_hpddm_define_subdomains" + + " -prefix_push pc_hpddm_levels_1_ " + + " -eps_nev 40 -eps_threshold 0.3 -st_share_sub_ksp -sub_pc_factor_mat_solver_type mumps -sub_pc_type cholesky " + + " -prefix_pop -prefix_push pc_hpddm_coarse_ " + + " -pc_type cholesky -pc_factor_mat_solver_type mumps " + + " -prefix_pop" + + " -prefix_pop"); + set(C, parent = Stokes); + + real[int] sol(Wh.ndof + Ph.ndof); + sol = Stokes^-1 * rhs; + Wh def2(u); + u[] = sol(0:Wh.ndof - 1); + Ph p; + p[] = sol(Wh.ndof:sol.n - 1); + plotMPI(Th, def2(u), [Vk, Vk], def2, real, cmm = "Global velocity (overlap = 1)"); + plotD(Th, p, cmm = "Global pressure (overlap = 1)"); +} diff --git a/examples/hpddm/transfer.edp b/examples/hpddm/transfer.edp index 44e532daf..3c409a260 100644 --- a/examples/hpddm/transfer.edp +++ b/examples/hpddm/transfer.edp @@ -23,8 +23,8 @@ int[int] fforder = [1]; border d(t=0,1) { x=0; y=1-t; }; ThNew = buildmesh(a(N)+b(N)+c(N)+d(N)); } - DmeshBuild(Th); - DmeshBuild(ThNew); + DmeshCreate(Th); + DmeshCreate(ThNew); fespace Vh(Th, Pk); fespace VhNew(ThNew, PkNew); Vh u = cos(8*pi*x) * sin(8*pi*y); @@ -49,8 +49,8 @@ int[int] fforder = [1]; if(mpirank == 0) ThNew = mmg3d(Th, hmax = 8.0e-2); broadcast(processor(0), ThNew); - DmeshBuild(Th); - DmeshBuild(ThNew); + DmeshCreate(Th); + DmeshCreate(ThNew); fespace Vh(Th, Pk); fespace VhNew(ThNew, PkNew); Vh def(u) = [cos(2*pi*x) * sin(2*pi*y) * cos(2*pi*z), cos(2*pi*z)]; diff --git a/examples/hpddm/transpose-solve-PETSc.edp b/examples/hpddm/transpose-solve-PETSc.edp index 4e1e5d7b8..c42d57edd 100644 --- a/examples/hpddm/transpose-solve-PETSc.edp +++ b/examples/hpddm/transpose-solve-PETSc.edp @@ -14,7 +14,10 @@ func Pk = P1; // finite element space mesh Th = square(getARGV("-global", 40), getARGV("-global", 40)); // global mesh -DmeshBuild(Th); +mesh ThBackup = Th; +int[int] n2o; +macro ThN2O()n2o// +DmeshCreate(Th); fespace Wh(Th, Pk); // local finite element space varf vPb(u, v) = intN(Th)(u * v + grad(u)' * grad(v) - 10 * dy(v) * u) + intN(Th)(v * x) + on(1, 2, 3, 4, u = 0.0); @@ -40,3 +43,13 @@ real local = u[].linfty; real norm; mpiAllReduce(local, norm, mpiCommWorld, mpiMAX); plotMPI(Th, v, Pk, def, real, cmm = "Global solution transposed problem (linfty of error = " + norm + ")"); +if(!NoGraphicWindow) { + fespace WhBackup(ThBackup, Pk); + int[int] rest = restrict(Wh, WhBackup, n2o); + u[] = v[]; + u[] .*= A.D; + WhBackup sol, solReduced; + for[i, v : rest] sol[][v] = u[][i]; + mpiReduce(sol[], solReduced[], processor(0), mpiSUM); + plot(solReduced, wait = 1); +} diff --git a/examples/plugin/cavityNewtonP3pnc.edp b/examples/plugin/cavityNewtonP3pnc.edp index c2fe6e470..6245c5e3c 100644 --- a/examples/plugin/cavityNewtonP3pnc.edp +++ b/examples/plugin/cavityNewtonP3pnc.edp @@ -96,7 +96,7 @@ varf vNS ([u1, u2, p], [v1, v2, q]) // DF(u)u - F(u) + on(3, u1=BCu1, u2=0) + on(1, 2, 4, u1=0, u2=0); ; - + real err =0; for (real re = reyini ; re <= reymax; re *= 2) { real lerr = 1; nu = 1./re; @@ -122,14 +122,18 @@ for (real re = reyini ; re <= reymax; re *= 2) { // Error b = u1[]-up1[]; - - cout << "iter = "<< i << ", err = " << b.l2 << ", rey = " << re << endl; + err = b.l2; + cout << "iter = "<< i << ", err = " << err << ", rey = " << re << endl; // Convergence criteria - if (b.l2 < 1e-4) break; + if (err < 1e-4) break; + if(err > 1e10) break; } + if(err > 1e10) break; + } - + if(err > 1e10) break; + // Stream-lines uu1 = u1; uu2 = u2; diff --git a/examples/plugin/test-ElementMixte.edp b/examples/plugin/test-ElementMixte.edp index 5583fbee2..9c0ab57ce 100644 --- a/examples/plugin/test-ElementMixte.edp +++ b/examples/plugin/test-ElementMixte.edp @@ -165,7 +165,7 @@ for( int j=0;j 1e-10) { w1[][i]+=1;cout << j << " " << i << " diff = " << w1[] << endl;} - cout << " ||w - (\delta_ij)_j= || " << w1[].linfty << endl; + if(w1[].linfty > 1e-10) { w1[][i]+=1;cout << j << " " << i << " diff = " << w1[] << " " << err1 << endl;} + cout << j << " ||w - (\delta_ij)_j= || " << w1[].linfty << endl; assert( w1[].linfty < 1e-10); NewMacro diff2(op) int2d(Th)(square(op(u1)-op(U1))+square(op(u2)-op(U2))) EndMacro - assert( diff2( 0+ )) ; + assert( diff2( )) ; assert( diff2(dx)) ; assert( diff2(dy)) ; } -assert(err1==0); +assert(abs(err1)<1e-10); Th=square(5,5,[10*(x+y/3),10*(y-x/3)]); V(u)=[0,0]; @@ -235,5 +235,8 @@ CheckRT(RT2,"RT2",x,y,pxy2,ccn,V,1) CheckRT(RT2Ortho,"RT2Ortho",-y,x,pxy2,cctau,V,1) CheckRT(BDM1,"BDM1",x,y,pxy0,ccn,V,1) CheckRT(BDM1Ortho,"BDM1Ortho",-y,x,pxy0,cctau,V,1) +CheckRT(BDM2,"BDM2",x,y,pxy1,ccn,V,1) +// buggus FH version 4.13 ??? in test !!!! +CheckRT(BDM2Ortho,"BDM2Ortho",-y,x,pxy1,cctau,V,1) } diff --git a/examples/tutorial/Convolution-Sample.edp b/examples/tutorial/Convolution-Sample.edp new file mode 100644 index 000000000..5941b4cae --- /dev/null +++ b/examples/tutorial/Convolution-Sample.edp @@ -0,0 +1,53 @@ +// ----- +// the convolution of fonction: int(Sh) ( fh(x-y) * phih(y) dy ), +// for regularisation +load "msh3" +load "Element_P3" +real eps = 0.02;// support of the function phih ... +func f = x >0 ; //y*x*x+y*y*y+h*tanh(val*(sin(5.0*y)-2.0*x)); + +meshL S1= segment(1); +meshL S10= segment(100); +fespace V1H(S1,P3HL);// P3 hermite .. function C1.. ??? +V1H [rho,drho]=[1-x,0];// fonction axi sur [0,1] +meshL Lh=segment(100); +fespace VLh(Lh,P1); +VLh xh=x,rhoh=rho; +// a trick to plot 1d function as arry .. +plot([xh[],rhoh[]],wait=1,cmm=" the rho function "); + +func r=dist(x,y); + +border csupport(t=0,2*pi){ x=eps*cos(t);y=eps*sin(t);} +mesh Sh=buildmesh(csupport(20)); +fespace Wh(Sh,P1); +Wh rr=min(r/eps,1.);// +real beta2h=1; +Wh phih= beta2h*rho(rr);// Ok version 5/10/2023. +cout << " int int2d(Sh)(phih) "< 0.9 && abs(sigma[][j] - intersection[0][i]) < 0.1) + value = 0.0; + } + } + } + } if(trueRestrict && mpiSize(comm) == size && removeZeros) { assert(level == 1); meshN ThIntersection; - fespace PhIntersectionPrivate(ThIntersection, P0); - PhIntersectionPrivate[int] recv(numberIntersection); - PhIntersectionPrivate[int] send(numberIntersection); mpiRequest[int] rq(2 * numberIntersection); - for(int i = 0; i < numberIntersection; ++i) { - ThIntersection = meshIntersection[i]; - Irecv(processor(intersection[0][i], comm, rq[i]), recv[i][]); - send[i] = khi; - Isend(processor(intersection[0][i], comm, rq[numberIntersection + i]), send[i][]); - } - meshName[0] = trunc(meshName[0], khi > 1.0e-6, label = 9999); - khi = khi; int[int] skip(0); - for(int k = 0; k < 2 * numberIntersection; ++k) { - int i = mpiWaitAny(rq); - if(i < numberIntersection) { + { + fespace PhIntersectionPrivate(ThIntersection, P0); + PhIntersectionPrivate[int] recv(numberIntersection); + PhIntersectionPrivate[int] send(numberIntersection); + for(int i = 0; i < numberIntersection; ++i) { ThIntersection = meshIntersection[i]; - PhIntersectionPrivate intersection = send[i] > 1.0e-6 && recv[i] > 1.0e-6; - if(intersection[].l2 > 1.0e-6) - meshIntersection[i] = trunc(meshIntersection[i], intersection > 1.0e-6, label = 9999); - else { - skip.resize(skip.n + 1); - skip[skip.n - 1] = i; + Irecv(processor(intersection[0][i], comm, rq[i]), recv[i][]); + send[i] = khi; + Isend(processor(intersection[0][i], comm, rq[numberIntersection + i]), send[i][]); + } + real[int] tmpRemove(meshName[0].nt); + int[int] n2oRemove; +IFMACRO(!privateDmesh#N2O) + meshName[0] = trunc(meshName[0], khi > 1.0e-6, label = 9999, new2old = n2oRemove); +ENDIFMACRO +IFMACRO(privateDmesh#N2O) + { + int[int] n2oFilter; + meshName[0] = trunc(meshName[0], khi > 1.0e-6, label = 9999, new2old = n2oRemove); + n2oFilter.resize(n2oRemove.n); + n2oFilter = privateDmesh#N2O(n2oRemove); + privateDmesh#N2O.resize(n2oFilter.n); + privateDmesh#N2O = n2oFilter; + } +ENDIFMACRO + khi = khi; + fespace PhLocalPrivate(meshName[0], P0); + for(int k = 0; k < 2 * numberIntersection; ++k) { + int i = mpiWaitAny(rq); + if(i < numberIntersection) { + ThIntersection = meshIntersection[i]; + PhIntersectionPrivate intersection = send[i] > 1.0e-6 && recv[i] > 1.0e-6; + if(intersection[].l2 > 1.0e-6) { + tmpRemove = 0.0; + tmpRemove(n2oIntersection[i]) = intersection[]; + PhLocalPrivate afterTrunc; + afterTrunc[] = tmpRemove(n2oRemove); + meshIntersection[i] = trunc(meshName[0], afterTrunc > 1.0e-6, label = 9999, new2old = n2oIntersection[i]); + } + else { + skip.resize(skip.n + 1); + skip[skip.n - 1] = i; + } } } } skip.sort; - intersection.resize(1 + numberIntersection - skip.n); + real[int][int] recv(numberIntersection); + real[int][int] send(numberIntersection); int j = 0; + for(int i = 0; i < numberIntersection; ++i) { + bool skipped = false; + if(j < skip.n) { + if(skip[j] == i) { + ++j; + skipped = true; + } + } + if(!skipped) { + recv[i].resize(meshIntersection[i].nv); + Irecv(processor(intersection[0][i], comm, rq[i]), recv[i]); + fespace VhIntersectionPrivate(meshIntersection[i], P1); + int[int] restriction = restrict(VhIntersectionPrivate, VhLocalPrivate, n2oIntersection[i]); + send[i].resize(meshIntersection[i].nv); + send[i] = khi[](restriction); + Isend(processor(intersection[0][i], comm, rq[numberIntersection + i]), send[i]); + } + } + mpiWaitAll(rq); + intersection.resize(1 + numberIntersection - skip.n); + j = 0; for(int i = 0; i < numberIntersection; ++i) { bool skipped = false; if(j < skip.n) { @@ -424,17 +485,20 @@ IFMACRO(!vectorialfe) fespace singleComponentWhPrivate(meshName[0], P); fespace WhIntersectionPrivate(meshIntersection[i], P); ENDIFMACRO - matrix meshName#R = interpolate(WhIntersectionPrivate, singleComponentWhPrivate); - meshName#R.thresholding(1.0e-10); - real[int] meshName#C; - int[int] meshName#I; - [meshName#I, intersection[1 + i - j], meshName#C] = meshName#R; - intersection[1 + i - j].resize(meshName#R.nbcoef); + intersection[1 + i - j] = restrict(WhIntersectionPrivate, singleComponentWhPrivate, n2oIntersection[i]); + fespace VhIntersectionPrivate(meshIntersection[i], P1); + int[int] restriction = restrict(VhIntersectionPrivate, VhLocalPrivate, n2oIntersection[i]); + privateDmesh#meshName#intersectionDef[1 + i - j].resize(VhLocalPrivate.ndof); + privateDmesh#meshName#intersectionDef[1 + i - j] = 0.0; + privateDmesh#meshName#intersectionDef[1 + i - j](restriction) = recv[i]; intersection[0][i - j] = intersection[0][i]; } } numberIntersection -= skip.n; intersection[0].resize(numberIntersection); + privateDmesh#meshName#intersectionDef.resize(1 + numberIntersection); + privateDmesh#meshName#intersectionDef[0].resize(numberIntersection); + for[i, value : intersection[0]] privateDmesh#meshName#intersectionDef[0][i] = value; if(fakeInterface != -111111 || overlap != 1) { PhGlobalPrivate suppPartition = khi > 0.1; AddLayers(globalName, suppPartition[], 1, phi[]); @@ -466,7 +530,7 @@ ENDIFMACRO fespace WhPartPrivate(meshName[level - 1], PkPart); WhPartPrivate defPart(func2vec); D[level - 1].resize(WhPartPrivate.ndof); - if((WhPartPrivate.ndof % meshName[level - 1].nt) == 0) { + if((WhPartPrivate.ndof % meshName[level - 1].nt) == 0 && WhPartPrivate.ndof != meshName[level - 1].nv) { int constant = rank; IFMACRO(privateCreateMat) fespace PhLocalPrivate(meshName[level - 1], P0); @@ -478,12 +542,10 @@ IFMACRO(!privateCreateMat) defPart(func2vec) = initPart(abs(part - constant) < 0.1); ENDIFMACRO } - else if(WhPartPrivate.ndof == meshName[level - 1].nv) { + else if(WhPartPrivate.ndof == meshName[level - 1].nv && (WhPartPrivate.ndof % meshName[level - 1].nt) != 0) func2vec[] = khi[]; - } - else { + else defPart(func2vec) = initPart(khi); - } D[level - 1] = func2vec[]; IFMACRO(privateDmesh#PostProcessD) privateDmesh#PostProcessD(D[level - 1]); @@ -643,6 +705,12 @@ NewMacro privateDmesh#CartesianPartitioning() EndMacro IFMACRO(ThName#RefinementFactor) exit(13); /* not currently implemented */ ENDIFMACRO +IFMACRO(removeZeros,true) +exit(103); /* not currently implemented */ +ENDIFMACRO +IFMACRO(trueRestrict,true) +exit(103); /* not currently implemented */ +ENDIFMACRO ENDIFMACRO IFMACRO(!ThName#RefinementFactor) int privateS = 1; @@ -650,9 +718,15 @@ ENDIFMACRO IFMACRO(ThName#RefinementFactor) int privateS = ThName#RefinementFactor; ENDIFMACRO +IFMACRO(ThName#Overlap) +meshN privateBorder; +buildOverlap(ThName, privateBorder, -111111, privateS, ThName#Overlap, intersection, privateDmesh#ThName#khi[0], P1, ThName#Comm, false); +ENDIFMACRO +IFMACRO(!ThName#Overlap) build(ThName, privateS, intersection, privateDmesh#ThName#khi[0], P1, ThName#Comm); ENDIFMACRO ENDIFMACRO +ENDIFMACRO } EndMacro @@ -750,11 +824,23 @@ ENDIFMACRO mpiRequest[int] rqRecv(neighborsName.n); for[i, value : neighborsName] Irecv(processor(value, ThName#Comm, rqRecv[i]), ThTab[i + (i >= between)]); +IFMACRO(!ThName#Overlap) ThTab[between] = trunc(ThName, gamma > 0.01, label = -111112); +ENDIFMACRO +IFMACRO(ThName#Overlap) + ThTab[between] = ThName; +ENDIFMACRO mpiRequest[int] rqSend(neighborsName.n); for[i, value : neighborsName] +IFMACRO(!ThName#Overlap) Isend(processor(value, ThName#Comm, rqSend[i]), ThTab[between]); +ENDIFMACRO +IFMACRO(ThName#Overlap) + Isend(processor(value, ThName#Comm, rqSend[i]), ThName); +ENDIFMACRO +IFMACRO(!ThName#Overlap) ThTab[neighborsName.n + 1] = trunc(ThName, gamma < 0.01, label = -111112); +ENDIFMACRO mpiWaitAll(rqRecv); IFMACRO(dimension,3) ThName = gluemesh(ThTab); @@ -775,7 +861,8 @@ ENDIFMACRO part(m:m + ThTab[i + 1].nt - 1) = neighborsName[i]; m += ThTab[i + 1].nt; } - part(m:ThName.nt - 1) = rank; + if(m != ThName.nt) + part(m:ThName.nt - 1) = rank; mpiWaitAll(rqSend); if(verbosity > 0) { mpiBarrier(ThName#Comm); @@ -792,7 +879,13 @@ ENDIFMACRO IFMACRO(ThName#N2O) NewMacro privateDmesh#N2O()ThName#N2O EndMacro ENDIFMACRO - buildWithPartitioning(ThName, part, 1, intersection, privateDmesh#ThName#khi[0], P1, ThName#Comm) +IFMACRO(ThName#Overlap) + meshN ThBorder; + buildOverlapWithPartitioning(ThName, ThBorder, part, -111111, 1, ThName#Overlap, intersection, privateDmesh#ThName#khi[0], P1, ThName#Comm, false); +ENDIFMACRO +IFMACRO(!ThName#Overlap) + buildWithPartitioning(ThName, part, 1, intersection, privateDmesh#ThName#khi[0], P1, ThName#Comm); +ENDIFMACRO } EndMacro macro initializeDmesh(NewName) @@ -845,7 +938,7 @@ NewMacro ThName#Comm()mpiCommWorld EndMacro ENDIFMACRO IFMACRO(!privateCreateMatCheckDmesh) if(ThName.nv != privateDmesh#ThName#khi[0].n || (privateDmesh#ThName#khi[1].n && ThName.nt != privateDmesh#ThName#khi[1].n)) { - DmeshBuild(ThName) + DmeshCreate(ThName) } ENDIFMACRO NewMacro privateCreateMat()1 EndMacro @@ -883,6 +976,12 @@ ENDIFMACRO ENDIFMACRO IFMACRO(ThName#PostProcessD) NewMacro privateDmesh#PostProcessD(arg)ThName#PostProcessD(arg) EndMacro +IFMACRO(removeZeros,true) + exit(105); /* not currently implemented */ +ENDIFMACRO +IFMACRO(trueRestrict,true) + exit(105); /* not currently implemented */ +ENDIFMACRO ENDIFMACRO partitionPrivate(ThTab, privateCreateMat, privateCreateMat, privateCreateMat, privateCreateMat, privateCreateMat, mpiRank(ThName#Comm), mpiSize(ThName#Comm), 1, 1, 1, privateCreateMat, DTab, PkName, intersection, ThName#Comm, -111111, ThName#PkPart, def, init, 1) } @@ -910,7 +1009,7 @@ ENDIFMACRO } ENDIFMACRO IFMACRO(!privateDmesh#ThName) -DmeshBuild(ThName) +DmeshCreate(ThName) { IFMACRO(!meshN) IFMACRO(dimension,3) @@ -938,7 +1037,7 @@ EndMacro macro createPartition(ThName, PartName, PkName) IFMACRO(!privateDmesh#ThName) -DmeshBuild(ThName) +DmeshCreate(ThName) ENDIFMACRO { NewMacro privateCreateMatCheckDmesh()1 EndMacro @@ -2160,7 +2259,6 @@ ENDIFMACRO cout.scientific << " --- distributed solution scattered (in " << mpiWtime() - timerScatter << ")" << endl; } }// EOM -ENDIFMACRO macro CoherentGlobalMesh(MatT, MatV, ThName, ThGlobal) { IFMACRO(!ThName#Comm) @@ -2216,7 +2314,7 @@ ENDIFMACRO macro DmeshSave(ThName, name)saveDmesh(ThName, name)// EOM macro DmeshLoad(ThName, name)loadDmesh(ThName, name)// EOM -macro DmeshBuild(ThName)buildDmesh(ThName)// EOM +macro DmeshCreate(ThName)buildDmesh(ThName)// EOM macro DmeshReconstruct(ThName)reconstructDmesh(ThName)// EOM macro DmeshReconstructWithNeighbors(ThName, neighborsName)reconstructDmeshWithNeighbors(ThName, neighborsName)// EOM macro DmeshInitialize(NewName)initializeDmesh(NewName)// EOM @@ -2230,3 +2328,4 @@ macro VecGather(ThName, comm, ThGatherName, Pk, uPrivate, uNew)gatherSolution(Th macro VecInterpolate(ThName, Pk, uPrivate, ThNew, PkNew, uNew)transfer(ThName, Pk, uPrivate, ThNew, PkNew, uNew)// EOM macro MatCreate(ThName, MatName, PkName)createMat(ThName, MatName, PkName)// EOM macro MatInterpolate(ThName, Pk, A, ThNew, PkNew, ANew, P)transferMat(ThName, Pk, A, ThNew, PkNew, ANew, P)// EOM +ENDIFMACRO diff --git a/plugin/mpi/Makefile.am b/plugin/mpi/Makefile.am index be5f9726c..1d4054809 100644 --- a/plugin/mpi/Makefile.am +++ b/plugin/mpi/Makefile.am @@ -25,8 +25,8 @@ if FFCS_MPIOK # compile nothing if no mpi ... LIST_COMPILE=@TOOL_DYLIB_mumps@ MPICG.$(DYLIB_SUFFIX) \ -mpi-cmaes.$(DYLIB_SUFFIX) @TOOL_DYLIB_parmetis@ @TOOL_DYLIB_parmmg@ @TOOL_DYLIB_hpddm@ @TOOL_DYLIB_bem@ -# RadiativeTransfer_htool is inclus in @TOOL_DYLIB_bem@ +mpi-cmaes.$(DYLIB_SUFFIX) @TOOL_DYLIB_parmetis@ @TOOL_DYLIB_parmmg@ @TOOL_DYLIB_hpddm@ @TOOL_DYLIB_bemtool@ +# RadiativeTransfer_htool is inclus in @TOOL_DYLIB_bemtool@ # FFCS - do not set loadpath to be able to run an external version of FF on the examples in this directory with # [[../../mkffref]] diff --git a/plugin/mpi/PETSc-code.hpp b/plugin/mpi/PETSc-code.hpp index d6a684170..4bffd4f29 100644 --- a/plugin/mpi/PETSc-code.hpp +++ b/plugin/mpi/PETSc-code.hpp @@ -302,7 +302,7 @@ void varfBem(const typename fes1::FESpace*& PUh, const typename fes2::FESpace*& typename Mesh2::RdHat pbt(1./(Mesh2::RdHat::d+1),1./(Mesh2::RdHat::d+1)); p2.reserve(3*m); - int mDofScalar = m/nnn; // computation of the dof of one component + int mDofScalar = m/nnn; // computation of the dof of one component for (int i=0; i(*dynamic_cast*>(y))), a(x) { assert(b && b->nargs); ffassert(FieldOfForm(b->largs, IsComplexType>::value) == IsComplexType>::value); - + #if defined(WITH_bemtool) && defined(WITH_htool) && defined(PETSC_HAVE_HTOOL) // Check the nbitem of inconnu and test in BemFormBilinear checkNbItemFEspacesInconnuAndTest(b->largs,b->N,b->M); @@ -1810,6 +1810,12 @@ namespace PETSc { } else if (r == atype< OpTrans< DistributedCSR< HpddmType > > >( )) { e_M[i][j] = e; t_M[i][j] = 2; + } else if (atype< KNM< PetscScalar >* >( )->CastingFrom(r)) { + e_M[i][j] = to< KNM< PetscScalar >* >(c_M); + t_M[i][j] = 9; + } else if (atype< Transpose< KNM< PetscScalar >* > >( )->CastingFrom(r)) { + e_M[i][j] = to< Transpose< KNM< PetscScalar >* > >(c_M); + t_M[i][j] = 10; } else if (atype< KN_< PetscScalar > >( )->CastingFrom(r)) { e_M[i][j] = to< KN_< PetscScalar > >(c_M); t_M[i][j] = 3; @@ -1843,11 +1849,7 @@ namespace PETSc { static E_F0* f(const basicAC_F0& args) { return new initCSRfromBlockMatrix(args, 0); } AnyType operator( )(Stack s) const { Dmat** exchange = new Dmat*[N * M](); - Mat* a = new Mat[N * M]( ); - std::vector< int > zeros; - zeros.reserve(std::min(N, M)); - std::vector> destroy; - destroy.reserve(N * M); + Mat* a = new Mat[N * M](); MPI_Comm comm1 = MPI_COMM_NULL, comm2 = MPI_COMM_NULL; for (int i = 0; i < N; ++i) { for (int j = 0; j < M; ++j) { @@ -1866,12 +1868,10 @@ namespace PETSc { } comm1 = comm2; exchange[i * M + j] = pt; + PetscObjectReference((PetscObject)a[i * M + j]); } else if (t == 2) { DistributedCSR< HpddmType >* pt = GetAny< DistributedCSR< HpddmType >* >(e_ij); - Mat B; - MatCreateHermitianTranspose(pt->_petsc, &B); - a[i * M + j] = B; - destroy.emplace_back(i, j); + MatCreateHermitianTranspose(pt->_petsc, a + i * M + j); exchange[i * M + j] = pt; PetscObjectGetComm((PetscObject)pt->_petsc, &comm2); if (comm1 != MPI_COMM_NULL) { @@ -1883,95 +1883,108 @@ namespace PETSc { } else if (t == 7) { PetscScalar r = GetAny< PetscScalar >(e_ij); if (std::abs(r) > 1.0e-16) { - Mat B; - MatCreateDense(comm1 != MPI_COMM_NULL ? comm1 : PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &B); + MatCreateDense(comm1 != MPI_COMM_NULL ? comm1 : PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, a + i * M + j); PetscScalar* array; - MatDenseGetArray(B, &array); + MatDenseGetArray(a[i * M + j], &array); if (array) array[0] = r; - MatDenseRestoreArray(B, &array); - a[i * M + j] = B; - destroy.emplace_back(i, j); - } else if (i == j) - zeros.emplace_back(i); - } else if (t == 3 || t == 4) { - KN< PetscScalar > x; - Mat B; - if (t == 3) x = GetAny< KN_< PetscScalar > >(e_ij); - else x = GetAny< Transpose< KN_< PetscScalar > > >(e_ij); - MatCreateDense(comm1 != MPI_COMM_NULL ? comm1 : PETSC_COMM_WORLD, x.n, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, &B); - PetscScalar* array; - MatDenseGetArray(B, &array); - if (array) std::copy_n(static_cast< PetscScalar* >(x), x.n, array); - MatDenseRestoreArray(B, &array); - if (t == 3) { - a[i * M + j] = B; + MatDenseRestoreArray(a[i * M + j], &array); } - else { - Mat C; - MatCreateHermitianTranspose(B, &C); - MatDestroy(&B); - a[i * M + j] = C; + } else if (t == 9 || t == 10) { + if (!a[i * M + j]) { + KNM< PetscScalar >* x; + if (t == 9) x = GetAny< KNM< PetscScalar >* >(e_ij); + else x = GetAny< Transpose< KNM< PetscScalar >* > >(e_ij); + PetscMPIInt rank; + MPI_Comm_rank(comm1 != MPI_COMM_NULL ? comm1 : PETSC_COMM_WORLD, &rank); + MatCreateDense(comm1 != MPI_COMM_NULL ? comm1 : PETSC_COMM_WORLD, x->N(), rank == 0 ? x->M() : 0, PETSC_DECIDE, x->M(), NULL, a + i * M + j); + PetscScalar* array; + MatDenseGetArray(a[i * M + j], &array); + if (array) std::copy_n(static_cast< PetscScalar* >(*x), x->N() * x->M(), array); + MatDenseRestoreArray(a[i * M + j], &array); + if (t == 10) { + Mat C; + MatCreateHermitianTranspose(a[i * M + j], &C); + MatDestroy(a + i * M + j); + a[i * M + j] = C; + } + if (i < M && j < N) { + Expression e = e_M[j][i]; + int u = t_M[j][i]; + if (e && ((t == 9 && u == 10) || (t == 10 && u == 9))) { + AnyType e_ji = (*e)(s); + KNM< PetscScalar >* y; + if (u == 9) y = GetAny< KNM< PetscScalar >* >(e_ji); + else y = GetAny< Transpose< KNM< PetscScalar >* > >(e_ji); + if (y == x) { + if (u == 9) { + MatHermitianTransposeGetMat(a[i * M + j], a + j * N + i); + PetscObjectReference((PetscObject)a[j * N + i]); + } + else + MatCreateHermitianTranspose(a[i * M + j], a + j * N + i); + } + } + } + } + } else if (t == 3 || t == 4) { + if (!a[i * M + j]) { + KN< PetscScalar > x; + if (t == 3) x = GetAny< KN_< PetscScalar > >(e_ij); + else x = GetAny< Transpose< KN_< PetscScalar > > >(e_ij); + MatCreateDense(comm1 != MPI_COMM_NULL ? comm1 : PETSC_COMM_WORLD, x.n, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, a + i * M + j); + PetscScalar* array; + MatDenseGetArray(a[i * M + j], &array); + if (array) std::copy_n(static_cast< PetscScalar* >(x), x.n, array); + MatDenseRestoreArray(a[i * M + j], &array); + if (t == 4) { + Mat C; + MatCreateHermitianTranspose(a[i * M + j], &C); + MatDestroy(a + i * M + j); + a[i * M + j] = C; + } + if (i < M && j < N) { + Expression e = e_M[j][i]; + int u = t_M[j][i]; + if (e && ((t == 3 && u == 4) || (t == 4 && u == 3))) { + AnyType e_ji = (*e)(s); + KN< PetscScalar > y; + if (u == 3) y = GetAny< KN_< PetscScalar > >(e_ji); + else y = GetAny< Transpose< KN_< PetscScalar > > >(e_ji); + y -= x; + double norm = y.linfty(); + MPI_Allreduce(MPI_IN_PLACE, &norm, 1, MPI_DOUBLE, MPI_MAX, comm1 != MPI_COMM_NULL ? comm1 : PETSC_COMM_WORLD); + if (norm < 1.0E-14) { + if (u == 3) { + MatHermitianTransposeGetMat(a[i * M + j], a + j * N + i); + PetscObjectReference((PetscObject)a[j * N + i]); + } + else + MatCreateHermitianTranspose(a[i * M + j], a + j * N + i); + } + } + } } - destroy.emplace_back(i, j); } else if (t == 8) { std::pair* p = GetAny*>(e_ij); ffassert(p && p->first->_petsc && p->second->_petsc); Mat mats[2] = { p->second->_petsc, p->first->_petsc }; - Mat C; - MatCreateComposite(PetscObjectComm((PetscObject)p->first->_petsc), 2, mats, &C); - MatCompositeSetType(C, MAT_COMPOSITE_MULTIPLICATIVE); - a[i * M + j] = C; + MatCreateComposite(PetscObjectComm((PetscObject)p->first->_petsc), 2, mats, a + i * M + j); + MatCompositeSetType(a[i * M + j], MAT_COMPOSITE_MULTIPLICATIVE); delete p; - destroy.emplace_back(i, j); } else { ExecError("Unknown type in submatrix"); } - } else if (i == j) - zeros.emplace_back(i); - } - } - for (int i = 0; i < zeros.size( ); ++i) { - int posX = -1, posY = -1; - for (int j = 0; j < M && posX == -1; ++j) { - if (j != zeros[i] && a[zeros[i] * M + j]) posX = j; - } - for (int j = 0; j < N && posY == -1; ++j) { - if (j != zeros[i] && a[zeros[i] + j * M]) posY = j; - } - if (posX == -1 && posY == -1) - ExecError("Zero row and zero column"); - else { - PetscInt x, X, y, Y; - if (posX != -1) { - MatGetSize(a[zeros[i] * M + posX], &X, &Y); - MatGetLocalSize(a[zeros[i] * M + posX], &x, &y); - } - if (posY != -1) { - MatGetSize(a[zeros[i] + posY * M], posX == -1 ? &X : NULL, &Y); - MatGetLocalSize(a[zeros[i] + posY * M], posX == -1 ? &x : NULL, &y); - if (posX == -1) { - x = y; - X = Y; - } - } else { - y = x; - Y = X; } - MatCreate(comm1 != MPI_COMM_NULL ? comm1 : PETSC_COMM_WORLD, a + zeros[i] * M + zeros[i]); - MatSetSizes(a[zeros[i] * M + zeros[i]], x, y, X, Y); - MatSetType(a[zeros[i] * M + zeros[i]], MATAIJ); - MatSeqAIJSetPreallocation(a[zeros[i] * M + zeros[i]], 0, NULL); - MatMPIAIJSetPreallocation(a[zeros[i] * M + zeros[i]], 0, NULL, 0, NULL); - MatAssemblyBegin(a[zeros[i] * M + zeros[i]], MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(a[zeros[i] * M + zeros[i]], MAT_FINAL_ASSEMBLY); - destroy.emplace_back(zeros[i], zeros[i]); } } Result sparse_mat = GetAny< Result >((*emat)(s)); if (sparse_mat->_petsc) sparse_mat->dtor( ); MatCreateNest(comm1 != MPI_COMM_NULL ? comm1 : PETSC_COMM_WORLD, N, NULL, M, NULL, a, &sparse_mat->_petsc); - for(std::pair p : destroy) - MatDestroy(a + p.first * M + p.second); + for (int i = 0; i < N; ++i) + for (int j = 0; j < M; ++j) { + if (a[i * M + j]) + PetscObjectDereference((PetscObject)a[i * M + j]); + } sparse_mat->_exchange = reinterpret_cast**>(exchange); delete[] a; return sparse_mat; @@ -2037,7 +2050,7 @@ namespace PETSc { {"schurPreconditioner", &typeid(KN< Matrice_Creuse< upscaled_type > >*)}, // 5 {"schurList", &typeid(KN< double >*)}, // 6 {"parent", &typeid(Type*)}, // 7 - {"nullspace", &typeid(KNM< upscaled_type >*)}, // 8 + {"nullspace", &typeid(KNM< upscaled_type >*)}, // 8 {"fieldsplit", &typeid(long)}, // 9 {"schurComplement", &typeid(KNM< upscaled_type >*)}, // 10 {"schur", &typeid(KN< Dmat >*)}, // 11 @@ -2390,12 +2403,24 @@ namespace PETSc { PetscInt M, N; MatNestGetSubMats(ptParent->_petsc, &M, &N, &mat); ffassert(M == N); + PCCompositeType type; + PCFieldSplitGetType(pc, &type); + bool keep = false; +#if PETSC_VERSION_GE(3, 20, 0) + if (type == PC_COMPOSITE_SCHUR) { + PCFieldSplitSchurPreType type; + PCFieldSplitGetSchurPre(pc, &type, NULL); + if (type == PC_FIELDSPLIT_SCHUR_PRE_SELF) + keep = true; + } +#endif for (int i = 0; i < M; ++i) { if (mat[i][i] == ptA->_petsc) { KSPDestroy(&ptA->_ksp); ptA->_ksp = subksp[i]; PetscObjectReference((PetscObject)subksp[i]); - KSPSetOperators(subksp[i], ptA->_petsc, ptA->_petsc); + if (!keep) + KSPSetOperators(subksp[i], ptA->_petsc, ptA->_petsc); break; } } @@ -2823,9 +2848,16 @@ namespace PETSc { }; template< class Type, unsigned short O, typename std::enable_if>::value>::type* = nullptr > AnyType view_dispatched(Type* const& ptA, std::string const& o, std::string* const& type, std::string* const& name, MPI_Comm const& comm) { + bool pop = false; + PetscViewer viewer = NULL; + if(type && type->compare("info_detail") == 0) { + ffassert(!O && (ptA->_petsc || ptA->_ksp)); + if(name) PetscViewerASCIIOpen(!O ? PetscObjectComm(ptA->_petsc ? (PetscObject)ptA->_petsc : (PetscObject)ptA->_ksp) : comm, name->c_str(), &viewer); + else viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)viewer)); + PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO_DETAIL); + pop = true; + } if (o.size() == 0 || o.compare("MAT") == 0) { - bool pop = false; - PetscViewer viewer = NULL; if(type) { if(type->compare("matlab") == 0) { if(name) PetscViewerASCIIOpen(!O ? PetscObjectComm((PetscObject)ptA->_petsc) : comm, name->c_str(), &viewer); @@ -2861,10 +2893,6 @@ namespace PETSc { MatCreate(comm, &ptA->_petsc); MatLoad(ptA->_petsc, viewer); } - if(pop) - PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ptA->_petsc))); - if(name) - PetscViewerDestroy(&viewer); } else { ffassert(!O); @@ -2874,10 +2902,14 @@ namespace PETSc { else if (o.compare("PC") == 0) { PC pc; KSPGetPC(ptA->_ksp, &pc); - PCView(pc, PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ptA->_ksp))); + PCView(pc, viewer ? viewer : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ptA->_ksp))); } } } + if(pop) + PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ptA->_petsc))); + if(name) + PetscViewerDestroy(&viewer); return 0L; } template< class Type, unsigned short O, typename std::enable_if>::value>::type* = nullptr > @@ -4708,12 +4740,16 @@ namespace PETSc { type = MATHERMITIANTRANSPOSEVIRTUAL; } } else MatGetSize(mat[i][j], &n, &m); - if(isType && (m > 1 || n == 1)) { - if(mpirank == 0) { - if(U) - *ptr++ = *out++; - else - *out++ = *ptr++; + if(isType && m >= n) { + PetscMPIInt rank; + MPI_Comm_rank(PetscObjectComm((PetscObject)mat[i][j]), &rank); + if(rank == 0) { + for(int m = 0; m < n; ++m) { + if(U) + *ptr++ = *out++; + else + *out++ = *ptr++; + } } break; } @@ -4779,8 +4815,7 @@ namespace PETSc { MatGetType((*t)._petsc, &type); PetscStrcmp(type, MATNEST, &isType); PetscScalar* p = reinterpret_cast(u->operator upscaled_type*()); - - if( (*t)._vector_global ){ + if ((*t)._vector_global) { for(int i = 0; i < u->n; ++i) p[i] = u->operator[](i); MPI_Allreduce(MPI_IN_PLACE, p, u->n, HPDDM::Wrapper::mpi_type(), MPI_SUM, PetscObjectComm((PetscObject)(*t)._petsc)); @@ -4790,9 +4825,8 @@ namespace PETSc { PetscInt M, N; PetscInt Mi = 0, Ni = 0; MatNestGetSubMats((*t)._petsc, &M, &N, &mat); - - PetscInt offset_cols=0; - PetscInt offset_local_cols=0; + PetscInt offset_cols = 0; + PetscInt offset_local_cols = 0; for(PetscInt j = 0; j < N; ++j) { // cols for(PetscInt i = 0; i < M; ++i) { // rows Ni = 0, Mi = 0; @@ -4800,29 +4834,27 @@ namespace PETSc { MatGetSize(mat[i][j], &Mi, &Ni); PetscInt cbegin; - PetscInt ni,mi; - MatGetLocalSize( mat[i][j], &mi, &ni); - MatGetOwnershipRangeColumn( mat[i][j], &cbegin, NULL); + PetscInt ni, mi; + MatGetLocalSize(mat[i][j], &mi, &ni); + MatGetOwnershipRangeColumn(mat[i][j], &cbegin, NULL); for(PetscInt i = 0; i < ni; ++i) - ptr[i+offset_local_cols] = p[i+cbegin+offset_cols]; + ptr[i + offset_local_cols] = p[i + cbegin + offset_cols]; offset_cols += Ni; offset_local_cols += ni; break; } // if mat[i][j] - } // loop rows } // loop cols } // if nested matrix - else{ + else { // case not nested matrix - PetscInt cbegin; - PetscInt n,m,N,M; + PetscInt cbegin, n, m, N, M; MatGetLocalSize( (*t)._petsc, &m, &n); MatGetOwnershipRangeColumn( (*t)._petsc, &cbegin, NULL); for(PetscInt i = 0; i < n; ++i) - ptr[i] = p[i+cbegin]; + ptr[i] = p[i + cbegin]; } VecRestoreArray(x, &ptr); } @@ -4857,7 +4889,7 @@ namespace PETSc { VecRestoreArray(y, &ptr); } } - if (t.A->_A) std::fill_n(out->operator upscaled_type*(), out->n, 0.0); + *out = K(); } else { VecSet(x, PetscScalar( )); Vec isVec; @@ -4897,7 +4929,7 @@ namespace PETSc { if (verbosity > 0 && mpirank == 0) cout << " --- system solved with PETSc (in " << MPI_Wtime( ) - timing << ")" << endl; p = reinterpret_cast(out->operator upscaled_type*()); - if((*t)._vector_global){ + if ((*t)._vector_global) { VecGetArray(y, &ptr); if(isType) { Mat** mat; @@ -4905,43 +4937,38 @@ namespace PETSc { PetscInt Mi = 0, Ni = 0; MatNestGetSubMats((*t)._petsc, &M, &N, &mat); - PetscInt offset_rows=0; - PetscInt offset_local_rows=0; + PetscInt offset_rows = 0; + PetscInt offset_local_rows = 0; for(PetscInt i = 0; i < M; ++i) { // rows for(PetscInt j = 0; j < N; ++j) { // cols - Ni = 0, Mi = 0; if(mat[i][j]) { MatGetSize(mat[i][j], &Mi, &Ni); PetscInt rbegin; PetscInt ni,mi; - MatGetLocalSize( mat[i][j], &mi, &ni); - MatGetOwnershipRange( mat[i][j], &rbegin, NULL); + MatGetLocalSize(mat[i][j], &mi, &ni); + MatGetOwnershipRange(mat[i][j], &rbegin, NULL); for(PetscInt i = 0; i < mi; ++i) - p[i+rbegin+offset_rows] = ptr[i+offset_local_rows]; + p[i + rbegin + offset_rows] = ptr[i + offset_local_rows]; offset_rows += Mi; offset_local_rows += mi; break; } // if mat[i][j] - } // loop rows } // loop cols } // if nested matrix - else{ + else { // case not nested matrix - PetscInt rbegin; - PetscInt n,m,N,M; - MatGetOwnershipRange( (*t)._petsc, &rbegin, NULL); - MatGetLocalSize( (*t)._petsc, &m, &n); - std::fill(p,p+u->n,0.); - + PetscInt rbegin, n, m, N, M; + MatGetOwnershipRange((*t)._petsc, &rbegin, NULL); + MatGetLocalSize((*t)._petsc, &m, &n); + std::fill_n(p, u->n, 0.0); for(PetscInt i = 0; i < m; ++i) - p[i+rbegin] = ptr[i]; + p[i + rbegin] = ptr[i]; } - VecRestoreArray(y, &ptr); MPI_Allreduce(MPI_IN_PLACE, p, u->n, HPDDM::Wrapper::mpi_type(), MPI_SUM, PetscObjectComm((PetscObject)(*t)._petsc)); } @@ -5015,8 +5042,7 @@ namespace PETSc { PetscBool isType; PetscScalar* p = reinterpret_cast(u->operator upscaled_type*()); PetscStrcmp(type, MATNEST, &isType); - - if((*t)._vector_global){ + if((*t)._vector_global) { MPI_Allreduce(MPI_IN_PLACE, p, u->n, HPDDM::Wrapper::mpi_type(), MPI_SUM, PetscObjectComm((PetscObject)(*t)._petsc)); // VecGetArray(x, &ptr) is defined before if(isType) { // case nested matrix @@ -5025,8 +5051,8 @@ namespace PETSc { PetscInt Mi = 0, Ni = 0; MatNestGetSubMats((*t)._petsc, &Mb, &Nb, &mat); - PetscInt offset_cols=0; - PetscInt offset_local_cols=0; + PetscInt offset_cols = 0; + PetscInt offset_local_cols = 0; for(PetscInt j = 0; j < Nb; ++j) { // cols for(PetscInt i = 0; i < Mb; ++i) { // rows Ni = 0, Mi = 0; @@ -5034,29 +5060,27 @@ namespace PETSc { MatGetSize(mat[i][j], &Mi, &Ni); PetscInt cbegin; - PetscInt ni,mi; - MatGetLocalSize( mat[i][j], &mi, &ni); - MatGetOwnershipRangeColumn( mat[i][j], &cbegin, NULL); + PetscInt ni, mi; + MatGetLocalSize(mat[i][j], &mi, &ni); + MatGetOwnershipRangeColumn(mat[i][j], &cbegin, NULL); for(PetscInt i = 0; i < ni; ++i) - ptr[i+offset_local_cols] = p[i+cbegin+offset_cols]; + ptr[i + offset_local_cols] = p[i + cbegin + offset_cols]; offset_cols += Ni; offset_local_cols += ni; break; } // if mat[i][j] - } // loop rows } // loop cols } // if nested matrix - else{ + else { // case not nested matrix - PetscInt cbegin; - PetscInt n,m; - MatGetLocalSize( (*t)._petsc, &m, &n); - MatGetOwnershipRangeColumn( (*t)._petsc, &cbegin, NULL); + PetscInt cbegin, n,m; + MatGetLocalSize((*t)._petsc, &m, &n); + MatGetOwnershipRangeColumn((*t)._petsc, &cbegin, NULL); for(PetscInt i = 0; i < n; ++i) - ptr[i] = p[i+cbegin]; + ptr[i] = p[i + cbegin]; } } else if (isType) { @@ -5098,8 +5122,8 @@ namespace PETSc { PetscInt Mi = 0, Ni = 0; MatNestGetSubMats((*t)._petsc, &Mb, &Nb, &mat); - PetscInt offset_rows=0; - PetscInt offset_local_rows=0; + PetscInt offset_rows = 0; + PetscInt offset_local_rows = 0; for(PetscInt i = 0; i < Mb; ++i) { // rows for(PetscInt j = 0; j < Nb; ++j) { // cols @@ -5107,32 +5131,28 @@ namespace PETSc { if(mat[i][j]) { MatGetSize(mat[i][j], &Mi, &Ni); - PetscInt rbegin; - PetscInt ni,mi; - MatGetLocalSize( mat[i][j], &mi, &ni); - MatGetOwnershipRange( mat[i][j], &rbegin, NULL); + PetscInt rbegin, ni, mi; + MatGetLocalSize(mat[i][j], &mi, &ni); + MatGetOwnershipRange(mat[i][j], &rbegin, NULL); for(PetscInt i = 0; i < mi; ++i) - p[i+rbegin+offset_rows] = ptr[i+offset_local_rows]; + p[i + rbegin + offset_rows] = ptr[i + offset_local_rows]; offset_rows += Mi; offset_local_rows += mi; break; } // if mat[i][j] - } // loop rows } // loop cols } // if nested matrix else{ // case not nested matrix - PetscInt rbegin; - PetscInt n,m; - MatGetOwnershipRange( (*t)._petsc, &rbegin, NULL); - MatGetLocalSize( (*t)._petsc, &m, &n); - std::fill(p,p+u->n,0.); - + PetscInt rbegin, n, m; + MatGetOwnershipRange((*t)._petsc, &rbegin, NULL); + MatGetLocalSize((*t)._petsc, &m, &n); + std::fill_n(p, u->n, 0.0); for(PetscInt i = 0; i < m; ++i) - p[i+rbegin] = ptr[i]; + p[i + rbegin] = ptr[i]; } VecRestoreArray(y, &ptr); @@ -5174,7 +5194,7 @@ namespace PETSc { u->operator[](i) = p[i]; } } - } + } VecRestoreArray(y, &ptr); VecDestroy(&y); } @@ -5863,17 +5883,16 @@ namespace PETSc { { typedef typename Call_CompositeFormBilinear::const_iterator const_iterator; int init; - + class Op : public E_F0mps { public: Call_CompositeFormBilinear *b; Expression a; int init; //AnyType operator()(Stack s) const; - Op(Expression aa,Expression bb,int initt) : b(new Call_CompositeFormBilinear(* dynamic_cast *>(bb))),a(aa),init(initt) - { + { assert(b && b->nargs); int NN = (int) b->euh->componentNbitem().size(); int MM = (int) b->evh->componentNbitem().size(); @@ -5909,36 +5928,36 @@ namespace PETSc { std::set jcols; ff_mat.CSR(); // transform the matrix to CSR format - + std::vector perm_row(ff_mat.n,-1); std::vector perm_col(ff_mat.m,-1); - + for (int ii=0; ii < ff_mat.n; ii++) { for (int la = ff_mat.p[ii]; la < ff_mat.p[ii+1]; la++) { - perm_row[ii] = 1; + perm_row[ii] = 1; perm_col[ff_mat.j[la]] = 1; if( la > ff_mat.p[ii]){ // if( !( ff_mat.j[la] > ff_mat.j[ff_mat.p[ii]]) ){ cerr << " The column index must be croissant :: Error " << ff_mat.j[la] << " " << ff_mat.j[ff_mat.p[ii]] < ff_mat.j[ff_mat.p[ii]] ); + ffassert( ff_mat.j[la] > ff_mat.j[ff_mat.p[ii]] ); } } } - + if(verbosity> 2) for(int ii=0; ii<10; ii++) cout << "perm_row["<((*b->euh)(stack)); pvectgenericfes * pVh= GetAny((*b->evh)(stack)); - ffassert( *pUh && *pVh ); + ffassert( *pUh && *pVh ); // Update is necessary when we get "pvectgenericfes" to take account a new mesh. (*pUh)->update(); (*pVh)->update(); @@ -6043,20 +6062,20 @@ namespace PETSc { KN VhNbItem = (*pVh)->vectOfNbitem(); // - const KNM> & block_largs=b->block_largs; + const KNM> & block_largs=b->block_largs; // check if we have a square matrix bool A_is_square= (void*)pUh == (void*)pVh || ((*pUh)->totalNbOfDF()) == ( (*pVh)->totalNbOfDF()) ; - // === simple check if A is symetrical === // + // === simple check if A is symetrical === // // voir avec les autres. - bool A_is_maybe_sym = (void*)pUh == (void*)pVh; + bool A_is_maybe_sym = (void*)pUh == (void*)pVh; // VF == true => VF type of Matrix //bool VF=isVF(b->block_largs); //=== used to set the solver ??? block matrix ??? ===/ bool VF = 0; - // set parameteer of the matrix :: + // set parameteer of the matrix : Data_Sparse_Solver ds; ds.factorize=0; ds.initmat=true; @@ -6065,14 +6084,14 @@ namespace PETSc { MPI_Comm comm = ds.commworld ? *static_cast< MPI_Comm* >(ds.commworld) : PETSC_COMM_WORLD; - // J'ai repris ce qu'il y avait. + // J'ai repris ce qu'il y avait. // PAC(e) :: Attention peut être pas compatible avec les matrices bloques. // A repenser :: surtout pour le parametre symetrique? on le met ce parametre à zéro pour l'instant. - // set ds.sym = 0 + // set ds.sym = 0 ds.sym = 0; if(verbosity>3) - cout << " === we consider the block matrix as a non symetric matrix === (to be change in the future)" << endl; + cout << " === we consider the block matrix as a non symetric matrix === (to be change in the future)" << endl; if (! A_is_square ) { @@ -6081,7 +6100,7 @@ namespace PETSc { // A quoi cela correspond?? Gestion du stack + autre WhereStackOfPtr2Free(stack)=new StackOfPtr2Free(stack);// FH aout 2007 - + PETSc::DistributedCSR< HpddmType > * Ares( GetAny*>((*a)(stack))); // test function (Vh) are the line @@ -6106,22 +6125,22 @@ namespace PETSc { cout << " " << endl; } // construction du block (i,j) - const list & b_largs=block_largs(i,j); + const list & b_largs=block_largs(i,j); // size of the block int N_block = UhNbOfDf[i]; int M_block = VhNbOfDf[j]; - // initialise the bem block and fem block + // initialise the bem block and fem block Matrice_Creuse Afem; Afem.init(); Afem.resize(M_block,N_block); int nsparseblocks = 0; Mat Abem = nullptr; - if(verbosity>2) cout << "size_block =" << b_largs.size() << endl; + if(verbosity>2) cout << "size_block =" << b_largs.size() << endl; if( b_largs.size()> 0){ - + // compute largs due BEM and FEM part list largs_FEM; list largs_BEM; @@ -6152,102 +6171,138 @@ namespace PETSc { if (VFBEM==1) ffassert (samemesh); - PETSc::DistributedCSR< HpddmType > * Abemblock = new PETSc::DistributedCSR< HpddmType >; + PETSc::DistributedCSR< HpddmType > *Abemblocki = nullptr, *Abemblock = nullptr; // block diagonal matrix - if( (*pUh)->typeFE[i] == 4 && (*pVh)->typeFE[j] == 4 ){ - ffassert( i==j ); // If not a block diagonal not coded yet - // MeshS --- MeshS - // ==== FESpace 3d Surf: inconnue et test === - const FESpaceS * PUh = (FESpaceS *) (*pUh)->vect[i]->getpVh(); - const FESpaceS * PVh = (FESpaceS *) (*pVh)->vect[j]->getpVh(); - - MatCreate(comm, &Abemblock->_petsc); - MatSetSizes(Abemblock->_petsc, PETSC_DECIDE, PETSC_DECIDE, PUh->NbOfDF, PUh->NbOfDF); - varfBem(PUh, PUh, 1, VFBEM, stack, b_largs_zz, ds, Abemblock); - - - } - else if( (*pUh)->typeFE[i] == 5 && (*pVh)->typeFE[j] == 5 ){ - ffassert( i==j ); // If not a block diagonal not coded yet - // MeshL --- MeshL - // ==== FESpace 3d Curve: inconnue et test === - const FESpaceL * PUh = (FESpaceL *) (*pUh)->vect[i]->getpVh(); - const FESpaceL * PVh = (FESpaceL *) (*pVh)->vect[j]->getpVh(); - - MatCreate(comm, &Abemblock->_petsc); - MatSetSizes(Abemblock->_petsc, PETSC_DECIDE, PETSC_DECIDE, PUh->NbOfDF, PUh->NbOfDF); - varfBem(PUh, PUh, 1, VFBEM, stack, b_largs_zz, ds, Abemblock); - - } - else{ - cerr << " BEM bilinear form " << endl; - cerr << " Block ("<< i <<" ,"<< j << ")" << endl; - cerr << " =: Pas prise en compte des FESpace inconnue de type := "<< typeFEtoString( (*pUh)->typeFE[i] ) << endl; - cerr << " =: avec des FESpace test de type := "<< typeFEtoString( (*pVh)->typeFE[j] ) << endl; - ffassert(0); + for (list::const_iterator ii=b_largs_zz.begin(); ii!=b_largs_zz.end(); ii++) { + list bi; + bi.clear(); + bi.push_back(*ii); + + if( (*pUh)->typeFE[i] == 4 && (*pVh)->typeFE[j] == 4 ){ + ffassert( i==j ); // If not a block diagonal not coded yet + // MeshS --- MeshS + // ==== FESpace 3d Surf: inconnue et test === + const FESpaceS * PUh = (FESpaceS *) (*pUh)->vect[i]->getpVh(); + const FESpaceS * PVh = (FESpaceS *) (*pVh)->vect[j]->getpVh(); + + Abemblocki = new PETSc::DistributedCSR< HpddmType >; + MatCreate(comm, &Abemblocki->_petsc); + MatSetSizes(Abemblocki->_petsc, PETSC_DECIDE, PETSC_DECIDE, PUh->NbOfDF, PUh->NbOfDF); + varfBem(PUh, PUh, 1, VFBEM, stack, bi, ds, Abemblocki); + } + else if( (*pUh)->typeFE[i] == 5 && (*pVh)->typeFE[j] == 5 ){ + ffassert( i==j ); // If not a block diagonal not coded yet + // MeshL --- MeshL + // ==== FESpace 3d Curve: inconnue et test === + const FESpaceL * PUh = (FESpaceL *) (*pUh)->vect[i]->getpVh(); + const FESpaceL * PVh = (FESpaceL *) (*pVh)->vect[j]->getpVh(); + + Abemblocki = new PETSc::DistributedCSR< HpddmType >; + MatCreate(comm, &Abemblocki->_petsc); + MatSetSizes(Abemblocki->_petsc, PETSC_DECIDE, PETSC_DECIDE, PUh->NbOfDF, PUh->NbOfDF); + varfBem(PUh, PUh, 1, VFBEM, stack, bi, ds, Abemblocki); + } + else{ + cerr << " BEM bilinear form " << endl; + cerr << " Block ("<< i <<" ,"<< j << ")" << endl; + cerr << " =: Pas prise en compte des FESpace inconnue de type := "<< typeFEtoString( (*pUh)->typeFE[i] ) << endl; + cerr << " =: avec des FESpace test de type := "<< typeFEtoString( (*pVh)->typeFE[j] ) << endl; + ffassert(0); + } + if (Abemblock == nullptr) { + Abemblock = Abemblocki; + } + else { + Mat mats[2] = { Abemblocki->_petsc, Abemblock->_petsc }; + Mat C; + MatCreateComposite(PetscObjectComm((PetscObject)Abemblock->_petsc), 2, mats, &C); + MatCompositeSetType(C, MAT_COMPOSITE_ADDITIVE); + MatDestroy(&Abemblock->_petsc); + Abemblock->_petsc = C; + delete Abemblocki; + } } Abem = Abemblock->_petsc; } else{ bool samemesh = (void*) (*pUh)->vect[i]->getppTh() == (void*) (*pVh)->vect[j]->getppTh(); // same Fem2D::Mesh +++ pot or kernel - - PETSc::DistributedCSR< HpddmType > * Abemblock = new PETSc::DistributedCSR< HpddmType >; - // block non diagonal matrix - if( (*pUh)->typeFE[i] == 5 && (*pVh)->typeFE[j] == 2 ){ - // case Uh[i] == MeshL et Vh[j] = Mesh2 - - if(verbosity >3) cout << " === creation de la matrice BEM pour un bloc non diagonal === " << endl; - const FESpaceL * PUh = (FESpaceL *) (*pUh)->vect[i]->getpVh(); + PETSc::DistributedCSR< HpddmType > *Abemblocki = nullptr, *Abemblock = nullptr; - MatCreate(comm, &Abemblock->_petsc); - MatSetSizes(Abemblock->_petsc, PETSC_DECIDE, PETSC_DECIDE, PUh->NbOfDF, PUh->NbOfDF); - varfBem(PUh, PUh, 1, VFBEM, stack, b_largs_zz, ds, Abemblock); - } - else if( (*pUh)->typeFE[i] == 2 && (*pVh)->typeFE[j] == 5 ){ - // case Uh[i] == Mesh2 et Vh[j] = MeshL + // block non diagonal matrix + for (list::const_iterator ii=b_largs_zz.begin(); ii!=b_largs_zz.end(); ii++) { + list bi; + bi.clear(); + bi.push_back(*ii); - if(verbosity >3) cout << " === creation de la matrice BEM pour un bloc non diagonal === " << endl; - const FESpaceL * PVh = (FESpaceL *) (*pVh)->vect[j]->getpVh(); + if( (*pUh)->typeFE[i] == 5 && (*pVh)->typeFE[j] == 2 ){ + // case Uh[i] == MeshL et Vh[j] = Mesh2 - MatCreate(comm, &Abemblock->_petsc); - MatSetSizes(Abemblock->_petsc, PETSC_DECIDE, PETSC_DECIDE, PVh->NbOfDF, PVh->NbOfDF); - varfBem(PVh, PVh, 1, VFBEM, stack, b_largs_zz, ds, Abemblock); - } - else{ - cerr << " BEM bilinear form " << endl; - cerr << " Block ("<< i <<" ,"<< j << ")" << endl; - cerr << " =: Pas prise en compte des FESpace inconnue de type := "<< typeFEtoString( (*pUh)->typeFE[i] ) << endl; - cerr << " =: avec des FESpace test de type := "<< typeFEtoString( (*pVh)->typeFE[j] ) << endl; - ffassert(0); + if(verbosity >3) cout << " === creation de la matrice BEM pour un bloc non diagonal === " << endl; + const FESpaceL * PUh = (FESpaceL *) (*pUh)->vect[i]->getpVh(); + + Abemblocki = new PETSc::DistributedCSR< HpddmType >; + MatCreate(comm, &Abemblocki->_petsc); + MatSetSizes(Abemblocki->_petsc, PETSC_DECIDE, PETSC_DECIDE, PUh->NbOfDF, PUh->NbOfDF); + varfBem(PUh, PUh, 1, VFBEM, stack, bi, ds, Abemblocki); + } + else if( (*pUh)->typeFE[i] == 2 && (*pVh)->typeFE[j] == 5 ){ + // case Uh[i] == Mesh2 et Vh[j] = MeshL + + if(verbosity >3) cout << " === creation de la matrice BEM pour un bloc non diagonal === " << endl; + const FESpaceL * PVh = (FESpaceL *) (*pVh)->vect[j]->getpVh(); + + Abemblocki = new PETSc::DistributedCSR< HpddmType >; + MatCreate(comm, &Abemblocki->_petsc); + MatSetSizes(Abemblocki->_petsc, PETSC_DECIDE, PETSC_DECIDE, PVh->NbOfDF, PVh->NbOfDF); + varfBem(PVh, PVh, 1, VFBEM, stack, bi, ds, Abemblocki); + } + else{ + cerr << " BEM bilinear form " << endl; + cerr << " Block ("<< i <<" ,"<< j << ")" << endl; + cerr << " =: Pas prise en compte des FESpace inconnue de type := "<< typeFEtoString( (*pUh)->typeFE[i] ) << endl; + cerr << " =: avec des FESpace test de type := "<< typeFEtoString( (*pVh)->typeFE[j] ) << endl; + ffassert(0); + } + if (Abemblock == nullptr) { + Abemblock = Abemblocki; + } + else { + Mat mats[2] = { Abemblocki->_petsc, Abemblock->_petsc }; + Mat C; + MatCreateComposite(PetscObjectComm((PetscObject)Abemblock->_petsc), 2, mats, &C); + MatCompositeSetType(C, MAT_COMPOSITE_ADDITIVE); + MatDestroy(&Abemblock->_petsc); + Abemblock->_petsc = C; + delete Abemblocki; + } } - - // creation de la matrice dense - + // creation de la matrice dense + // BEM matrix is constructed with different FESpace ffassert( (*pUh)->vect[i]->getpVh() != (*pVh)->vect[j]->getpVh() ) ; - + if( (*pUh)->typeFE[i] == 5 && (*pVh)->typeFE[j] == 2 ){ - // case Uh[i] == MeshL et Vh[j] = Mesh2 + // case Uh[i] == MeshL et Vh[j] = Mesh2 const FESpaceL * PUh = (FESpaceL *) (*pUh)->vect[i]->getpVh(); const FESpace * PVh = (FESpace *) (*pVh)->vect[j]->getpVh(); // construction of the matrix of interpolation - + // The transpose is in the build matrix now bool transpose_MI = false; Matrice_Creuse * MI_BBB = PETSC_buildMatrixInterpolationForCompositeFESpace( PUh, PVh, transpose_MI ); MatriceMorse * mr =MI_BBB->pHM(); - - // Transform Real Matrix in Complex Matrix + + // Transform Real Matrix in Complex Matrix MatriceMorse * mA = new MatriceMorse(mr->n,mr->m,0,0); // we divide by mpisize because the interpolation matrix is computed in sequential int size; MPI_Comm_size(comm, &size); *mr *= 1.0/size; *mA = *mr; - + // transform intrerpolation matrix to matrix IS Mat matIS_MI; ff_createMatIS( *mA, matIS_MI, comm); @@ -6257,8 +6312,8 @@ namespace PETSc { if (std::is_same< PetscScalar, PetscReal >::value) MatCreateTranspose(matIS_MI, &mAAT); else MatCreateHermitianTranspose(matIS_MI, &mAAT); MatDestroy(&matIS_MI); - - // create composite + + // create composite Mat mats[2] = { Abemblock->_petsc , mAAT}; Mat C; MatCreateComposite(PetscObjectComm((PetscObject)Abemblock->_petsc), 2, mats, &C); @@ -6309,19 +6364,19 @@ namespace PETSc { ffassert(0); } - + } #endif } if( largs_FEM.size() >0){ const list & b_largs_zz = largs_FEM; - varfToCompositeBlockLinearSystemALLCASE_pfes( i, j, (*pUh)->typeFE[i], (*pVh)->typeFE[j], + varfToCompositeBlockLinearSystemALLCASE_pfes( i, j, (*pUh)->typeFE[i], (*pVh)->typeFE[j], 0, 0, (*pUh)->vect[i], (*pVh)->vect[j], true, false, ds.sym, ds.tgv, ds.commworld, - b_largs_zz, stack, + b_largs_zz, stack, 0, 0, Afem.pHM()); - + nsparseblocks++; } } @@ -6329,13 +6384,13 @@ namespace PETSc { if(b_largs.size()> 0){ Afem.pHM()->half = ds.sym; - + Mat matIS; if(nsparseblocks){ ff_createMatIS( *(Afem.pHM()), matIS, comm); Afem.destroy(); } - + if (Abem != nullptr) { if (!nsparseblocks) { a[j * maxJVh + i] = Abem; @@ -6349,17 +6404,15 @@ namespace PETSc { } } else { - a[j * maxJVh + i] = matIS; + a[j * maxJVh + i] = matIS; } destroy.emplace_back(j, i); } - + } // end loop j offsetMatrixUh += UhNbOfDf[i]; } // end loop i - - // - if(verbosity>3) cout << "Ares->_vector_global=" << Ares->_vector_global << endl; + if(verbosity>3) cout << "Ares->_vector_global=" << Ares->_vector_global << endl; if( NpUh==1 && maxJVh==1 ){ Ares->_petsc = a[0]; diff --git a/plugin/mpi/common_hpddm.hpp b/plugin/mpi/common_hpddm.hpp index 4ae61286a..7c8c4a622 100644 --- a/plugin/mpi/common_hpddm.hpp +++ b/plugin/mpi/common_hpddm.hpp @@ -592,6 +592,7 @@ static void Init_Common() { if(!Global.Find("savevtk").NotNull()) { Global.Add("savevtk", "(", new OneOperatorCode ); Global.Add("savevtk", "(", new OneOperatorCode ); + Global.Add("savevtk", "(", new OneOperatorCode< VTK_WriteMeshT_Op< MeshS > >); } if(!Global.Find("periodicity").NotNull()) { Global.Add("periodicity", "(", new OneOperator3_*, KN< KN< long > >*, KN< double >*>(periodicity)); diff --git a/plugin/seq/Element_Mixte.cpp b/plugin/seq/Element_Mixte.cpp index 8b00fa71c..fa1739f47 100644 --- a/plugin/seq/Element_Mixte.cpp +++ b/plugin/seq/Element_Mixte.cpp @@ -1197,6 +1197,265 @@ void TypeOfFE_RTdc::Pi_h_alpha(const baseFElement & K,KN_ & v) const v[k++]=-signe*E.x; } } + +// Add jun 2023 FH... in test +class TypeOfFE_BDM2_2d : public TypeOfFE { +public: + static int Data[]; + static double Pi_h_coef[]; + bool Ortho; + const QuadratureFormular1d &QFE; + const QuadratureFormular &QFT; + TypeOfFE_BDM2_2d(bool ortho) + : TypeOfFE(12, 2, Data, 1, 1, + 6*(3*3+7) + , // nb coef mat interpole + 3*3+7, // nb P interpolation + 0), + QFE(QF_GaussLegendre3),QFT(QuadratureFormular_T_5), Ortho(ortho) { + int kkk = 0, i = 0; + // point on edge + for (int e = 0; e < 3; ++e) { + for (int p = 0; p < QFE.n; ++p) { + R2 A(TriangleHat[VerticesOfTriangularEdge[e][0]]); + R2 B(TriangleHat[VerticesOfTriangularEdge[e][1]]); + // 3 dof per edge and 2 componante + pij_alpha[kkk++] = IPJ(3 * e, i, 0); + pij_alpha[kkk++] = IPJ(3 * e, i, 1); + pij_alpha[kkk++] = IPJ(3 * e + 1, i, 0); + pij_alpha[kkk++] = IPJ(3 * e + 1, i, 1); + pij_alpha[kkk++] = IPJ(3 * e + 2, i, 0); + pij_alpha[kkk++] = IPJ(3 * e + 2, i, 1); + + P_Pi_h[i++] = B * (QFE[p].x) + A * (1. - QFE[p].x); // X=0 => A X=1 => B; + } + } + // kkk = 3*QFE.n*6 (QFE.n == 5) + // point on triangle (QFT.n == 7) + // kkk += QFT.n * 6 + for (int p = 0; p < QFT.n; ++p) { + // 3 dof per edge and 2 componante + pij_alpha[kkk++] = IPJ(9, i, 0); + pij_alpha[kkk++] = IPJ(9, i, 1); + pij_alpha[kkk++] = IPJ(9 + 1, i, 0); + pij_alpha[kkk++] = IPJ(9 + 1, i, 1); + pij_alpha[kkk++] = IPJ(9 + 2, i, 0); + pij_alpha[kkk++] = IPJ(9 + 2, i, 1); + + P_Pi_h[i++] = QFT[p]; // X=0 => A X=1 => B; + } + if(verbosity>99) + { + cout << "\n kkk " << kkk << " " << this->pij_alpha.N( ) << " /"<< QFE.n << " " << QFT.n <P_Pi_h.N( ) << endl; + } + ffassert(kkk == this->pij_alpha.N( )); + ffassert(i == this->P_Pi_h.N( )); + } + + void Pi_h_alpha(const baseFElement &K, + KN_< double > &v) const + { // compute the coef of interpolation ... + const Triangle &T(K.T); + int k = 0; + double c1[][3] = {{9, -18, 3} /* 0 */, {-18, 84, -18} /* 1 */, {3, -18, 9} /* 2 */}; + // QuadratureFormular1d QFe=QF_GaussLegendre3; + for (int i = 0; i < 3; i++) { + R2 E(Ortho ? T.Edge(i) : -T.Edge(i).perp( )); + R s = T.EdgeOrientation(i); + + for (int p = 0; p < QFE.n; ++p) { + R l1 = QFE[p].x, l2 = 1 - QFE[p].x; + + R p0 = l1 ; + R p1 = l1 * l2; + R p2 = l2 ; + R sa = QFE[p].a*s; + R cc2 = sa * p0; // + R cc1 = sa * p1; // + R cc0 = sa * p2; // + if (s < 0) { + swap(cc0, cc2); + } + + v[k++] = cc0 * E.x; + v[k++] = cc0 * E.y; + v[k++] = cc1 * E.x; + v[k++] = cc1 * E.y; + v[k++] = cc2 * E.x; + v[k++] = cc2 * E.y; + } + } + // cout << " " << k << " == " << this->pij_alpha.N( ) << " " << QFE.n << endl; + // maintenant les 3 dof interne ... + for (int p = 0; p < QFT.n; ++p) { + R2 Q=T(QFT[p]); + for(int e=0; e<3;++e) + { + R a = QFT[p].a; + R2 P= (Q-(R2) T[e]).perp()/2; + if(Ortho) + P=P.perp( );// correction 4 sep 2023 FH.. + v[k++] = a * P.x ; + v[k++] = a * P.y ; + } + } + // cout << " " << k << " == " << this->pij_alpha.N( ) << " " << QFT.n <pij_alpha.N( )); + + } + void FB(const bool *whatd, const Mesh &Th, const Triangle &K, const RdHat &PHat, + RNMK_ &val) const; + +}; + +// ENDOFCLASS TypeOfFE_PkEdge +int TypeOfFE_BDM2_2d::Data[] = {3, 3, 3, 4,4,4, 5,5,5, 6,6,6, // support on what + 0, 1, 2, 0,1,2, 0,1,2, 0,1,2, // df on node + 0, 0, 0, 1,1,1, 2,2,2, 3,3,3, // th node of df + 0, 0, 0, 0,0,0, 0,0,0, 0,0,0, // df previou FE + 0, 1, 2, 3,4,5, 6,7,8, 9,10,11 , // which df on prev + 0,0,0,0, 12, 12}; +void TypeOfFE_BDM2_2d::FB(const bool *whatd, const Mesh &, const Triangle &K, const RdHat &PHat, + RNMK_ &val) const { + R2 X = K(PHat); + R2 Q[] = {R2(K[0]), R2(K[1]), R2(K[2])}; + R det = K.area*2; + + R2 AB(Q[0],Q[1]),AC(Q[0],Q[2]),Piola[2]={R2(AB.x,AC.x)/det,R2(AB.y,AC.y)/det}; + // cout << " Piola : " << Piola[0] << " , " << Piola[1] << " " << det << endl; + R l0 = 1 - PHat.x - PHat.y, l1 = PHat.x, l2 = PHat.y; + R L[6] = {l0, l1, l2, l1*l2, l2*l0, l0*l1}; + R2 Dl[3] = {K.H(0), K.H(1), K.H(2)}; + R2 DL[6] = {Dl[0],Dl[1],Dl[2],l2*Dl[1]+l1*Dl[2], l2*Dl[0]+l0*Dl[2],l0*Dl[1]+l1*Dl[0] }; + R2 O; + R2 DDLx[6] = {O,O,O,Dl[2].x*Dl[1]+Dl[1].x*Dl[2], Dl[2].x*Dl[0]+Dl[0].x*Dl[2],Dl[1].x*Dl[2]+Dl[2].x*Dl[1] }; + R2 DDLy[6] = {O,O,O,Dl[2].y*Dl[1]+Dl[1].y*Dl[2], Dl[2].y*Dl[0]+Dl[0].y*Dl[2],Dl[1].y*Dl[2]+Dl[2].y*Dl[1] }; + + double C1[12][12] = { + {0, 0, 9, 0, 0, 3, -18, -12, 0, 0, -18, 0}, + {0, 0, -30, 0, 0, -30, 90, 90, 0, 30, 30, 0}, + {0, 0, 3, 0, 0, 9, -12, -18, 0, -18, 0, 0}, + + {9, 0, 0, 0, 3, -3, 0, 0, -30, 12, -18, 0}, + {-30, 0, 0, 0, -30, 30, 30, -30, 180, -90, 30, 0}, + {3, 0, 0, 0, 9, -9, -18, 18, -30, 18, 0, 0}, + + {0, -9, 3, -3, 0, 0, 0, 0, 0, 18, -12, 30}, + {0, 30, -30, 30, 0, 0, 30, -30, 0, -30, 90, -180}, + {0, -3, 9, -9, 0, 0, -18, 18, 0, 0, -18, 30}, + + {0, 0, 0, 0, 0, 0, -36, 36, 0, -12, 12, 0}, + {0, 0, 0, 0, 0, 0, -12, 12, 0, -36, 12, 0}, + {0, 0, 0, 0, 0, 0, -12, 12, 0, -12, 36, 0}}; + + + assert(val.N( ) >= 12); + assert(val.M( ) == 2); + + val = 0; + int p[12]={0,1,2,3,4,5,6,7,8,9,10,11}; + double oe[4]={K.EdgeOrientation(0),-K.EdgeOrientation(1),K.EdgeOrientation(2),1.}; + if (oe[0] < 0) swap(p[0],p[2]); + if (oe[1] < 0) swap(p[3],p[5]); + if (oe[2] < 0) swap(p[6],p[8]); + + if (whatd[op_id]) { + R2 fb[12]; + for( int df =0; df <12; ++df) + { + R2 f; + for( int i =0; i <6; ++i)// + { + R2 fC(C1[df][2*i] , C1[df][2*i+1]); + f += L[i]*fC; + } + f*=oe[df/3];// signe du dof * fb[df] + fb[df]= R2((Piola[0],f),(Piola[1],f)); + if(Ortho) + fb[df]= fb[df].perp(); + } + + for( int df =0; df <12; ++df) + { + val(df, 0, op_id) = fb[p[df]].x; + val(df, 1, op_id) = fb[p[df]].y; + } + } + + if (whatd[op_dx] || whatd[op_dy]) { + R2 fbx[12],fby[12];// Dx, Dy of f + for( int df =0; df <12; ++df) + { + R2 fx,fy; + for( int i =0; i <6; ++i)// + { + R2 fC(C1[df][2*i] , C1[df][2*i+1]); + // cout << df << " "<< i <<" " << R2(ci,cii)<< endl; + fx += DL[i].x * fC; + fy += DL[i].y * fC; + } + fx*=oe[df/3];// signe ;du dof * fb[df] + fy*=oe[df/3];// signe du dof * fb[df] + R2 gx= R2((Piola[0],fx),(Piola[1],fx)); + R2 gy= R2((Piola[0],fy),(Piola[1],fy)); + if(Ortho) gx=gx.perp(),gy=gy.perp(); + fbx[df]= gx; + fby[df]= gy; + // cout << df << " -- gx : " << gx << " gy :" << gy << endl; + // for(int k=9; k<12;++k) cout<< k << " " << fb[k] << endl; + } + if (whatd[op_dx]) { + for( int df =0; df <12; ++df) + { + val(df, 0, op_dx) = fbx[p[df]].x; + val(df, 1, op_dx) = fbx[p[df]].y; + } + + } + if (whatd[op_dy]) { + for( int df =0; df <12; ++df) + { + val(df, 0, op_dy) = fby[p[df]].x; + val(df, 1, op_dy) = fby[p[df]].y; + } + + } + } + if ( whatd[op_dxx] || whatd[op_dyy] || whatd[op_dxy]) { + R2 fbxx[12],fbyy[12],fbxy[12];// Dx, Dy of f + ffassert(0); // To Do + if (whatd[op_dyy]) { + for( int df =0; df <12; ++df) + { + val(df, 0, op_dyy) = fbyy[p[df]].x; + val(df, 1, op_dyy) = fbyy[p[df]].y; + } + + } + if (whatd[op_dxx]) { + for( int df =0; df <12; ++df) + { + val(df, 0, op_dxx) = fbxx[p[df]].x; + val(df, 1, op_dxx) = fbxx[p[df]].y; + } + + } + if (whatd[op_dxy]) { + for( int df =0; df <12; ++df) + { + val(df, 0, op_dxy) = fbxy[p[df]].x; + val(df, 1, op_dxy) = fbxy[p[df]].y; + } + + } + } + + + } + + + // a static variable to add the finite element to freefem++ static TypeOfFE_RTdc Elm_TypeOfFE_RT0dc_2d; // RT0dc @@ -1217,5 +1476,10 @@ void TypeOfFE_RTdc::Pi_h_alpha(const baseFElement & K,KN_ & v) const static AddNewFE Elm__TypeOfFE_RT2_2dOrtho("RT2Ortho", &Elm_TypeOfFE_RT2_2dOrtho); static AddNewFE Elm__TypeOfFE_BDM1_2d("BDM1", &Elm_TypeOfFE_BDM1_2d); static AddNewFE Elm__TypeOfFE_BDM1_2dOrtho("BDM1Ortho", &Elm_TypeOfFE_BDM1_2dOrtho); +static TypeOfFE_BDM2_2d Elm_TypeOfFE_BDM2_2d(false); // BDM1 +static TypeOfFE_BDM2_2d Elm_TypeOfFE_BDM2_2dOrtho(true); // BDM1ortho +static AddNewFE Elm__TypeOfFE_BDM2_2d("BDM2", &Elm_TypeOfFE_BDM2_2d); +static AddNewFE Elm__TypeOfFE_BDM2_2dOrtho("BDM2Ortho", &Elm_TypeOfFE_BDM2_2dOrtho); + } // namespace Fem2D // --- fin -- diff --git a/plugin/seq/Element_P3pnc_3d.cpp b/plugin/seq/Element_P3pnc_3d.cpp index f1ff40355..82d67f7fe 100644 --- a/plugin/seq/Element_P3pnc_3d.cpp +++ b/plugin/seq/Element_P3pnc_3d.cpp @@ -65,7 +65,6 @@ TypeOfFE_P3pnc_3d::TypeOfFE_P3pnc_3d( ) static R3 Pt[] = {R3(0., 0., 0.), R3(1., 0., 0.), R3(0., 1., 0.), R3(0., 0., 1.)}; // 4 ref tetrahedron vertices static const int nvfo[4][3] ={{1,2,3}, {0,2,3},{0,1,3},{ 0,1,2}}; - cout << "TypeOfFE_P3pnc_3d QFf exact:"<< QFf.exact << ", QFt exact " <pInterpolation[i] = ipt; // pk in (13.1) this->cInterpolation[i] = 0; // jk in (13.1) @@ -115,8 +112,6 @@ TypeOfFE_P3pnc_3d::TypeOfFE_P3pnc_3d( ) int kfK = nvfo[f][(kf+2)%3];// vertex dans K opose a kf int kfK1 = nvfo[f][(kf+1)%3];// vertex dans K .. int dof = 6*f+kf+3; - if(qq==0) - cout << " dof " << dof << " " << f << " "<< kfK << endl; { this->pInterpolation[i] = ipt; // pk in (13.1) this->cInterpolation[i] = 0; // jk in (13.1) @@ -126,7 +121,6 @@ TypeOfFE_P3pnc_3d::TypeOfFE_P3pnc_3d( ) } } // end loop on face - cout << " i =" << i << " " <PtInterpolation[ipt].toBary(ll); @@ -138,7 +132,6 @@ TypeOfFE_P3pnc_3d::TypeOfFE_P3pnc_3d( ) this->coefInterpolation[i] = QFt[q].a*ll[kt]; // alfak: we will fill them with 'set' (below) } } - cout << i << "== " << this->NbcoefforInterpolation << " " << QFt.n*4 + QFf.n*6*4 << endl;; ffassert(i==this->NbcoefforInterpolation); } @@ -211,7 +204,7 @@ void TypeOfFE_P3pnc_3d::set(const Mesh &Th, const Element &K, InterpolationMatri double ll[4]; // dans Khat for (int q = 0; q < QFt.n; ++q, ++ip) - { int ipt=ip; + { int ipt = nump ? nump[ip] : ip;// Correct aug 2023 FH. double ll[4]; // dans Khat M.P[ipt].toBary(ll);// M.coef[i++] = QFt[q].a*ll[0]; @@ -221,6 +214,7 @@ void TypeOfFE_P3pnc_3d::set(const Mesh &Th, const Element &K, InterpolationMatri } // ffassert(i== ncoef+ocoef); + //if(verbosity>1000) cout << " M = @@@@ \n"<< M << endl; } void TypeOfFE_P3pnc_3d::FB(const What_d whatd, const Mesh &Th, const Mesh3::Element &K, diff --git a/plugin/seq/Element_P4.cpp b/plugin/seq/Element_P4.cpp index f404dc6ec..715ed25a9 100644 --- a/plugin/seq/Element_P4.cpp +++ b/plugin/seq/Element_P4.cpp @@ -873,6 +873,7 @@ GTypeOfFE< MeshS > &Elm_P4_S(P4_S); static void init( ) { + if(verbosity && mpirank ==0 ) cout << " load : P4 "; AddNewFE("P4", &P4LagrangeP4); static ListOfTFE FE_P4("P4", &P4LagrangeP4); // to add P4 in list of Common FE diff --git a/plugin/seq/iovtk.cpp b/plugin/seq/iovtk.cpp index 122eee8c6..5630dcbb8 100644 --- a/plugin/seq/iovtk.cpp +++ b/plugin/seq/iovtk.cpp @@ -5558,7 +5558,11 @@ class VTK_WriteMeshT_Op : public E_F0mps { }; vector< Expression2 > l; - static const int n_name_param = 7; +#ifndef COMMON_HPDDM_PARALLEL_IO + static const int n_name_param = 8; +#else + static const int n_name_param = 9; +#endif static basicAC_F0::name_and_type name_param[]; Expression nargs[n_name_param]; long arg(int i, Stack stack, long a) const { @@ -5672,7 +5676,11 @@ basicAC_F0::name_and_type VTK_WriteMeshT_Op< MeshS >::name_param[] = { {"floatmesh", &typeid(bool)}, {"floatsol", &typeid(bool)}, {"bin", &typeid(bool)}, - {"swap", &typeid(bool)}}; + {"swap", &typeid(bool)}, +#ifdef COMMON_HPDDM_PARALLEL_IO + {"communicator", &typeid(pcommworld)}, +#endif + {"append", &typeid(bool)}}; template<> basicAC_F0::name_and_type VTK_WriteMeshT_Op< MeshL >::name_param[] = { @@ -5683,7 +5691,10 @@ basicAC_F0::name_and_type VTK_WriteMeshT_Op< MeshL >::name_param[] = { {"floatmesh", &typeid(bool)}, {"floatsol", &typeid(bool)}, {"bin", &typeid(bool)}, - {"swap", &typeid(bool)}}; +#ifdef COMMON_HPDDM_PARALLEL_IO + {"communicator", &typeid(pcommworld)}, +#endif + {"append", &typeid(bool)}}; template< class MMesh > void VTK_WRITE_MESHT(const string &filename, FILE *fp, const MMesh &Th, bool binary, int datasize, @@ -5993,6 +6004,15 @@ AnyType VTK_WriteMeshT_Op< MMesh >::operator( )(Stack stack) const { int datasize = floatmesh ? sizeof(float) : sizeof(double); int datasizeSol = floatsol ? sizeof(float) : sizeof(double); +#ifdef COMMON_HPDDM_PARALLEL_IO + int time, size; + std::string base_filename; + std::string CellDataArrayForPvtu = ""; + std::string PointDataArrayForPvtu = ""; + parallelIO(pffname, nargs[7] ? (MPI_Comm *)GetAny< pcommworld >((*nargs[7])(stack)) : 0, + nargs[8] && GetAny< bool >((*nargs[8])(stack)), time, size, base_filename); +#endif + int iii = 0; if (nargs[0]) { // char *data = newcopy(dataname->c_str()); @@ -6151,9 +6171,15 @@ AnyType VTK_WriteMeshT_Op< MMesh >::operator( )(Stack stack) const { for (int ii = 0; ii < nbofsol; ii++) { if (order[ii] == 1) { if (datasize == sizeof(float)) { +#ifdef COMMON_HPDDM_PARALLEL_IO + PointDataArrayForPvtu = PointDataArrayForPvtu + " \n"; +#endif VTU_DATA_ARRAY(fp, "Float32", nameofuser[ii], l[ii].nbfloat, binary); l[ii].writesolutionP1_float(fp, Th, stack, binary, swap, 1); } else if (datasize == sizeof(double)) { +#ifdef COMMON_HPDDM_PARALLEL_IO + PointDataArrayForPvtu = PointDataArrayForPvtu + " \n"; +#endif VTU_DATA_ARRAY(fp, "Float64", nameofuser[ii], l[ii].nbfloat, binary); l[ii].writesolutionP1_double(fp, Th, stack, binary, swap, 1); } @@ -6172,9 +6198,15 @@ AnyType VTK_WriteMeshT_Op< MMesh >::operator( )(Stack stack) const { for (int ii = 0; ii < nbofsol; ii++) { if (order[ii] == 0) { if (datasize == sizeof(float)) { +#ifdef COMMON_HPDDM_PARALLEL_IO + CellDataArrayForPvtu = CellDataArrayForPvtu + " \n"; +#endif VTU_DATA_ARRAY(fp, "Float32", nameofuser[ii], l[ii].nbfloat, binary); l[ii].writesolutionP0_float(fp, Th, stack, surface, binary, swap, 1); } else if (datasize == sizeof(double)) { +#ifdef COMMON_HPDDM_PARALLEL_IO + CellDataArrayForPvtu = CellDataArrayForPvtu + " \n"; +#endif VTU_DATA_ARRAY(fp, "Float64", nameofuser[ii], l[ii].nbfloat, binary); l[ii].writesolutionP0_double(fp, Th, stack, surface, binary, swap, 1); } @@ -6189,6 +6221,10 @@ AnyType VTK_WriteMeshT_Op< MMesh >::operator( )(Stack stack) const { fprintf(fp, "\n"); fprintf(fp, "\n"); fprintf(fp, "\n"); +#ifdef COMMON_HPDDM_PARALLEL_IO + PvtuWriter(pffname, size, time, base_filename, CellDataArrayForPvtu, PointDataArrayForPvtu); + PvdWriter(pffname, size, time, base_filename, CellDataArrayForPvtu, PointDataArrayForPvtu); +#endif } else { cout << " iovtk extension file is not correct (" << VTK_FILE << " != 1 or 2 ) " << endl; ExecError(" iovtk : extension file"); @@ -6872,8 +6908,8 @@ static void Load_Init( ) { if (!Global.Find("savevtk").NotNull( )) { Global.Add("savevtk", "(", new OneOperatorCode< VTK_WriteMesh_Op >); Global.Add("savevtk", "(", new OneOperatorCode< VTK_WriteMesh3_Op >); + Global.Add("savevtk", "(", new OneOperatorCode< VTK_WriteMeshT_Op< MeshS > >); } - Global.Add("savevtk", "(", new OneOperatorCode< VTK_WriteMeshT_Op< MeshS > >); Global.Add("savevtk", "(", new OneOperatorCode< VTK_WriteMeshT_Op< MeshL > >); Global.Add("vtkload", "(", new VTK_LoadMesh); Global.Add("vtkload3", "(", new VTK_LoadMesh3); diff --git a/plugin/seq/metis.cpp b/plugin/seq/metis.cpp index da7bcb6be..c4b1b5805 100644 --- a/plugin/seq/metis.cpp +++ b/plugin/seq/metis.cpp @@ -97,7 +97,7 @@ KN< R > *partmetis( KN< R > *const &part, FESPACE *const &pVh, long const &lpart } #endif - } else epart = 0; + } else epart = static_cast(0); part->resize(nv); *part = npart; return part; @@ -154,7 +154,7 @@ KN< R > *partmetis(Stack s, KN< R > *const &part, Mesh *const &pTh, long const & } #endif - } else epart = 0; + } else epart = static_cast(0); part->resize(nt); *part = epart; return part; @@ -184,7 +184,7 @@ KN< long > *partmetisd(Stack s, KN< long > *const &part, Mesh *const &pTh, long printf(" %d-way Edge-Cut: %7d, Balance: %5.2f\n", nparts, nve, ComputeElementBalance(nt, nparts, epart)); #endif - } else epart = 0; + } else epart = static_cast(0); part->resize(nt); *part = epart; return part; diff --git a/plugin/seq/mmg.cpp b/plugin/seq/mmg.cpp index eeeb33c3e..b9f193104 100644 --- a/plugin/seq/mmg.cpp +++ b/plugin/seq/mmg.cpp @@ -10,6 +10,44 @@ using namespace Fem2D; template int ffmesh_to_MMG5_pMesh(const ffmesh &, MMG5_pMesh&); +template<> +int ffmesh_to_MMG5_pMesh(const Mesh &Th, MMG5_pMesh& mesh) { + int nVertices = Th.nv; + int nTriangles = Th.nt; + int nQuadrilaterals = 0; + int nEdges = Th.neb; + + if ( MMG2D_Set_meshSize(mesh,nVertices,nTriangles, + nQuadrilaterals,nEdges) != 1 ) { + exit(EXIT_FAILURE); + } + + for (int k = 0; k < Th.nv; k++) { + if ( MMG2D_Set_vertex(mesh,Th.vertices[k].x,Th.vertices[k].y, + Th.vertices[k].lab, k+1) != 1 ) { + exit(EXIT_FAILURE); + } + } + + for (int k = 0; k < Th.nt; k++) { + const Triangle &K(Th[k]); + if ( MMG2D_Set_triangle(mesh,Th.operator()(K[0])+1,Th.operator()(K[1])+1, + Th.operator()(K[2])+1,K.lab,k+1) != 1 ) { + exit(EXIT_FAILURE); + } + } + + for (int k = 0; k < Th.neb; k++) { + const Mesh::BorderElement &K(Th.be(k)); + if ( MMG2D_Set_edge(mesh,Th.operator()(K[0])+1,Th.operator()(K[1])+1, + K.lab,k+1) != 1 ) { + exit(EXIT_FAILURE); + } + } + + return 0; +} + template<> int ffmesh_to_MMG5_pMesh(const Mesh3 &Th, MMG5_pMesh& mesh) { int nVertices = Th.nv; @@ -26,7 +64,7 @@ int ffmesh_to_MMG5_pMesh(const Mesh3 &Th, MMG5_pMesh& mesh) { for (int k = 0; k < Th.nv; k++) { if ( MMG3D_Set_vertex(mesh,Th.vertices[k].x,Th.vertices[k].y, - Th.vertices[k].z, Th.vertices[k].lab, k+1) != 1 ) { + Th.vertices[k].z, Th.vertices[k].lab, k+1) != 1 ) { exit(EXIT_FAILURE); } } @@ -65,7 +103,7 @@ int ffmesh_to_MMG5_pMesh(const MeshS &Th, MMG5_pMesh& mesh) { for (int k = 0; k < Th.nv; k++) { if ( MMGS_Set_vertex(mesh,Th.vertices[k].x,Th.vertices[k].y, - Th.vertices[k].z, Th.vertices[k].lab, k+1) != 1 ) { + Th.vertices[k].z, Th.vertices[k].lab, k+1) != 1 ) { exit(EXIT_FAILURE); } } @@ -80,7 +118,7 @@ int ffmesh_to_MMG5_pMesh(const MeshS &Th, MMG5_pMesh& mesh) { for (int k = 0; k < Th.nbe; k++) { const BoundaryEdgeS &K(Th.be(k)); - if ( MMG3D_Set_edge(mesh,Th.operator()(K[0])+1,Th.operator()(K[1])+1, + if ( MMGS_Set_edge(mesh,Th.operator()(K[0])+1,Th.operator()(K[1])+1, K.lab,k+1) != 1 ) { exit(EXIT_FAILURE); } @@ -91,6 +129,76 @@ int ffmesh_to_MMG5_pMesh(const MeshS &Th, MMG5_pMesh& mesh) { template int MMG5_pMesh_to_ffmesh(const MMG5_pMesh&, ffmesh *&); +template<> +int MMG5_pMesh_to_ffmesh(const MMG5_pMesh& mesh, Mesh *&T_TH) { + int ier; + + int nVertices = 0; + int nTriangles = 0; + int nEdges = 0; + + if ( MMG2D_Get_meshSize(mesh,&nVertices,&nTriangles,NULL, + &nEdges) !=1 ) { + ier = MMG5_STRONGFAILURE; + } + + Vertex *v = new Vertex[nVertices]; + Triangle *t = new Triangle[nTriangles]; + Triangle *tt = t; + Mesh::BorderElement *b = new Mesh::BorderElement[nEdges]; + Mesh::BorderElement *bb = b; + int k; + + int corner, required; + + for (k = 0; k < nVertices; k++) { + if ( MMG2D_Get_vertex(mesh,&(v[k].x),&(v[k].y), + &(v[k].lab),&(corner),&(required)) != 1 ) { + cout << "Unable to get mesh vertex " << k << endl; + ier = MMG5_STRONGFAILURE; + } + } + + for ( k=0; kset(v, iv[0], iv[1], iv[2], lab); + } + + for ( k=0; kset(v, iv[0], iv[1], lab); + } + + T_TH = new Mesh(nVertices, nTriangles, nEdges, v, t, b); + + if (verbosity > 1) { + cout << "transformation maillage --> mesh " << endl; + cout << "vertices =" << nVertices << endl; + cout << "triangles =" << nTriangles << endl; + cout << "edges =" << nEdges << endl; + cout << "T_TH" << T_TH->nv << " " << T_TH->nt << " " << T_TH->neb << endl; + } + + return 0; +} + template<> int MMG5_pMesh_to_ffmesh(const MMG5_pMesh& mesh, Mesh3 *&T_TH3) { int ier; @@ -101,7 +209,7 @@ int MMG5_pMesh_to_ffmesh(const MMG5_pMesh& mesh, Mesh3 *&T_TH3) { int nEdges = 0; if ( MMG3D_Get_meshSize(mesh,&nVertices,&nTetrahedra,NULL,&nTriangles,NULL, - &nEdges) !=1 ) { + &nEdges) !=1 ) { ier = MMG5_STRONGFAILURE; } @@ -170,7 +278,7 @@ int MMG5_pMesh_to_ffmesh(const MMG5_pMesh& mesh, MeshS *&T_TH3) { int nTriangles = 0; int nEdges = 0; - if ( MMGS_Get_meshSize(mesh,&nVertices,&nTriangles,&nEdges) !=1 ) { + if ( MMGS_Get_meshSize(mesh,&nVertices,&nTriangles,&nEdges) !=1 ) { ier = MMG5_STRONGFAILURE; } @@ -234,7 +342,7 @@ template class mmg_Op : public E_F0mps { public: Expression eTh; - static const int n_name_param = std::is_same::value ? 28 : 21; + static const int n_name_param = std::is_same::value ? 30 : 24; static basicAC_F0::name_and_type name_param[]; Expression nargs[n_name_param]; @@ -291,6 +399,7 @@ basicAC_F0::name_and_type mmg_Op::name_param[] = { //{"nreg" , &typeid(bool)},/*!< [0/1], Enable normal regularization */ {"renum" , &typeid(bool)},/*!< [1/0], Turn on/off point relocation with Scotch */ {"anisosize" , &typeid(bool)},/*!< [1/0], Turn on/off anisotropic metric creation when no metric is provided */ +{"nosizreq" , &typeid(bool)},/*!< [0/1], Allow/avoid overwritten of sizes at required points (advanced usage) */ {"octree" , &typeid(long)},/*!< [n], Specify the max number of points per PROctree cell (DELAUNAY) */ {"angleDetection" , &typeid(double)},/*!< [val], Value for angle detection */ {"hmin" , &typeid(double)},/*!< [val], Minimal mesh size */ @@ -298,6 +407,7 @@ basicAC_F0::name_and_type mmg_Op::name_param[] = { {"hsiz" , &typeid(double)},/*!< [val], Constant mesh size */ {"hausd" , &typeid(double)},/*!< [val], Control global Hausdorff distance (on all the boundary surfaces of the mesh) */ {"hgrad" , &typeid(double)},/*!< [val], Control gradation */ +{"hgradreq" , &typeid(double)},/*!< [val], Control required gradation */ {"ls" , &typeid(double)},/*!< [val], Value of level-set */ //{"rmc" , &typeid(double)},/*!< [-1/val], Remove small connex componants in level-set mode */ {"requiredTriangle" , &typeid(KN*)},/*!< [val], References of surfaces with required triangles */ @@ -305,8 +415,8 @@ basicAC_F0::name_and_type mmg_Op::name_param[] = { {"localParameter" , &typeid(KNM*)}/*!< [val], Local parameters on given surfaces */ }; -template<> -basicAC_F0::name_and_type mmg_Op::name_param[] = { +template +basicAC_F0::name_and_type mmg_Op::name_param[] = { {"metric" , &typeid(KN< double > *)}, {"verbose" , &typeid(long)},/*!< [-1..10], Tune level of verbosity */ {"mem" , &typeid(long)},/*!< [n/-1], Set memory size to n Mbytes or keep the default value */ @@ -314,18 +424,20 @@ basicAC_F0::name_and_type mmg_Op::name_param[] = { {"angle" , &typeid(bool)},/*!< [1/0], Turn on/off angle detection */ {"iso" , &typeid(bool)},/*!< [1/0], Level-set meshing */ {"keepRef" , &typeid(bool)},/*!< [1/0], Preserve the initial domain references in level-set mode */ -//{"optim" , &typeid(bool)},/*!< [1/0], Optimize mesh keeping its initial edge sizes */ +{"localParameter" , &typeid(KNM*)},/*!< [val], Local parameters on given surfaces */ {"noinsert" , &typeid(bool)},/*!< [1/0], Avoid/allow point insertion */ {"noswap" , &typeid(bool)},/*!< [1/0], Avoid/allow edge or face flipping */ {"nomove" , &typeid(bool)},/*!< [1/0], Avoid/allow point relocation */ {"nreg" , &typeid(bool)},/*!< [0/1], Disabled/enabled normal regularization */ {"renum" , &typeid(bool)},/*!< [1/0], Turn on/off point relocation with Scotch */ +{"nosizreq" , &typeid(bool)},/*!< [0/1], Allow/avoid overwritten of sizes at required points (advanced usage) */ {"angleDetection" , &typeid(double)},/*!< [val], Value for angle detection */ {"hmin" , &typeid(double)},/*!< [val], Minimal mesh size */ {"hmax" , &typeid(double)},/*!< [val], Maximal mesh size */ {"hsiz" , &typeid(double)},/*!< [val], Constant mesh size */ {"hausd" , &typeid(double)},/*!< [val], Control global Hausdorff distance (on all the boundary surfaces of the mesh) */ {"hgrad" , &typeid(double)},/*!< [val], Control gradation */ +{"hgradreq" , &typeid(double)},/*!< [val], Control required gradation */ {"ls" , &typeid(double)},/*!< [val], Value of level-set */ {"requiredEdge" , &typeid(KN*)},/*!< [val], References of boundaries with required edges */ {"requiredVertex" , &typeid(KN*)}/*!< [val], Indices of required vertices */ @@ -340,6 +452,163 @@ class mmg_ff : public OneOperator { E_F0 *code(const basicAC_F0 &args) const { return new mmg_Op(args, t[0]->CastTo(args[0])); } }; +template<> +AnyType mmg_Op::operator( )(Stack stack) const { + // initialisation + MeshPoint *mp(MeshPointStack(stack)), mps = *mp; + Mesh *pTh = GetAny< Mesh * >((*eTh)(stack)); + + ffassert(pTh); + Mesh &Th = *pTh; + int nv = Th.nv; + int nt = Th.nt; + int nbe = Th.neb; + + KN< double > *pmetric = 0; + KN< long > *prequiredEdge = 0; + KN< long > *prequiredVertex = 0; + KNM< double > *plocalParameter = 0; + + if (nargs[0]) { + pmetric = GetAny< KN< double > * >((*nargs[0])(stack)); + } + if (nargs[22]) { + prequiredEdge = GetAny< KN< long > * >((*nargs[22])(stack)); + } + if (nargs[23]) { + prequiredVertex = GetAny< KN< long > * >((*nargs[23])(stack)); + } + if (nargs[7]) { + plocalParameter = GetAny< KNM< double > * >((*nargs[7])(stack)); + } + + MMG5_pMesh mesh; + MMG5_pSol sol; + + mesh = nullptr; + sol = nullptr; + + MMG2D_Init_mesh(MMG5_ARG_start, + MMG5_ARG_ppMesh,&mesh,MMG5_ARG_ppMet,&sol, + MMG5_ARG_end); + + ffmesh_to_MMG5_pMesh(Th, mesh); + + if (pmetric && pmetric->N( ) > 0) { + const KN< double > &metric = *pmetric; + if (metric.N( ) == Th.nv) { + if ( MMG2D_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Scalar) != 1 ) { + printf("Unable to allocate the metric array.\n"); + exit(EXIT_FAILURE); + } + if ( MMG2D_Set_scalarSols(sol,metric) != 1 ) { + printf("Unable to set metric.\n"); + exit(EXIT_FAILURE); + } + } + else { + if ( MMG2D_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Tensor) != 1 ) { + printf("Unable to allocate the metric array.\n"); + exit(EXIT_FAILURE); + } + static const int perm[3] = {0, 1, 2}; + for (int k=0; kN( ) > 0) { + const KN< long > &requiredEdge = *prequiredEdge; + std::sort(requiredEdge + 0, requiredEdge + requiredEdge.N()); + int na; + if ( MMG2D_Get_meshSize(mesh,NULL,NULL,NULL,&na) !=1 ) { + exit(EXIT_FAILURE); + } + for (int k=1; k<=na; k++) { + int ref, dummy; + if ( MMG2D_Get_edge(mesh, &dummy, &dummy, &ref, + &dummy, &dummy) != 1 ) { + exit(EXIT_FAILURE); + } + if (std::binary_search(requiredEdge + 0, requiredEdge + requiredEdge.N(), ref)) { + if ( MMG2D_Set_requiredEdge(mesh,k) != 1 ) { + exit(EXIT_FAILURE); + } + } + } + } + if (prequiredVertex && prequiredVertex->N( ) > 0) { + const KN< long > &requiredVertex = *prequiredVertex; + for (int k=0; kM( ) > 0) { + const KNM< double > &localParameter = *plocalParameter; + ffassert(localParameter.N() == 4); + if ( MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_numberOfLocalParam,localParameter.M()) != 1 ) { + exit(EXIT_FAILURE); + } + for(int j = 0; j < localParameter.M(); ++j) { + if ( MMG2D_Set_localParameter(mesh,sol,MMG5_Edg,localParameter(0,j),localParameter(1,j),localParameter(2,j),localParameter(3,j)) != 1 ) { + exit(EXIT_FAILURE); + } + } + } + + long iso=0L; + + int i=1; + if (nargs[i]) MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_verbose, arg(i,stack,0L)); i++; /*!< [-1..10], Tune level of verbosity */ + if (nargs[i]) MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_mem, arg(i,stack,0L)); i++; /*!< [n/-1], Set memory size to n Mbytes or keep the default value */ + if (nargs[i]) MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_debug, arg(i,stack,false)); i++; /*!< [1/0], Turn on/off debug mode */ + if (nargs[i]) MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_angle, arg(i,stack,false)); i++; /*!< [1/0], Turn on/off angle detection */ + if (nargs[i]) {iso = arg(i,stack,false); MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_iso,iso);} i++; /*!< [1/0], Level-set meshing */ + i++; /*!< [1/0], Preserve the initial domain references in level-set mode */ + i++; /*!< [val], Local parameters on given surfaces */ + if (nargs[i]) MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_noinsert, arg(i,stack,false)); i++; /*!< [1/0], Avoid/allow point insertion */ + if (nargs[i]) MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_noswap, arg(i,stack,false)); i++; /*!< [1/0], Avoid/allow edge or face flipping */ + if (nargs[i]) MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_nomove, arg(i,stack,false)); i++; /*!< [1/0], Avoid/allow point relocation */ + if (nargs[i]) MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_nreg, arg(i,stack,false)); i++; /*!< [0/1], Disabled/enabled normal regularization */ + i++; /*!< [1/0], Turn on/off point relocation with Scotch */ + if (nargs[i]) MMG2D_Set_iparameter(mesh,sol,MMG2D_IPARAM_nosizreq, arg(i,stack,false)); i++; /*!< [0/1], Allow/avoid overwritten of sizes at required points (advanced usage) */ + if (nargs[i]) MMG2D_Set_dparameter(mesh,sol,MMG2D_DPARAM_angleDetection,arg(i,stack,0.)); i++; /*!< [val], Value for angle detection */ + if (nargs[i]) MMG2D_Set_dparameter(mesh,sol,MMG2D_DPARAM_hmin, arg(i,stack,0.)); i++; /*!< [val], Minimal mesh size */ + if (nargs[i]) MMG2D_Set_dparameter(mesh,sol,MMG2D_DPARAM_hmax, arg(i,stack,0.)); i++; /*!< [val], Maximal mesh size */ + if (nargs[i]) MMG2D_Set_dparameter(mesh,sol,MMG2D_DPARAM_hsiz, arg(i,stack,0.)); i++; /*!< [val], Constant mesh size */ + if (nargs[i]) MMG2D_Set_dparameter(mesh,sol,MMG2D_DPARAM_hausd, arg(i,stack,0.)); i++; /*!< [val], Control global Hausdorff distance (on all the boundary surfaces of the mesh) */ + if (nargs[i]) MMG2D_Set_dparameter(mesh,sol,MMG2D_DPARAM_hgrad, arg(i,stack,0.)); i++; /*!< [val], Control gradation */ + if (nargs[i]) MMG2D_Set_dparameter(mesh,sol,MMG2D_DPARAM_hgradreq, arg(i,stack,0.)); i++; /*!< [val], Control required gradation */ + if (nargs[i]) MMG2D_Set_dparameter(mesh,sol,MMG2D_DPARAM_ls, arg(i,stack,0.)); i++; /*!< [val], Value of level-set */ + + int ier; + if (!iso) + ier = MMG2D_mmg2dlib(mesh,sol); + else + ier = MMG2D_mmg2dls(mesh,sol,NULL); + + Mesh *Th_T = nullptr; + + MMG5_pMesh_to_ffmesh(mesh,Th_T); + + MMG2D_Free_all(MMG5_ARG_start, + MMG5_ARG_ppMesh,&mesh,MMG5_ARG_ppMet,&sol, + MMG5_ARG_end); + + R2 Pn, Px; + Th_T->BoundingBox(Pn, Px); + Th_T->quadtree = new Fem2D::FQuadTree(Th_T, Pn, Px, Th_T->nv); + + Add2StackOfPtr2FreeRC(stack, Th_T); + return Th_T; +} + template<> AnyType mmg_Op::operator( )(Stack stack) const { // initialisation @@ -360,14 +629,14 @@ AnyType mmg_Op::operator( )(Stack stack) const { if (nargs[0]) { pmetric = GetAny< KN< double > * >((*nargs[0])(stack)); } - if (nargs[25]) { - prequiredTriangle = GetAny< KN< long > * >((*nargs[25])(stack)); + if (nargs[27]) { + prequiredTriangle = GetAny< KN< long > * >((*nargs[27])(stack)); } - if (nargs[26]) { - prequiredVertex = GetAny< KN< long > * >((*nargs[26])(stack)); + if (nargs[28]) { + prequiredVertex = GetAny< KN< long > * >((*nargs[28])(stack)); } - if (nargs[27]) { - plocalParameter = GetAny< KNM< double > * >((*nargs[27])(stack)); + if (nargs[29]) { + plocalParameter = GetAny< KNM< double > * >((*nargs[29])(stack)); } MMG5_pMesh mesh; @@ -383,28 +652,28 @@ AnyType mmg_Op::operator( )(Stack stack) const { MMG5_ARG_end); ffmesh_to_MMG5_pMesh(Th, mesh); - + if (pmetric && pmetric->N( ) > 0) { const KN< double > &metric = *pmetric; if (metric.N( ) == Th.nv) { - if ( MMG3D_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Scalar) != 1 ) { + if ( MMG3D_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Scalar) != 1 ) { printf("Unable to allocate the metric array.\n"); exit(EXIT_FAILURE); } - if ( MMG3D_Set_scalarSols(sol,metric) != 1 ) { + if ( MMG3D_Set_scalarSols(sol,metric) != 1 ) { printf("Unable to set metric.\n"); exit(EXIT_FAILURE); } } else { - if ( MMG3D_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Tensor) != 1 ) { + if ( MMG3D_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Tensor) != 1 ) { printf("Unable to allocate the metric array.\n"); exit(EXIT_FAILURE); } static const int perm[6] = {0, 1, 3, 2, 4, 5}; for (int k=0; k::operator( )(Stack stack) const { //if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_nreg, arg(i,stack,false)); i++; /*!< [0/1], Enable normal regularization */ if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_renum, arg(i,stack,false)); i++; /*!< [1/0], Turn on/off point relocation with Scotch */ if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_anisosize, arg(i,stack,false)); i++; /*!< [1/0], Turn on/off anisotropic metric creation when no metric is provided */ + if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_nosizreq, arg(i,stack,false)); i++; /*!< [0/1], Allow/avoid overwritten of sizes at required points (advanced usage) */ if (nargs[i]) MMG3D_Set_iparameter(mesh,sol,MMG3D_IPARAM_octree, arg(i,stack,0L)); i++; /*!< [n], Specify the max number of points per PROctree cell (DELAUNAY) */ if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_angleDetection,arg(i,stack,0.)); i++; /*!< [val], Value for angle detection */ if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_hmin, arg(i,stack,0.)); i++; /*!< [val], Minimal mesh size */ @@ -477,6 +747,7 @@ AnyType mmg_Op::operator( )(Stack stack) const { if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_hsiz, arg(i,stack,0.)); i++; /*!< [val], Constant mesh size */ if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_hausd, arg(i,stack,0.)); i++; /*!< [val], Control global Hausdorff distance (on all the boundary surfaces of the mesh) */ if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_hgrad, arg(i,stack,0.)); i++; /*!< [val], Control gradation */ + if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_hgradreq, arg(i,stack,0.)); i++; /*!< [val], Control required gradation */ if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_ls, arg(i,stack,0.)); i++; /*!< [val], Value of level-set */ //if (nargs[i]) MMG3D_Set_dparameter(mesh,sol,MMG3D_DPARAM_rmc, arg(i,stack,0.)); i++; /*!< [-1/val], Remove small connex componants in level-set mode */ @@ -487,9 +758,9 @@ AnyType mmg_Op::operator( )(Stack stack) const { ier = MMG3D_mmg3dlib(mesh,sol); else ier = MMG3D_mmg3dls(mesh,sol,met); - + Mesh3 *Th_T = nullptr; - + MMG5_pMesh_to_ffmesh(mesh,Th_T); MMG3D_Free_all(MMG5_ARG_start, @@ -497,7 +768,7 @@ AnyType mmg_Op::operator( )(Stack stack) const { MMG5_ARG_end); Th_T->BuildGTree(); - + Add2StackOfPtr2FreeRC(stack, Th_T); return Th_T; } @@ -521,11 +792,11 @@ AnyType mmg_Op::operator( )(Stack stack) const { if (nargs[0]) { pmetric = GetAny< KN< double > * >((*nargs[0])(stack)); } - if (nargs[19]) { - prequiredEdge = GetAny< KN< long > * >((*nargs[19])(stack)); + if (nargs[21]) { + prequiredEdge = GetAny< KN< long > * >((*nargs[21])(stack)); } - if (nargs[20]) { - prequiredVertex = GetAny< KN< long > * >((*nargs[20])(stack)); + if (nargs[22]) { + prequiredVertex = GetAny< KN< long > * >((*nargs[22])(stack)); } MMG5_pMesh mesh; @@ -539,28 +810,28 @@ AnyType mmg_Op::operator( )(Stack stack) const { MMG5_ARG_end); ffmesh_to_MMG5_pMesh(Th, mesh); - + if (pmetric && pmetric->N( ) > 0) { const KN< double > &metric = *pmetric; if (metric.N( ) == Th.nv) { - if ( MMGS_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Scalar) != 1 ) { + if ( MMGS_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Scalar) != 1 ) { printf("Unable to allocate the metric array.\n"); exit(EXIT_FAILURE); } - if ( MMGS_Set_scalarSols(sol,metric) != 1 ) { + if ( MMGS_Set_scalarSols(sol,metric) != 1 ) { printf("Unable to set metric.\n"); exit(EXIT_FAILURE); } } else { - if ( MMGS_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Tensor) != 1 ) { + if ( MMGS_Set_solSize(mesh,sol,MMG5_Vertex,Th.nv,MMG5_Tensor) != 1 ) { printf("Unable to allocate the metric array.\n"); exit(EXIT_FAILURE); } static const int perm[6] = {0, 1, 3, 2, 4, 5}; for (int k=0; k::operator( )(Stack stack) const { exit(EXIT_FAILURE); } if (std::binary_search(requiredEdge + 0, requiredEdge + requiredEdge.N(), ref)) { - if ( MMG3D_Set_requiredEdge(mesh,k) != 1 ) { + if ( MMGS_Set_requiredEdge(mesh,k) != 1 ) { exit(EXIT_FAILURE); } } @@ -605,18 +876,20 @@ AnyType mmg_Op::operator( )(Stack stack) const { if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_angle, arg(i,stack,false)); i++; /*!< [1/0], Turn on/off angle detection */ if (nargs[i]) {iso = arg(i,stack,false); MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_iso,iso);} i++; /*!< [1/0], Level-set meshing */ if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_keepRef, arg(i,stack,false)); i++; /*!< [1/0], Preserve the initial domain references in level-set mode */ - //if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_optim, arg(i,stack,false)); i++; /*!< [1/0], Optimize mesh keeping its initial edge sizes */ + i++; /*!< [val], Local parameters on given surfaces */ if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_noinsert, arg(i,stack,false)); i++; /*!< [1/0], Avoid/allow point insertion */ if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_noswap, arg(i,stack,false)); i++; /*!< [1/0], Avoid/allow edge or face flipping */ if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_nomove, arg(i,stack,false)); i++; /*!< [1/0], Avoid/allow point relocation */ if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_nreg, arg(i,stack,false)); i++; /*!< [0/1], Disabled/enabled normal regularization */ if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_renum, arg(i,stack,false)); i++; /*!< [1/0], Turn on/off point relocation with Scotch */ + if (nargs[i]) MMGS_Set_iparameter(mesh,sol,MMGS_IPARAM_nosizreq, arg(i,stack,false)); i++; /*!< [0/1], Allow/avoid overwritten of sizes at required points (advanced usage) */ if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_angleDetection,arg(i,stack,0.)); i++; /*!< [val], Value for angle detection */ if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hmin, arg(i,stack,0.)); i++; /*!< [val], Minimal mesh size */ if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hmax, arg(i,stack,0.)); i++; /*!< [val], Maximal mesh size */ if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hsiz, arg(i,stack,0.)); i++; /*!< [val], Constant mesh size */ if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hausd, arg(i,stack,0.)); i++; /*!< [val], Control global Hausdorff distance (on all the boundary surfaces of the mesh) */ if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hgrad, arg(i,stack,0.)); i++; /*!< [val], Control gradation */ + if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_hgradreq, arg(i,stack,0.)); i++; /*!< [val], Control required gradation */ if (nargs[i]) MMGS_Set_dparameter(mesh,sol,MMGS_DPARAM_ls, arg(i,stack,0.)); i++; /*!< [val], Value of level-set */ int ier; @@ -626,7 +899,7 @@ AnyType mmg_Op::operator( )(Stack stack) const { ier = MMGS_mmgsls(mesh,sol,NULL); MeshS *Th_T = nullptr; - + MMG5_pMesh_to_ffmesh(mesh,Th_T); MMGS_Free_all(MMG5_ARG_start, @@ -634,7 +907,7 @@ AnyType mmg_Op::operator( )(Stack stack) const { MMG5_ARG_end); Th_T->BuildGTree(); - + Add2StackOfPtr2FreeRC(stack, Th_T); return Th_T; } @@ -645,6 +918,7 @@ static void Load_Init( ) { } Global.Add("mmg3d", "(", new mmg_ff); + Global.Add("mmg2d", "(", new mmg_ff); Global.Add("mmgs", "(", new mmg_ff); } diff --git a/plugin/seq/msh3.cpp b/plugin/seq/msh3.cpp index 78206391b..f6e414f95 100644 --- a/plugin/seq/msh3.cpp +++ b/plugin/seq/msh3.cpp @@ -9442,7 +9442,7 @@ AnyType Line_Op::operator( )(Stack stack) const { int iv[2]; iv[0] = i, iv[1] = i + 1; int lab = 0; - (tt++)->set(v, iv, lab); + (tt++)->set(v, iv, reg);// correct 6 oct 2023 FH.. } B *bb = b; int ibeg[1], iend[1]; diff --git a/plugin/seq/plotPDF.cpp b/plugin/seq/plotPDF.cpp index 49f95334e..f47baa1e9 100644 --- a/plugin/seq/plotPDF.cpp +++ b/plugin/seq/plotPDF.cpp @@ -3457,9 +3457,9 @@ void plot_vector2flow( std::stringstream &Content, const Fem2D::Mesh &Th, (fmax - fmin)/nbfill; if( varrow ) - std::cout << "plotPDF(): Option 'varrow' is not implmeneted yet." << endl; + std::cout << "plotPDF(): Option 'varrow' is not implemented yet." << endl; if( nbarrow != 0 ) - std::cout << "plotPDF(): Option 'nbarrow' is not implmeneted yet." << endl; + std::cout << "plotPDF(): Option 'nbarrow' is not implemented yet." << endl; std::stringstream &st = Content; st.str(""); diff --git a/readme/README_MAC.md b/readme/README_MAC.md index 846007fa7..401404b7e 100644 --- a/readme/README_MAC.md +++ b/readme/README_MAC.md @@ -19,6 +19,18 @@ # How to compile FreeFem++ on MacOSX +## Sep 2023 + +the mac version is compile wit the script + bin/compile-version-MacOS gitbranch path-where-to-compile configure-parameter + + configure-parameter compilation: + generic for --enable-genenic + debug for --enable-debug + optim for --enable-optim + + and exact : + bin/compile-version-MacOS master $HOME/tmp generic ## Feb 2022 The Version 4.10 is compiled Under Monteray (Apple M1: arm64e) and on (Apple : x86 ) @@ -26,11 +38,11 @@ The Version 4.10 is compiled Under Monteray (Apple M1: arm64e) and on (Apple : x ### for Apple M1 (arm64e) Add with home brew -1) install brew see +1. install brew see [install Homebrew]{https://brew.sh/index_fr} -2) install this brew package to compile FreeFEM++ +2. install this brew package to compile FreeFEM++ ```bash brew install \ @@ -52,7 +64,7 @@ gd isl libunistring openssl@1.1 zip \ gdbm jansson libvmaf p11-kit zstd ``` -3) download the sources of FreeFEM++ with git +3. download the sources of FreeFEM++ with git ```bash git clone git@github.com:FreeFem/FreeFem-sources.git ff++ @@ -62,7 +74,7 @@ or the develop version ```bash git clone -b develop git@github.com:FreeFem/FreeFem-sources.git ff++ ``` -4) Compile FreeFem++ without petsc and MPI +4. Compile FreeFem++ without petsc and MPI ```bash cd ff++ @@ -72,7 +84,7 @@ or the develop version sudo make install ``` -4) Compile FreeFem++ with PETSc and MPI +4. Compile FreeFem++ with PETSc and MPI ```bash cd ff++ diff --git a/src/Graphics/Pcrgraph.cpp b/src/Graphics/Pcrgraph.cpp index db4d0a759..b70f883bf 100644 --- a/src/Graphics/Pcrgraph.cpp +++ b/src/Graphics/Pcrgraph.cpp @@ -909,7 +909,7 @@ void openSVG(const char *filename ) } else { - cerr << " Err openning SVG file " << fsvg << endl; + cerr << " Err opening SVG file " << fsvg << endl; } return; } diff --git a/src/Graphics/Xrgraph.cpp b/src/Graphics/Xrgraph.cpp index 2ab95c2c8..2c93e6c7e 100644 --- a/src/Graphics/Xrgraph.cpp +++ b/src/Graphics/Xrgraph.cpp @@ -1291,7 +1291,7 @@ void openSVG(const char *filename ) } else { - cerr << " Err openning SVG file " << fsvg << endl; + cerr << " Err opening SVG file " << fsvg << endl; } return; } diff --git a/src/Graphics/glrgraph.hpp b/src/Graphics/glrgraph.hpp index dd5dfce68..4444e08ae 100644 --- a/src/Graphics/glrgraph.hpp +++ b/src/Graphics/glrgraph.hpp @@ -1718,7 +1718,7 @@ void openSVG(const char *filename ) } else { - cerr << " Err openning SVG file " << fsvg << endl; + cerr << " Err opening SVG file " << fsvg << endl; } return; } diff --git a/src/Graphics/macrgraf.cpp b/src/Graphics/macrgraf.cpp index a420cd451..f0e004aa1 100644 --- a/src/Graphics/macrgraf.cpp +++ b/src/Graphics/macrgraf.cpp @@ -1268,7 +1268,7 @@ void openSVG(const char *filename ) } else { - cerr << " Err openning SVG file " << fsvg << endl; + cerr << " Err opening SVG file " << fsvg << endl; } return; } diff --git a/src/Graphics/sansrgraph.cpp b/src/Graphics/sansrgraph.cpp index d7f1a8152..852ad380b 100644 --- a/src/Graphics/sansrgraph.cpp +++ b/src/Graphics/sansrgraph.cpp @@ -751,7 +751,7 @@ void openSVG(const char *filename ) } else { - cerr << " Err openning SVG file " << fsvg << endl; + cerr << " Err opening SVG file " << fsvg << endl; } return; } diff --git a/src/femlib/CG.cpp b/src/femlib/CG.cpp index 162148ea3..4a0024960 100644 --- a/src/femlib/CG.cpp +++ b/src/femlib/CG.cpp @@ -228,7 +228,7 @@ int ConjugueGradient(CGMatVirt &A, // fonction et pointeur assert( !(gCg != gCg) ); // verif if NaN => bad Matrix.. if(gCg < 1e-30) { if(niveauimpression) - std::cout << " GC: on a converge on 0 iteration ||g||_C^2" << gCg << std::endl; + std::cout << " GC: converge after 0 iteration ||g||_C^2" << gCg << std::endl; nret = 2; nbitermax=0; } @@ -247,7 +247,7 @@ int ConjugueGradient(CGMatVirt &A, // fonction et pointeur if(gCg < eps2) // We have converged ... { if(niveauimpression) - std::cout << " GC: converge in " < +class OneBinaryOperatorCode : public OneOperator{ + typedef typename T::result_type R; + typedef typename T::first_argument_type A; + typedef typename T::second_argument_type B; +/* + class Op : public E_F0 { + typedef typename T::result_type Result; + Expression a,b; + public: + AnyType operator()(Stack s) const + {return SetAny(static_cast(::f( GetAny((*a)(s)) , + GetAny((*b)(s)) + );} + Op(Expression aa,Expression bb,Expression cc) : a(aa),b(bb),c(cc) {} + bool MeshIndependent() const { + return a->MeshIndependent() && b->MeshIndependent() && c->MeshIndependent();} + }; +*/ + public: + E_F0 * code(const basicAC_F0 & args) const + { return new CODE(t[0]->CastTo(args[0]),t[1]->CastTo(args[1]));} + OneBinaryOperatorCode(): + OneOperator(map_type[typeid(R).name()], + map_type[typeid(A).name()], + map_type[typeid(B).name()] + ) {} +}; template class OneTernaryOperator : public OneOperator{ typedef typename T::result_type R; typedef typename T::first_argument_type A; typedef typename T::second_argument_type B; typedef typename T::third_argument_type C; - +/* class Op : public E_F0 { - typedef typename C::result_type Result; + typedef typename T::result_type Result; Expression a,b,c; public: AnyType operator()(Stack s) const @@ -2359,7 +2387,7 @@ class OneTernaryOperator : public OneOperator{ bool MeshIndependent() const { return a->MeshIndependent() && b->MeshIndependent() && c->MeshIndependent();} }; - +*/ public: E_F0 * code(const basicAC_F0 & args) const { return new CODE(t[0]->CastTo(args[0]),t[1]->CastTo(args[1]),t[2]->CastTo(args[2]));} @@ -2377,9 +2405,9 @@ class OneQuadOperator : public OneOperator{ typedef typename T::second_argument_type B; typedef typename T::third_argument_type C; typedef typename T::fourth_argument_type D; - +/* class Op : public E_F0 { - typedef typename C::result_type Result; + typedef typename T::result_type Result; Expression a,b,c,d; public: AnyType operator()(Stack s) const @@ -2392,7 +2420,7 @@ class OneQuadOperator : public OneOperator{ bool MeshIndependent() const { return a->MeshIndependent() && b->MeshIndependent() && c->MeshIndependent() && d->MeshIndependent();} }; - + */ public: E_F0 * code(const basicAC_F0 & args) const { return new CODE(t[0]->CastTo(args[0]),t[1]->CastTo(args[1]),t[2]->CastTo(args[2]),t[3]->CastTo(args[3]));} @@ -2616,11 +2644,11 @@ class OneBinaryOperator : public OneOperator{ CompileError( " They are used Named parameter "); return new Op(t0->CastTo(args[0]),t1->CastTo(args[1]));} - OneBinaryOperator(): + OneBinaryOperator(int ppref=0): OneOperator(map_type[typeid(R).name()],map_type[typeid(A).name()],map_type[typeid(B).name()]), t0(t[0]), t1(t[1]) - {pref = SameType::OK ;} + {pref = ppref+SameType::OK ;} OneBinaryOperator(aType tt0,aType tt1): OneOperator(map_type[typeid(R).name()], diff --git a/src/fflib/Operator.hpp b/src/fflib/Operator.hpp index 63fcc4f7d..f08d9b385 100644 --- a/src/fflib/Operator.hpp +++ b/src/fflib/Operator.hpp @@ -834,6 +834,14 @@ struct PrintPnd { static ostream* f(ostream* const & a,const A & b) { if(verbosity>9999) cout << "PrintPnd: " << b << endl; *a << *b; return a;} }; +template +struct PrintPn { + using first_argument_type = ostream*; + using second_argument_type = A; + using result_type = ostream*; + static ostream* f(ostream* const & a,const A & b) + { if(verbosity>9999) cout << "PrintPn: " << b << endl; *a << b; return a;} +}; diff --git a/src/fflib/array_long.cpp b/src/fflib/array_long.cpp index 3ce862a58..c4269db77 100644 --- a/src/fflib/array_long.cpp +++ b/src/fflib/array_long.cpp @@ -164,13 +164,23 @@ void initArrayOperatorlong() // ArrayDCL(); Dcl_TypeandPtr_, KN *>(0, 0, 0, ::Destroy >, ::ClearReturnKK_, KN_ >, ::ClearReturnpKK >); + map_type[typeid(KN_).name()]->AddCast( + new E_F1_funcT,KN*>(UnRefpKN,KN_ > ) ); + atype* >()->Add("[", "", new OneOperator2_*, long >(get_elements)); TheOperators->Add("<-", new OneOperator2_ *, KN *, long>(&set_init_init), new InitArrayfromArray*, true>); + TheOperators->Add("<-", new OneBinaryOperator ,KN_ > > );//, + // new OneOperator2_ *,KN *,KN *>(&set_initp)); // Add FH nov 2023 + TheOperators->Add("=", new OneBinaryOperator ,KN_ > > ); //, + // new OneBinaryOperator ,KN* > > ); // Add FH nov 2023 + map_type_of_map[make_pair(atype(), atype())] = atype*>(); // vector [string] Add *>("n", ".", new OneOperator1 *>(get_n)); + Add >("n", ".", new OneOperator1 >(get__n)); extern KN *pkarg; Global.New("ARGV", CPValue >(*pkarg));// add FH mars 2010 Global.Add("toZarray", "(", new OneOperator_2KN_); + Global.Add("toSarray", "(", new OneOperator_2KN_); TheOperators->Add("=", new OneBinaryOperator, Inv_KN_long> >); TheOperators->Add("<-", new OneBinaryOperator*, Inv_KN_long> >); @@ -196,7 +206,15 @@ void initArrayOperatorlong() TheOperators->Add("<-", new InitMapfromArray*, string *, string*, true> ); - TheOperators->Add("<<", new OneBinaryOperator * > >); // add may 2018 FH + TheOperators->Add("<<", new OneBinaryOperator * > >(50)); // add may 2018 FH, coor + TheOperators->Add("<<", new OneBinaryOperator > >); // add now 2023 FH and change priorite to remove ambiguity + + // sep 2023 FH for G. Sadaka . + Dcl_Type< Resize > > (); + Add *>("resize",".",new OneOperator1< Resize >,KN *>(to_Resize)); + Add > >("(","",new OneOperator2_ *,Resize > , long >(resize1)); + + } // void xxxx() { diff --git a/src/fflib/array_tlp.hpp b/src/fflib/array_tlp.hpp index ce2102e03..2e3ef80a6 100644 --- a/src/fflib/array_tlp.hpp +++ b/src/fflib/array_tlp.hpp @@ -705,8 +705,13 @@ class SetArrayofKNfromKN : public OneOperator { template long get_MyMap_n(MyMap *p) {return p ? (p->m ? p->m->size():0):0 ;} template long get_n(KN * p){ return p->N();}// +template long get__n(KN_ p){ return p.N();}// + template long get_n(KNM * p){ return p->N();}// template long get_m(KNM * p){ return p->M();}// +template long get__n(KNM_ p){ return p.N();}// +template long get__m(KNM_ p){ return p.M();}// + template K get_max(KN * p){ return p->max();}// template K get_min(KN * p){ return p->min();}// @@ -1172,7 +1177,7 @@ template KNM * set_mat(KNM * a,KNM * b ){ *a=*b; return a;} -template +template class OneOperator_2KN_ : public OneOperator {public: class Op : public E_F0 { public: @@ -1182,13 +1187,13 @@ class OneOperator_2KN_ : public OneOperator {public: Op( const E_Array &bb) : N(bb.size()), tab(new Expression[N]) { for(int i=0;i()->CastTo( bb[i]); + tab[i]=atype()->CastTo( bb[i]); } AnyType operator()(Stack s) const { K * p = Add2StackOfPtr2FreeA(s,new K[N]); // mark to be delete .. KN_ A(p,N); // FH: Correct jan 2019 bug is: stupide A what delete if KN type for(int i=0;i( (*tab[i])(s)); + A[i]= GetAny( (*tab[i])(s)); return SetAny >(A);} }; E_F0 * code(const basicAC_F0 & a) const @@ -1544,8 +1549,11 @@ void ArrayOperator() Add *>("<-","(",new InitArrayfromArray*,true>); Add *>("<-","(",new InitMatfromAArray); Add *>("n",".",new OneOperator1 *>(get_n)); + Add >("n",".",new OneOperator1 >(get__n)); Add *>("n",".",new OneOperator1 *>(get_n)); Add *>("m",".",new OneOperator1 *>(get_m)); + Add >("n",".",new OneOperator1 >(get__n)); + Add >("m",".",new OneOperator1 >(get__m)); //ajout ars 2012 FH Add > *>("n",".",new OneOperator1 > *>(get_n)); Add > *>("n",".",new OneOperator1 > *>(get_n)); @@ -1558,7 +1566,6 @@ void ArrayOperator() Add > > >("(","",new OneOperator2_ > *,Resize > > , long >(resize1)); Add > > >("(","",new OneOperator2_ > *,Resize > > , long >(resize1)); Dcl_Type >(); - // AddOpeqarray("="); TheOperators->Add("=", new InitArrayfromArray*,false>(10)); diff --git a/src/fflib/compositeFESpace.cpp b/src/fflib/compositeFESpace.cpp index 8266c6fec..574f46990 100644 --- a/src/fflib/compositeFESpace.cpp +++ b/src/fflib/compositeFESpace.cpp @@ -1270,10 +1270,6 @@ void separateFEMpartBemPart(const list & largs, list &largs_FEM, lis // largs_BEM.push_back(*ii); } } - if(nbBEM >1){ - cerr << "Two BEM operator in a sub-matrix defined by a time product of the same two FESpace is not allowed" << endl; - ffassert(0); - } #endif #endif @@ -1671,33 +1667,39 @@ template< class R> HashMatrix * varfBemToBlockDenseMatrix(const int& iUh, const int &jVh, const int &typeUh, const int &typeVh, const generic_v_fes * LLUh, const generic_v_fes * LLVh, const list & b_largs_zz, Stack stack, Data_Bem_Solver &dsbem){ - // creation of the H-matrix - HMatrixVirt ** Hmat = new HMatrixVirt *(); - - // construction of the H-matrix - - // IF LLVh is FESpace(Mesh2) ==> LLVh == LLUh (by hypothesis) - // Remove this comment in the future. + HashMatrix * hm_A_block = nullptr; + for (list::const_iterator i=b_largs_zz.begin(); i!=b_largs_zz.end(); i++) { + list bi; + bi.clear(); + bi.push_back(*i); - //cout << "iUh=" << iUh << " , " << "jVh=" << jVh << endl; - //cout << LLUh << " " << LLVh << endl; + // creation of the H-matrix + HMatrixVirt ** Hmat = new HMatrixVirt *(); - varfBemToCompositeBlockLinearSystem_hmat( iUh, jVh, typeUh, typeVh, - LLUh, LLVh, b_largs_zz, stack, dsbem,Hmat); + // construction of the H-matrix + varfBemToCompositeBlockLinearSystem_hmat( iUh, jVh, typeUh, typeVh, + LLUh, LLVh, bi, stack, dsbem,Hmat); - // creation of dense matrix - KNM* M= HMatrixVirtToDense< KNM, R >(Hmat); - HashMatrix * hm_A_block = new HashMatrix(*M); + // creation of dense matrix + KNM* M= HMatrixVirtToDense< KNM, R >(Hmat); - // delete dense matrix - M->destroy(); - delete M; + HashMatrix * hm_i = new HashMatrix(*M); + if (!hm_A_block) + hm_A_block = hm_i; + else { + hm_A_block->Add(hm_i, R(1.0), false, 0, 0); + delete hm_i; + } - // delete Hmat - if( *Hmat) - delete *Hmat; - delete Hmat; + // delete dense matrix + M->destroy(); + delete M; + // delete Hmat + if( *Hmat) + delete *Hmat; + delete Hmat; + } return hm_A_block; } diff --git a/src/fflib/lgfem.cpp b/src/fflib/lgfem.cpp index 69b28ab2f..839b961d5 100644 --- a/src/fflib/lgfem.cpp +++ b/src/fflib/lgfem.cpp @@ -4828,6 +4828,8 @@ AnyType Plot::operator( )(Stack s) const { cadre(xx0, xx1, yy0, yy1); } getcadre(xx0, xx1, yy0, yy1); + if (verbosity > 2) + cout << " getcadre bb= xmin =" << xx0 << ", max =" << xx1 << ", ymin = " << yy0 << ", ymax = " << yy1 << endl; const R ccoeff = coeff; bool plotting = true; // drawing part ------------------------------ @@ -5621,6 +5623,7 @@ AnyType ClearReturn(Stack stack, const AnyType &a) { return m; } template long get_n(KN * p){ return p->N();}// +template long get__n(KN_ p){ return p.N();}// template< class R > void DclTypeMatrix( ) { @@ -5828,7 +5831,37 @@ R3 toR3(const R& a,const R& b,const R &c){return R3(a,b,c);} R3 toR3(const R3& a,const R3& b){return R3(a,b);} R3 NElement(lgBoundaryEdge const & a) { return R3(a.NBoundaryElement()); }// add Jan 2022 +template +void Add_u_init_array(int ii=0) +{ + typedef Complex C; + typedef double R; + + TheOperators->Add( + "<-", + new init_FE_eqarray< FFset3< R, vfes, RNM_VirtualMatrix< R >::plusAx > >(10), + new init_FE_eqarray< FFset3< R, vfes, RNM_VirtualMatrix< R >::solveAxeqb > >(10), + new init_FE_eqarray< FFset3< R, vfes, RNM_VirtualMatrix< R >::solveAtxeqb > >(10), + new init_FE_eqarray< FFset3< R, vfes, RNM_VirtualMatrix< R >::plusAtx > >(10), + new init_FE_eqarray< FF_L_args< R, vfes, Call_FormLinear< v_fes > > >(10) + + ); + TheOperators->Add( + "<-", + new init_FE_eqarray< FFset3< C, vfes, RNM_VirtualMatrix< C >::plusAx > >(10), + new init_FE_eqarray< FFset3< C, vfes, RNM_VirtualMatrix< C >::solveAxeqb > >(10), + new init_FE_eqarray< FFset3< C, vfes, RNM_VirtualMatrix< C >::solveAtxeqb > >(10), + new init_FE_eqarray< FFset3< C, vfes, RNM_VirtualMatrix< C >::plusAtx > >(10) + + ); + + if( ii) + TheOperators->Add( + "<-", + new init_FE_eqarray< FFset3< R, vfes, KN_< R > > >(10), + new init_FE_eqarray< FFset3< C, vfes, KN_< C > > >(10)); +} void init_lgfem( ) { if (verbosity && (mpirank == 0)) cout << "lg_fem "; #ifdef HAVE_CADNA @@ -6551,7 +6584,8 @@ void init_lgfem( ) { new init_FE_eqarray< FF_L_args< Complex, v_fes3, Call_FormLinear< v_fes3 > > >( 10) // 3D volume ); - + Add_u_init_array(1); + Add_u_init_array(1); // Fin // Fin TheOperators->Add( diff --git a/src/fflib/lgmat.cpp b/src/fflib/lgmat.cpp index bac584750..cce234b6b 100644 --- a/src/fflib/lgmat.cpp +++ b/src/fflib/lgmat.cpp @@ -3692,6 +3692,43 @@ MatriceMorse * removeHalf(MatriceMorse & A,long half,double tol) return r; } + +template + KN > *Matrix2KKN(KN > *pij,Matrice_Creuse * pA,bool trans) + { + ffassert(pij); + if(init) pij->init(); + MatriceCreuse * pa=pA->A; + MatriceMorse *pma= dynamic_cast* > (pa); + if( pma ) { + long n = trans? pma->n:pma->m; + pij->resize(n); + KN > & ij=*pij; + int * pp, *pj; + R * pc; + if( trans) + pma->CSC(pp,pj,pc); + else + pma->CSR(pp,pj,pc); + + for( int i=0; i + KN > *Matrix2KKN(KN > *pij,Matrice_Creuse * pA) +{ return Matrix2KKN (pij,pA,false);} +template +KN > *Matrix2KKN(KN > *pij,Matrice_Creuse_Transpose pAt) +{ return Matrix2KKN(pij,pAt,true);} template newpMatrice_Creuse removeHalf(Stack stack,Matrice_Creuse *const & pA,long const & half,const double & tol) @@ -3991,6 +4028,15 @@ void init_lgmat() Global.Add("removeHalf", "(", new OneOperator3s_ ,Matrice_Creuse * ,long,double>(removeHalf)); Global.Add("removeHalf", "(", new OneOperator4s_ * ,Matrice_Creuse * ,long,double>(removeHalf)); + TheOperators->Add("<-", + new OneOperator2 >*,KN >*,Matrice_Creuse * >(Matrix2KKN )); + TheOperators->Add("<-", + new OneOperator2 >*,KN >*,Matrice_Creuse_Transpose >(Matrix2KKN )); + TheOperators->Add("=", + new OneOperator2 >*,KN >*,Matrice_Creuse * >(Matrix2KKN )); + TheOperators->Add("=", + new OneOperator2 >*,KN >*,Matrice_Creuse_Transpose >(Matrix2KKN )); + } int Data_Sparse_Solver_version() { return VDATASPARSESOLVER;} diff --git a/src/fflib/lgmesh3.cpp b/src/fflib/lgmesh3.cpp index d03dd66a2..4083e164c 100644 --- a/src/fflib/lgmesh3.cpp +++ b/src/fflib/lgmesh3.cpp @@ -2494,6 +2494,50 @@ class Op4_pfL2K : public quad_function *,int>,R,R,R,K> { pu }; }; +template +class Op3_pfL2K : public ternary_function *,int>,R,R,K> { public: + + + class Op : public E_F0mps { public: + Expression a,b,c; + Op(Expression aa,Expression bb,Expression cc) + : a(aa),b(bb),c(cc){} + AnyType operator()(Stack s) const + { + + R xx(GetAny((*b)(s))); + R yy(GetAny((*c)(s))); + R zz=R(); + MeshPoint & mp = *MeshPointStack(s),mps=mp; + mp.set(xx,yy,zz); + AnyType ret = pfLr2R(s,(*a)(s)); + mp=mps; + return ret;} + + }; +}; +template +class Op2_pfL2K : public binary_function *,int>,R,K> { public: + + + class Op : public E_F0mps { public: + Expression a,b; + Op(Expression aa,Expression bb) + : a(aa),b(bb){} + AnyType operator()(Stack s) const + { + + R xx(GetAny((*b)(s))); + R yy=R(); + R zz=R(); + MeshPoint & mp = *MeshPointStack(s),mps=mp; + mp.set(xx,yy,zz); + AnyType ret = pfLr2R(s,(*a)(s)); + mp=mps; + return ret;} + + }; +}; template @@ -2894,7 +2938,11 @@ void init_lgmesh3() { Add("[]",".",new OneOperator1 *,pfLc>(pf3r2vect)); Add("(","",new OneQuadOperator,Op4_pfL2K::Op> ); Add("(","",new OneQuadOperator,Op4_pfL2K::Op> ); - + Add("(","",new OneTernaryOperator,Op3_pfL2K::Op> ); + Add("(","",new OneTernaryOperator,Op3_pfL2K::Op> ); + Add("(","",new OneBinaryOperatorCode,Op2_pfL2K::Op> ); + Add("(","",new OneBinaryOperatorCode,Op2_pfL2K::Op> ); + Add("(","",new OneQuadOperator,Op4_K2R::Op> ); // Add("(","",new OneTernaryOperator,Op3_K2R::Op> ); // FH stupide Add("(","",new OneQuadOperator,Op4_K2R::Op> ); diff --git a/src/fflib/problem.cpp b/src/fflib/problem.cpp index 5180c2ceb..bb559d892 100755 --- a/src/fflib/problem.cpp +++ b/src/fflib/problem.cpp @@ -223,7 +223,7 @@ template pair_stack_double * bs=static_cast(bstack); Stack stack= bs->first; double binside = *bs->second; // truc FH pour fluide de grad2 (decentrage bizard) - ffassert(mat.onFace); // Finite Volume or discontinuous Galerkine + ffassert(mat.onFace); // Finite Volume or discontinuous Galerkin ffassert(ie>=0 && ie < 4); // int on Face MeshPoint mp= *MeshPointStack(stack); R ** copt = Stack_Ptr(stack,ElemMatPtrOffset); @@ -445,7 +445,7 @@ template pair_stack_double * bs=static_cast(bstack); Stack stack= bs->first; double binside = *bs->second; // truc FH pour fluide de grad2 (decentrage bizard) - assert(mat.onFace); // Finite Volume or discontinuous Galerkine + assert(mat.onFace); // Finite Volume or discontinuous Galerkin assert(ie>=0 && ie < 3); // int on edge MeshPoint mp= *MeshPointStack(stack); R ** copt = Stack_Ptr(stack,ElemMatPtrOffset); @@ -612,7 +612,7 @@ template pair_stack_double * bs=static_cast(bstack); Stack stack= bs->first; double binside = *bs->second; // truc FH pour fluide de grad2 (decentrage bizard) - assert(mat.onFace); // Finite Volume or discontinous Galerkine + assert(mat.onFace); // Finite Volume or discontinuous Galerkin assert(ie>=0 && ie < 3); // int on edge MeshPoint mp= *MeshPointStack(stack); R ** copt = Stack_Ptr(stack,ElemMatPtrOffset); @@ -9450,6 +9450,19 @@ void Element_rhs(const Mesh3 & ThI,const Mesh3::Element & KI, const FESpace3 & ffassert(0); } +// creating an instance of Element_rhs +// case 3D + +template +void Element_rhs3(Expression const * const mapt,const Mesh3 & ThI,const Mesh3::Element & KI, const FESpace3 & Vh, + const LOperaD &Op,double * p,void * vstack,KN_ & B, + const GQuadratureFormular & FI = QuadratureFormular_Tet_5,bool alledges=false,bool intmortar=false, + R3 *Q=0,int optim=1) +{ + // add mapt case in 3d here !!!!!!!! FH + ffassert(0); +} + // generic template for AssembleVarForm template bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESpace2 & Vh,bool sym, @@ -10635,7 +10648,9 @@ bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESp else if (CDomainOfIntegration::intallVFedges==kind) cout << " -- boundary int all VF edges nQP: ("<< FIT.n << ")," ; else cout << " -- int 3d (nQP: "<< FIV.n << " ) in " ; } - if( di.withmap()) { ExecError(" no map in the case (5)??");} + Expression const * const mapt=*di.mapt?di.mapt:0 ; + + // if( di.withmap()) { ExecError(" no map in the case (5)??");} // if(di.islevelset()) InternalError("So no levelset integration type on this case (3)"); if(di.islevelset() && (CDomainOfIntegration::int2d!=kind) && (CDomainOfIntegration::int3d!=kind) ) @@ -10713,6 +10728,7 @@ bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESp if(VF) InternalError(" no jump or average in int1d of RHS"); if(di.islevelset()) // init on level set (of RHS) { + ffassert(mapt==0); double uset = HUGE_VAL; R3 Q[4]; KN phi(ThI.nv);phi=uset; @@ -10780,8 +10796,11 @@ bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESp int ie,i =ThI.BoundaryElement(e,ie); if ( sameMesh) Element_rhs(Vh[i],ie,Th.be(e).lab,*l->l,buf,stack,*B,FIT,false,useopt); + else if(!mapt) + Element_rhs(ThI,ThI[i],Vh,ie,Th.be(e).lab,*l->l,buf,stack,*B,FIT,false,useopt); else - Element_rhs(ThI,ThI[i],Vh,ie,Th.be(e).lab,*l->l,buf,stack,*B,FIT,false,useopt); + ffassert(0); // DO do FH... + if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006 clean Ptr } } @@ -10834,8 +10853,11 @@ bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESp { if ( sameMesh ) Element_rhs(Vh[t],*l->l,buf,stack,*B,FIV,useopt); - else + else if(!mapt) Element_rhs(ThI,ThI[t],Vh,*l->l,buf,stack,*B,FIV,useopt); + else ffassert(0); + // Element_rhs3(mapt,ThI,ThI[t],Vh,*l->l,buf,stack,*B,FIV,useopt); + if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006 clean Ptr } @@ -10851,13 +10873,17 @@ bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESp { if ( sameMesh ) Element_rhs(Vh[i],*l->l,buf,stack,*B,FIV,useopt); - else + else if(!mapt) Element_rhs(ThI,ThI[i],Vh,*l->l,buf,stack,*B,FIV,useopt); + else + Element_rhs3(mapt,ThI,ThI[i],Vh,*l->l,buf,stack,*B,FIV,useopt);// to do !!!!! + // Element_rhs(mapt,ThI,ThI[i],Vh,*l->l,buf,stack,*B,FIV,useopt); + if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006 clean Ptr }} } else if(kind==CDomainOfIntegration::intallfaces ) { - + assert(mapt==0); if(VF) { if ( !sameMesh) InternalError(" no jump or average in intallfaces for RHS in not samemesh"); @@ -10894,9 +10920,12 @@ bool AssembleVarForm(Stack stack,const MMesh & Th,const FESpace1 & Uh,const FESp // if face on bord get the lab ??? if ( sameMesh) Element_rhs(Vh[i],ie,lab,*l->l,buf,stack,*B,FIT,false,useopt); - else + else if(!mapt) Element_rhs(ThI,ThI[i],Vh,ie,lab,*l->l,buf,stack,*B,FIT,false,useopt); - if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006 clean Ptr + else + ffassert(mapt==0);// + //Element_rhs(mapt,ThI,ThI[i],Vh,ie,lab,*l->l,buf,stack,*B,FIT,false,useopt); + if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006 clean Ptr } } @@ -11981,7 +12010,7 @@ AnyType Problem::eval(Stack stack,Data * data,CountPointer100000) )) { cout << " ######## We forget of deleting " << NbPtr + if (verbosity && (NbPtr || (stu1>100000) )) { cout << " ######## unfreed pointers " << NbPtr << " Nb pointer, " << lg1-lg0 << "Bytes " << " , mpirank " << mpirank << ", memory leak ="<< stu1 << endl;} return 0;} ; diff --git a/src/mpi/parallelempi.cpp b/src/mpi/parallelempi.cpp index 57b76b070..61dcef006 100644 --- a/src/mpi/parallelempi.cpp +++ b/src/mpi/parallelempi.cpp @@ -114,6 +114,13 @@ class DoOnWaitMPI_Request; map ToDoOnWaitMPI_Request; void GetPendingWait() ; +// change the priority to remove ambiguity between KN * and KN_ +template<> struct SameType *> { static const int OK=100;};// +template<> struct SameType *> { static const int OK=100;};// +template<> struct SameType *> { static const int OK=100;};// +template<> struct SameType *> { static const int OK=100;};// +template<> struct SameType *> { static const int OK=100;};// +template<> struct SameType *> { static const int OK=100;};// template struct MPI_TYPE { static MPI_Datatype TYPE(){ return MPI_BYTE; } }; template<> struct MPI_TYPE { static MPI_Datatype TYPE(){ return MPI_LONG; } }; @@ -323,6 +330,18 @@ struct MPIrank { return WSend((R *)a, n, who, MPI_TAG* >::TAG, comm, rq); } + template + long Send (const KN_ &aa) const { + const KN_ & a = aa; + ffassert(&a); + int n = a.N(); + CheckContigueKN(aa); + if (verbosity > 99) + cout << " .... _ " << who << " >> " << &a << " " << a.N() << " " << a.step << " " << MPI_TAG* >::TAG + << " from " << mpirank << " " << (R *)a << endl; + return WSend((R *)a, n, who, MPI_TAG* >::TAG, comm, rq); + } + template const MPIrank & Bcast (const KN &a) const { assert(&a); @@ -340,12 +359,12 @@ struct MPIrank { assert(&a); CheckContigueKNM(a); - if (verbosity > 99) + if (verbosity > 399) cout << " ---- " << who << " >> " << &a << " " << a.N() << "x" << a.M() << " " << MPI_TAG* >::TAG << " from " << mpirank << " " << (R *)a << endl; int n = a.N()*a.M(); long ll = WRecv((R *)a, n, who, MPI_TAG* >::TAG, comm, rq); - if (verbosity > 99) + if (verbosity > 399) cout << " ++++ " << who << " >> " << &a << " " << a.N() << "x" << a.M() << " " << MPI_TAG* >::TAG << " from " << mpirank << " " << (R *)a << endl; ffassert(a.N()*a.M() == n); @@ -358,11 +377,22 @@ struct MPIrank { ffassert(&a); int n = a.N()*a.M(); CheckContigueKNM(*aa); - if (verbosity > 99) + if (verbosity > 399) cout << " .... " << who << " >> " << &a << " " << a.N() << "x" << a.M() << " " << a.step << " " << MPI_TAG* >::TAG << " from " << mpirank << " " << (R *)a << endl; return WSend((R *)a, n, who, MPI_TAG* >::TAG, comm, rq); } + template + long Send (const KNM_ &aa) const { + const KNM_ & a = aa; + ffassert(&a); + int n = a.N()*a.M(); + CheckContigueKNM(aa); + if (verbosity > 399) + cout << " .... _ " << who << " >> " << &a << " " << a.N() << "x" << a.M() << " " << a.step << " " << MPI_TAG* >::TAG + << " from " << mpirank << " " << (R *)a << endl; + return WSend((R *)a, n, who, MPI_TAG* >::TAG, comm, rq); + } template const MPIrank & Bcast (const KNM &a) const { @@ -2754,6 +2784,12 @@ void f_init_lgparallele() Global.Add("Send","(", new OneBinaryOperator *> >); Global.Add("Send","(", new OneBinaryOperator *> >); Global.Add("Send","(", new OneBinaryOperator *> >); + Global.Add("Send","(", new OneBinaryOperator> >); + Global.Add("Send","(", new OneBinaryOperator > >); + Global.Add("Send","(", new OneBinaryOperator > >); + Global.Add("Send","(", new OneBinaryOperator > >); + Global.Add("Send","(", new OneBinaryOperator > >); + Global.Add("Send","(", new OneBinaryOperator > >); Global.Add("Send","(", new OneBinaryOperator >); Global.Add("Send","(", new OneBinaryOperator >); Global.Add("Send","(", new OneBinaryOperator >); @@ -2765,11 +2801,17 @@ void f_init_lgparallele() Global.Add("Isend","(", new OneBinaryOperator >); Global.Add("Isend","(", new OneBinaryOperator >); Global.Add("Isend","(", new OneBinaryOperator *> >); + Global.Add("Isend","(", new OneBinaryOperator > >); Global.Add("Isend","(", new OneBinaryOperator *> >); + Global.Add("Isend","(", new OneBinaryOperator> >); Global.Add("Isend","(", new OneBinaryOperator *> >); + Global.Add("Isend","(", new OneBinaryOperator> >); Global.Add("Isend","(", new OneBinaryOperator *> >); Global.Add("Isend","(", new OneBinaryOperator *> >); Global.Add("Isend","(", new OneBinaryOperator *> >); + Global.Add("Isend","(", new OneBinaryOperator> >); + Global.Add("Isend","(", new OneBinaryOperator> >); + Global.Add("Isend","(", new OneBinaryOperator> >); Global.Add("Isend","(", new OneBinaryOperator >); Global.Add("Isend","(", new OneBinaryOperator >); Global.Add("Isend","(", new OneBinaryOperator >);