diff --git a/CMakeLists.txt b/CMakeLists.txt index 134c6f5..7bf6924 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,7 +25,7 @@ cmake_minimum_required(VERSION 3.21.0) project( feelpp-project - VERSION 3.0.0 ) + VERSION 3.0.1 ) set(EXTRA_VERSION "") set(PROJECT_SHORTNAME "fp") diff --git a/README.adoc b/README.adoc index f9704bf..6894fc2 100644 --- a/README.adoc +++ b/README.adoc @@ -2,7 +2,7 @@ :cpp: C++ :project: feelpp-project -= {feelpp} Project Template += {feelpp} Template Project Christophe Prud'homme v2: diff --git a/README.md b/README.md index f9704bf..3c6aa03 100644 --- a/README.md +++ b/README.md @@ -1,38 +1,31 @@ -:feelpp: Feel++ -:cpp: C++ -:project: feelpp-project +# Feel++ Template Project +Author: Christophe Prud'homme ![GitHub](https://github.com/prudhomm) +Version: v2 -= {feelpp} Project Template -Christophe Prud'homme -v2: +![CI](https://github.com/feelpp/feelpp-project/workflows/CI/badge.svg) -image:https://github.com/feelpp/feelpp-project/workflows/CI/badge.svg[CI] +This repository provides a basic starting point for a Feel++ application including: -This repository provides a basic starting point for a {feelpp} application including: +- [x] Feel++ applications in C++ to use Feel++ and Feel++ toolboxes in `src` +- [x] Documentation using asciidoc and antora +- [x] Python Feel++ notebooks that can be downloaded from the documentation +- [x] Continuous integration including tests for the C++ applications +- [x] Docker image generation for the project -- [x] {feelpp} applications in {cpp} to use {feelpp} and {feelpp} toolboxes in `src` -- [x] documentation using asciidoc and antora -- [x] python {feelpp} notebooks that can be downloaded from the documentation -- [x] continuous integration including tests for the {cpp} applications -- [x] docker image generation for the project +The documentation for feelpp-project is available [here](https://feelpp.github.io/feelpp-project), and you can build on it for your project by enabling the [github pages](https://docs.github.com/en/pages) for your repository. -The documentation for feelpp-project is available at link:https://feelpp.github.io/feelpp-project[here] and you can build on it for your project by enabling the link:https://docs.github.com/en/pages[github pages] for your repository. +## Renaming the project -== Renaming the project +By default, the project is named `feelpp-project` if you cloned the repository `feelpp/feelpp-project`. However, if you used the previous repository as a template, then the project is renamed using the name of the repository using the script `rename.sh` at the initialization of the repository. If the name does not suit you, you can change it again using the script `rename.sh` and providing the new name as an argument. -By default the project is named `feelpp-project` if you cloned the repository `feelpp/feelpp-project`. -However if you used the previous repository as a template, then the project is renamed using the name of the repository using the script `rename.sh` at the initialization of the repository. -If the name does not suit you, you can change it again using the script `rename.sh` and providing the new name as argument. +> **Warning**: The script `rename.sh` will rename the project. However, some URLs might not be set properly if you rename the project yourself. You need to check the following files: `docs/site.yml` and `docs/package.json` and fix the URLs after the rename process is done. -WARNING: the script `rename.sh` will rename the project however some url might be set properly if you rename the project yourself. You need to check the following files: `docs/site.yml` and `docs/package.json` and fix the urls after the rename process is done. +## Updating the project version -== Updating the {project} version +The version of the project is defined in the files `CMakeLists.txt`, `docs/antora.yml`, and `docs/package.json`. You need to update with the same version in all files. -The version of the project is defined in the files `CMakeLists.txt`, `docs/antora.yml` and `docs/package.json`. -You need to update with the same version in all files. +## Release process -== Release process - -- [x] update the version in CMakeLists.txt -- [x] update the version in docs/antora.yml -- [x] commit the changes with the message "Release vx.y.z". At this point the CI will generate the docker image and push it to docker hub \ No newline at end of file +- [x] Update the version in `CMakeLists.txt` +- [x] Update the version in `docs/antora.yml` +- [x] Commit the changes with the message "Release vx.y.z". At this point, the CI will generate the docker image and push it to Docker Hub. diff --git a/docs/modules/ROOT/pages/examples/fin.adoc b/docs/modules/ROOT/pages/examples/fin.adoc index 8f258de..331e217 100644 --- a/docs/modules/ROOT/pages/examples/fin.adoc +++ b/docs/modules/ROOT/pages/examples/fin.adoc @@ -6,6 +6,8 @@ Christophe Prud'homme :page-illustration: fin2d-4.png :description: We simulate the heat transfer in a thermal fin +== Problem statement + We consider the problem of designing a thermal fin to effectively remove heat from a surface. The two-dimensional fin, shown in Figure below, consists of a vertical central post and four horizontal subfins; the fin conducts heat from a prescribed uniform flux source at the root, stem:[\Gamma_{\text {root }}], through the large-surface-area subfins to surrounding flowing air. The fin is characterized by a five-component parameter vector, or input, stem:[\mu_{=}\left(\mu_1, \mu_2, \ldots, \mu_5\right)], where stem:[\mu_i=k^i, i=1, \ldots, 4], and stem:[\mu_5=\mathrm{Bi} ; \mu] may take on any value in a specified design set stem:[D \subset \mathbb{R}^5]. .Geometry of 2D thermal fin @@ -16,9 +18,12 @@ Here stem:[k^i] is the thermal conductivity of the ith subfin (normalized relati We are interested in the design of this thermal fin, and we thus need to look at certain outputs or cost-functionals of the temperature as a function of stem:[\mu]. We choose for our output stem:[T_{\text {root }}], the average steady-state temperature of the fin root normalized by the prescribed heat flux into the fin root. The particular output chosen relates directly to the cooling efficiency of the fin lower values of stem:[T_{\text {root }}] imply better thermal performance. The steadystate temperature distribution within the fin, stem:[u(\mu)], is governed by the elliptic partial differential equation [stem] ++++ --k^i \Delta u^i=0 \text { in } \Omega^i, i=0, \ldots, 4, +\rho_i C_i \frac{\partial u^i}{\partial t} -k^i \Delta u^i=0 \text { in } \Omega^i, i=0, \ldots, 4, ++++ -where stem:[\Delta] is the Laplacian operator, and stem:[u_i] refers to the restriction of stem:[u] to stem:[\Omega^i] . Here stem:[\Omega^i] is the region of the fin with conductivity stem:[k^i, i=0, \ldots, 4: \Omega^0] is thus the central post, and stem:[\Omega^i, i=1, \ldots, 4], corresponds to the four subfins. The entire fin domain is denoted stem:[\Omega\left(\bar{\Omega}=\cup_{i=0}^4 \bar{\Omega}^i\right)]; the boundary stem:[\Omega] is denoted stem:[\Gamma]. We must also ensure continuity of temperature and heat flux at the conductivity discontinuity interfaces stem:[\Gamma_{\text {int }}^i \equiv \partial \Omega^0 \cap \partial \Omega^i, i=1, \ldots, 4], where stem:[\partial \Omega^i] denotes the boundary of stem:[\Omega^i], we have on stem:[\Gamma_{\text {int }}^i i=1, \ldots, 4] : +where stem:[\Delta] is the Laplacian operator, and stem:[u_i] refers to the restriction of stem:[u] to stem:[\Omega^i]. +Here stem:[\Omega^i] is the region of the fin with conductivity stem:[k^i, i=0, \ldots, 4] and volumetric heat capacity stem:[(\rho C)_i, i=0, \cdots, 4]: stem:[\Omega^0] is thus the central post, and stem:[\Omega^i, i=1, \ldots, 4], corresponds to the four subfins. + +The entire fin domain is denoted stem:[\Omega\left(\bar{\Omega}=\cup_{i=0}^4 \bar{\Omega}^i\right)]; the boundary stem:[\Omega] is denoted stem:[\Gamma]. We must also ensure continuity of temperature and heat flux at the conductivity discontinuity interfaces stem:[\Gamma_{\text {int }}^i \equiv \partial \Omega^0 \cap \partial \Omega^i, i=1, \ldots, 4], where stem:[\partial \Omega^i] denotes the boundary of stem:[\Omega^i], we have on stem:[\Gamma_{\text {int }}^i i=1, \ldots, 4] : [stem] ++++ \begin{aligned} diff --git a/package-lock.json b/package-lock.json index a50b66b..b10e762 100644 --- a/package-lock.json +++ b/package-lock.json @@ -11,6 +11,7 @@ "dependencies": { "@antora/cli": "^3.1", "@antora/collector-extension": "^1.0.0-alpha.3", + "@antora/site-generator": "^3.1", "@antora/site-generator-default": "^3.1", "@asciidoctor/core": "^2.2.6", "@feelpp/asciidoctor-extensions": "1.0.0-rc.8", diff --git a/package.json b/package.json index 56df588..002793c 100644 --- a/package.json +++ b/package.json @@ -29,6 +29,7 @@ "dependencies": { "@antora/cli": "^3.1", "@antora/collector-extension": "^1.0.0-alpha.3", + "@antora/site-generator": "^3.1", "@antora/site-generator-default": "^3.1", "@asciidoctor/core": "^2.2.6", "asciidoctor": "^2.2.6", diff --git a/pyproject.toml b/pyproject.toml index ef2d81e..ae49b4c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,8 +5,8 @@ build-backend = "scikit_build_core.build" [project] name = "feelpp-project" -version = "3.0.0" -description="Feelpp Template Project" +version = "3.0.1" +description="Feel++ Template Project" readme = "README.md" authors = [ { name = "Christophe Prud'homme", email = "christophe.prudhomme@cemosis.fr" }, diff --git a/rename.sh b/rename.sh index 1c53379..494ee8d 100644 --- a/rename.sh +++ b/rename.sh @@ -70,6 +70,8 @@ else projecttitle=$(create_title ${projectname}) shortprojectname=$( create_short_name $projectname ) appname="app" + projectname_="${projectname//-/_}" + fi echo " project name: $projectname" @@ -89,27 +91,31 @@ if [ "$1" == "-r" ]; then fi export projectname +export projectname_ export projecttitle export shortprojectname export appname -files=( "README.adoc" "CMakeLists.txt" "src/CMakeLists.txt" +files=( "README.adoc" "README.md" "CMakeLists.txt" "src/CMakeLists.txt" "src/.tests.laplacian" "src/.tests.toolbox" "site.yml" "docs/antora.yml" "package.json" "package-lock.json" + "pyproject.toml" "docs/modules/ROOT/pages/index.adoc" ".github/workflows/ci.yml" ) for i in "${files[@]}" do echo "processing renaming in $i ...." perl -077pi.bak -e 's/feelpp-project/$ENV{'projectname'}/sg' $i + perl -077pi.bak -e 's/feelpp_project/$ENV{'projectname_'}/sg' $i perl -077pi.bak -e 's/Feel\+\+ Template Project/$ENV{'projecttitle'}/sg' $i perl -077pi.bak -e 's/\"fp\"/"$ENV{'shortprojectname'}"/sg' $i perl -077pi.bak -e 's/myapp/$ENV{'appname'}/sg' $i done +git mv src/feelpp_project src/$projectname_ diff --git a/requirements.txt b/requirements.txt index 7721efc..5bdb4cf 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,6 +3,7 @@ build scikit-build >= 0.8.1 scikit_build_core pyproject-metadata +ipykernel numpy scipy pandas diff --git a/src/_laplacian.cpp b/src/_laplacian.cpp index ee7f19c..d9fdc09 100644 --- a/src/_laplacian.cpp +++ b/src/_laplacian.cpp @@ -90,7 +90,7 @@ PYBIND11_MODULE(_laplacian, m ) if (import_mpi4py()<0) return ; m.doc() = fmt::format("Python bindings for Laplacian class" ); // Optional module docstring laplacian_inst<2,1>(m); -// laplacian_inst<2,2>(m); + laplacian_inst<2,2>(m); // laplacian_inst<2,3>(m); // laplacian_inst<3,1>(m); // laplacian_inst<3,2>(m); diff --git a/src/feelpp_project/__init__.py b/src/feelpp_project/__init__.py index 23d63ad..3f2ff2d 100644 --- a/src/feelpp_project/__init__.py +++ b/src/feelpp_project/__init__.py @@ -1,4 +1,4 @@ -print(f"Welcome in feelpp_project!") + diff --git a/src/feelpp_project/laplacian.py b/src/feelpp_project/laplacian.py index 2304032..d360442 100644 --- a/src/feelpp_project/laplacian.py +++ b/src/feelpp_project/laplacian.py @@ -4,6 +4,7 @@ _laps = { 'laplacian(2,1)': Laplacian2DP1, + 'laplacian(2,2)': Laplacian2DP2, } diff --git a/src/src/.tests.laplacian b/src/src/.tests.laplacian deleted file mode 100644 index 8db35e6..0000000 --- a/src/src/.tests.laplacian +++ /dev/null @@ -1 +0,0 @@ -laplacian-1 --config-file ${CMAKE_CURRENT_SOURCE_DIR}/cases/laplacian/fin/fin1/fin2d.cfg \ No newline at end of file diff --git a/src/src/.tests.toolbox b/src/src/.tests.toolbox deleted file mode 100644 index 0f090ec..0000000 --- a/src/src/.tests.toolbox +++ /dev/null @@ -1 +0,0 @@ -toolbox-1 --config-file ${CMAKE_CURRENT_SOURCE_DIR}/cases/electric/2d.cfg diff --git a/src/src/cases/electric/2d.cfg b/src/src/cases/electric/2d.cfg deleted file mode 100644 index 555064e..0000000 --- a/src/src/cases/electric/2d.cfg +++ /dev/null @@ -1,15 +0,0 @@ -directory=feelpp-project/update - -myexpr=-(x*x+y*y)*sin(x*y):x:y - -[electric] -mesh.filename=$cfgdir/2d.geo -filename=$cfgdir/2d.json - -verbose_solvertimer=1 -pc-type=gamg -ksp-monitor=1 - -[electric.gmsh] -hsize=0.1 -partition=1 diff --git a/src/src/cases/electric/2d.geo b/src/src/cases/electric/2d.geo deleted file mode 100644 index ada4f9f..0000000 --- a/src/src/cases/electric/2d.geo +++ /dev/null @@ -1,17 +0,0 @@ -h=0.1; - -Point(1) = {0,0,0,h}; -Point(2) = {1,0,0,h}; -Point(3) = {0,1,0,h}; -Point(4) = {1,1,0,h}; - -Line(1) = {1,2}; -Line(2) = {2,4}; -Line(3) = {4,3}; -Line(4) = {3,1}; - -Line Loop(1) = {1,2,3,4}; -Plane Surface(1) = {1}; - -Physical Surface("omega") = {1}; -Physical Line("gamma") = {1,2,3,4}; diff --git a/src/src/cases/electric/2d.json b/src/src/cases/electric/2d.json deleted file mode 100644 index f278893..0000000 --- a/src/src/cases/electric/2d.json +++ /dev/null @@ -1,52 +0,0 @@ -// -*- mode: javascript -*- -{ - "Name":"update", - "ShortName":"update", - "Materials": - { - "omega": - { - "sigma":"-1" - } - }, - "BoundaryConditions": - { - "electric-potential": - { - "Dirichlet": - { - "gamma": - { - "expr":"sin(x*y):x:y" - } - } - } - }, - "PostProcess": - { - "Exports": - { - "fields":["electric-potential","electric-field","pid"] - }, - "Measures": - { - "Norm": - { - "potential-error": - { - "type": "L2-error", - "field": "electric-potential", - "solution": "sin(x*y):x:y", - "markers": "omega" - }, - "field-error": - { - "type": "L2-error", - "field": "electric-field", - "solution": "{y*cos(x*y),x*cos(x*y)}:x:y", - "markers": "omega" - } - } - } - } -} diff --git a/src/src/cases/laplacian/fin/fin.geo b/src/src/cases/laplacian/fin/fin.geo deleted file mode 100644 index 5519a9f..0000000 --- a/src/src/cases/laplacian/fin/fin.geo +++ /dev/null @@ -1,105 +0,0 @@ -SetFactory("OpenCASCADE"); -//+ -h = 0.05; - -Nfins=4; -Lfins=2.5; -dim=2; - -N = DefineNumber[ Nfins, Name "Parameters/N" ]; -L = DefineNumber[ Lfins, Name "Parameters/L" ]; -t = DefineNumber[ 0.25, Name "Parameters/t" ]; -d = DefineNumber[ 0.75, Name "Parameters/d" ]; -hmax=0.5; -Mesh.CharacteristicLengthMax = hmax; - -If ( dim == 2 ) - -//+ -Rectangle(1) = {0, 0, 0, 1, N*(d+t)+t, 0}; -For r In {1:N} - Printf("fin = %g",r); - Rectangle(r+1) = {-L, r*(d+t), 0, 2*L+1, t, 0}; -EndFor - -S[]=BooleanFragments{ Surface{1}; Delete; }{ Surface{2:N+1}; Delete;}; -Characteristic Length{ PointsOf{ Surface{ : }; } } = h; - -If ( Nfins==1 ) - -EndIf -//+ -// Physical Surface (Sprintf("fin-%g",r)) = r+1; - -Physical Surface("Post") = {1:2*N}; -finid = 1; -For r In {2*N+1:#S[]:2 } - Printf("Fin number %g %g ", r, finid); - Physical Surface(Sprintf("Fin_%g",finid)) = {r,r+1}; - finid += 1; -EndFor - -bdy[] = CombinedBoundary { Surface{:}; }; - -Physical Curve("Gamma_ext") = bdy[0]; - -For ii In { 1 : (#bdy[]-1) } - If (bdy[ii] != 4) - Printf("boundary number %g = %g", ii, Abs(bdy[ii])); - Physical Curve("Gamma_ext") += Abs(bdy[ii]); - Else - Physical Curve("Gamma_root") = bdy[ii]; - EndIf -EndFor - -// Nfins = 1 2 cavities -If ( Nfins == 1 ) - Physical Curve("Cavity_1_1") = {1}; - Physical Curve("Cavity_1_2") = {10}; - - Physical Curve("Cavity_2_1") = {3}; - Physical Curve("Cavity_2_2") = {13}; -EndIf -ElseIf ( dim == 3 ) - -//+ -Box(1) = {0, 0, 0, 1, 1, N*(d+t)+t}; -For r In {1:N} - Printf("fin = %g",r); - Box(r+1) = {-L, -L, r*(d+t), 2*L+1, 2*L+1, t}; -EndFor - -S[]=BooleanFragments{ Volume{1}; Delete; }{ Volume{2:N+1}; Delete;}; -Characteristic Length{ PointsOf{ Surface{ : }; } } = h; - -//+ -// Physical Volume (Sprintf("fin-%g",r)) = r+1; - -Physical Volume("Post") = {1:2*N}; -finid = 1; -For r In {2*N+1:#S[]:1 } - Printf("Fin number %g %g ", r, finid); - Physical Volume(Sprintf("Fin_%g",finid)) = { r }; - finid += 1; -EndFor - -bdy[] = CombinedBoundary { Volume{:}; }; - -Physical Surface("Gamma_ext") = bdy[0]; - -For ii In { 1 : (#bdy[]-1) } - If (bdy[ii] != 5) - Printf("boundary number %g = %g", ii, Abs(bdy[ii])); - Physical Surface("Gamma_ext") += Abs(bdy[ii]); - Else - Physical Surface("Gamma_root") = bdy[ii]; - EndIf -EndFor - -// Mesh 2; -// Show "*"; - - -EndIf -// Mesh 2; -// Show "*"; diff --git a/src/src/cases/laplacian/fin/fin1/fin2d.cfg b/src/src/cases/laplacian/fin/fin1/fin2d.cfg deleted file mode 100644 index 781299d..0000000 --- a/src/src/cases/laplacian/fin/fin1/fin2d.cfg +++ /dev/null @@ -1,10 +0,0 @@ -specs=$cfgdir/../fin2d.json - -[gmsh] -geo-variables-list=Nfins=1:dim=2 - -[bdf] -order=2 -time-step=0.5 -time-final=10 -steady=1 \ No newline at end of file diff --git a/src/src/cases/laplacian/fin/fin2d.json b/src/src/cases/laplacian/fin/fin2d.json deleted file mode 100644 index a8e2a3a..0000000 --- a/src/src/cases/laplacian/fin/fin2d.json +++ /dev/null @@ -1,102 +0,0 @@ -{ - "Name": "Thermalfin 2D", - "ShortName": "thermalfin2d", - "Models": { - "laplacian": { - "name": "omega", - "Materials": [ - "Post", - "Fin_1", - "Fin_2", - "Fin_3", - "Fin_4" - ] - } - }, - "Meshes": { - "laplacian": { - "Import": { - "filename": "$cfgdir/../fin.geo", - "partition": 0 - } - } - }, - "Spaces": { - "laplacian": { - "Domain": { - //"marker": ["Post"] - "levelset": "x > 0:x:y" - } - } - }, - "Time": - { - "laplacian" :{ - "steady": false, - "order" : 1, - "start": 0.0, - "end": 1.0, - "step": 0.1 - } - }, - "Materials": { - "Post": { - "k": "1", - "Cp": "1", - "rho": "1" - }, - "Fin_1": { - "k": "1", - "Cp": "1", - "rho": "1" - }, - "Fin_2": { - "k": "1", - "Cp": "1", - "rho": "1" - }, - "Fin_3": { - "k": "1", - "Cp": "1", - "rho": "1" - }, - "Fin_4": { - "k": "1", - "Cp": "1", - "rho": "1" - } - }, - "InitialConditions": { - "laplacian": { - "temperature": { - "Expression": { - "Tini": { - "markers": [ - "Fin_1", - "Fin_2", - "Fin_3", - "Fin_4", - "Post" - ], - "expr": "0" - } - } - } - } - }, - "BoundaryConditions": { - "laplacian": { - "flux": { - "Gamma_root": { - "expr": "1" - } - }, - "convective_laplacian_flux": { - "Gamma_ext": { - "h": "0.1", - "Text": "0" - } - } - } - } -} diff --git a/src/src/laplacian.py b/src/src/laplacian.py deleted file mode 100644 index 9da84af..0000000 --- a/src/src/laplacian.py +++ /dev/null @@ -1,287 +0,0 @@ -import numpy as np -import os -import pyvista as pv -from xvfbwrapper import Xvfb -import sys -import feelpp -import feelpp.toolboxes.core as tb -from feelpp.toolboxes.cfpdes import * -import pandas as pd -from xvfbwrapper import Xvfb -import pyvista as pv -import numpy as np -import plotly.express as px -from plotly.subplots import make_subplots -import itertools - -def generateGeometry(filename, dim=2, hsize=0.1): - """create gmsh mesh - - Args: - filename (str): name of the file - dim (int): dimension of the mesh - hsize (float): mesh size - """ - geo = """SetFactory("OpenCASCADE"); - h={}; - dim={}; - """.format(hsize, dim) - if dim == 2: - geo += """ - Rectangle(1) = {0, 0, 0, 1, 1, 0}; - Characteristic Length{ PointsOf{ Surface{1}; } } = h; - Physical Curve("Gamma_D") = {1,2,3,4}; - Physical Surface("Omega") = {1}; - """ - elif dim == 3: - geo += """ - Box(1) = {0, 0, 0, 1, 1, 1}; - Characteristic Length{ PointsOf{ Volume{1}; } } = h; - Physical Surface("Gamma_D") = {1,2,3,4,5,6}; - Physical Volume("Omega") = {1}; - """ - with open(filename, 'w') as f: - # Write the string to the file - f.write(geo) - - -def getMesh(filename, hsize=0.05, dim=2, verbose=False): - """create mesh - - Args: - filename (str): name of the file - hsize (float): mesh size - dim (int): dimension of the mesh - verbose (bool): verbose mode - """ - import os - for ext in [".msh", ".geo"]: - f = os.path.splitext(filename)[0]+ext - # print(f) - if os.path.exists(f): - os.remove(f) - if verbose: - print( - f"generate mesh {filename} with hsize={hsize} and dimension={dim}") - generateGeometry(filename=filename, dim=dim, hsize=hsize) - mesh = feelpp.load(feelpp.mesh(dim=dim, realdim=dim), filename, hsize) - return mesh - - -def laplacian(hsize, json, dim=2, verbose=False): - """solve the laplacian problem - - Args: - hsize (_type_): mesh size - json (_type_): json data - dim (int, optional): dimension. Defaults to 2. - verbose (bool, optional): verbosity level. Defaults to False. - - Returns: - dict: measures - """ - if verbose: - print(f"Solving the laplacian problem for hsize = {hsize}...") - laplacian = cfpdes(dim=dim, keyword=f"cfpdes-{dim}d") - laplacian.setMesh( - getMesh(f"omega-{dim}.geo", hsize=hsize, dim=dim, verbose=verbose)) - laplacian.setModelProperties(json) - laplacian.init(buildModelAlgebraicFactory=True) - laplacian.printAndSaveInfo() - laplacian.solve() - laplacian.exportResults() - measures = laplacian.postProcessMeasures().values() - - return measures - - -def pv_get_mesh(mesh_path): - reader = pv.get_reader(mesh_path) - mesh = reader.read() - return mesh - - -def pv_plot(mesh, field, clim=None, cmap='viridis', cpos='xy', show_scalar_bar=True, show_edges=True): - mesh.plot(scalars=field, clim=clim, cmap=cmap, cpos=cpos, - show_scalar_bar=show_scalar_bar, show_edges=show_edges) - - -def myplots(dim=2, field="cfpdes.laplace.u", factor=1, cmap='viridis'): - mesh = pv_get_mesh(f"cfpdes-{dim}d.exports/Export.case") # <4> - pv_plot(mesh, field) # <5> - pl = pv.Plotter() - contours = mesh[0].contour() - pl.add_mesh(mesh[0], opacity=0.85) - pl.add_mesh(contours, color="white", line_width=5, - render_lines_as_tubes=True) - pl.show() - if dim == 2: - warped = mesh[0].warp_by_scalar(field, factor=factor) - warped.plot(cmap=cmap, show_scalar_bar=False, show_edges=True) - else: - slices = mesh.slice_orthogonal(x=0.2, y=0.4,z=.6) - slices.plot() - - -def runLaplacianPk(df, model, verbose=False): - """generate the Pk case - - Args: - order (int, optional): order of the basis. Defaults to 1. - """ - meas = dict() - dim, order, json = model - for h in df['h']: - m = laplacian(hsize=h, json=json, dim=dim, verbose=verbose) - for norm in ['L2', 'H1']: - meas.setdefault(f'P{order}-Norm_laplace_{norm}-error', []) - meas[f'P{order}-Norm_laplace_{norm}-error'].append( - m.pop(f'Norm_laplace_{norm}-error')) - df = df.assign(**meas) - for norm in ['L2', 'H1']: - df[f'P{order}-laplace_{norm}-convergence-rate'] = np.log2(df[f'P{order}-Norm_laplace_{norm}-error'].shift( - ) / df[f'P{order}-Norm_laplace_{norm}-error']) / np.log2(df['h'].shift() / df['h']) - return df - - -def runConvergenceAnalysis(json, dim=2, hs=[0.1, 0.05, 0.025, 0.0125], orders=[1, 2], verbose=False): - df = pd.DataFrame({'h': hs}) - for order in orders: - df = runLaplacianPk(df=df, model=[dim, order, json( - dim=dim, order=order)], verbose=verbose) - print(df.to_markdown()) # <1> - return df - - -df = runConvergenceAnalysis(json=laplacian_json, dim=2, verbose=True) -df3d = runConvergenceAnalysis(json=laplacian_json, dim=3, hs=[ - 0.1, 0.05, 0.03], orders=[1], verbose=True) - - -def plot_convergence(df, dim, orders=[1, 2]): - fig = px.line(df, x="h", y=[f'P{order}-Norm_laplace_{norm}-error' for order, - norm in list(itertools.product(orders, ['L2', 'H1']))]) - fig.update_xaxes(title_text="h", type="log") - fig.update_yaxes(title_text="Error", type="log") - for order, norm in list(itertools.product(orders, ['L2', 'H1'])): - fig.update_traces(name=f'P{order} - {norm} error - {df[f"P{order}-laplace_{norm}-convergence-rate"].iloc[-1]:.2f}', selector=dict( - name=f'P{order}-Norm_laplace_{norm}-error')) - fig.update_layout( - title=f"Convergence rate for the {dim}D Laplacian problem", - autosize=False, - width=900, - height=900, - ) - fig.show() - -sys.argv = ["feelpp_cfpdes_poisson"] -e = feelpp.Environment(sys.argv, - opts=tb.toolboxes_options( - "coefficient-form-pdes", "cfpdes"), - config=feelpp.globalRepository("cfpdes-poisson-homogeneous-dirichlet")) - -# generate 2D abd 3D meshes -for dim in [2, 3]: - mesh = getMesh(f"omega-{dim}d.geo", hsize=0.1, dim=dim, verbose=True) - - -laplacian_json = lambda order, dim=2, name="u": { - "Name": "Laplacian", - "ShortName": "Laplacian", - "Models": - { - f"cfpdes-{dim}d": - { - "equations": "laplace" - }, - "laplace": { - "setup": { - "unknown": { - "basis": f"Pch{order}", - "name": f"{name}", - "symbol": "u" - }, - "coefficients": { - "c": "1", - - "f": "8*pi*pi*sin(2*pi*x)*sin(2*pi*y):x:y" if dim == 2 else "12*pi*pi*sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z):x:y:z" - } - } - } - }, - "Materials": - { - "Omega": - { - "markers": ["Omega"] - } - }, - "BoundaryConditions": - { - "laplace": - { - "Dirichlet": - { - "g": - { - "markers": ["Gamma_D"], - "expr": "0" - } - } - } - }, - "PostProcess": - { - f"cfpdes-{dim}d": - { - "Exports": - { - "fields": ["all"], - "expr": { - "u_exact": "sin(2*pi*x)*sin(2*pi*y):x:y" if dim == 2 else "sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z):x:y:z", - "grad_u_exact": "{2*pi*cos(2*pi*x)*sin(2*pi*y),2*pi*sin(2*pi*x)*cos(2*pi*y)}:x:y" if dim == 2 else "{2*pi*cos(2*pi*x)*sin(2*pi*y)*sin(2*pi*z),2*pi*sin(2*pi*x)*cos(2*pi*y)*sin(2*pi*z),2*pi*sin(2*pi*x)*sin(2*pi*y)*cos(2*pi*z)}:x:y:z" - } - }, - "Measures": - { - "Norm": - { - "laplace": - { - "type": ["L2-error", "H1-error"], - "field": f"laplace.{name}", - "solution": "sin(2*pi*x)*sin(2*pi*y):x:y" if dim == 2 else "sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z):x:y:z", - "grad_solution": "{2*pi*cos(2*pi*x)*sin(2*pi*y),2*pi*sin(2*pi*x)*cos(2*pi*y)}:x:y" if dim == 2 else "{2*pi*cos(2*pi*x)*sin(2*pi*y)*sin(2*pi*z),2*pi*sin(2*pi*x)*cos(2*pi*y)*sin(2*pi*z),2*pi*sin(2*pi*x)*sin(2*pi*y)*cos(2*pi*z)}:x:y:z", - "markers": "Omega", - "quad": 6 - } - } - } - } - } -} -# simulate the laplacian problem for 2D and 3D -for dim in [2, 3]: - with open(f'laplacian-{dim}d.json', 'w') as f: - # Write the string to the file - import json - f.write(json.dumps(laplacian_json(dim=dim, order=1), indent=1)) - # execute the laplacian problem using P1 basis on a mesh of the unit square of size 0.1 - laplacian(hsize=0.1, json=laplacian_json( - order=1, dim=dim), dim=dim, verbose=True) -# execute the laplacian problem using P2 basis on a mesh of the unit square of size 0.1 -# laplacian(hsize=0.025,json=laplacian_json(dim=2,order=2),dim=2,verbose=True) - -vdisplay = Xvfb() -vdisplay.start() - -pv.set_jupyter_backend('panel') - -myplots(dim=2,factor=0.5) - -myplots(dim=3,factor=0.5) - - -plot_convergence(df,dim=2) - -plot_convergence(df3d,dim=3,orders=[1]) \ No newline at end of file diff --git a/src/src/toolbox.cpp b/src/src/toolbox.cpp deleted file mode 100644 index 763fba4..0000000 --- a/src/src/toolbox.cpp +++ /dev/null @@ -1,53 +0,0 @@ -#include -#include - -int main(int argc, char** argv) -{ - using namespace Feel; - po::options_description myoptions("my options"); - myoptions.add( toolboxes_options("electric") ); - myoptions.add_options() - ( "myexpr", po::value()->default_value("3*x*y:x:y"), "my expression") - ; - - Environment env( _argc=argc, _argv=argv, - _desc=myoptions, - _about=about(_name="update_toolbox", - _author="Feel++ Consortium", - _email="feelpp-devel@feelpp.org")); - - typedef FeelModels::Electric< Simplex<2,1>, - Lagrange<1, Scalar,Continuous,PointSetFekete> > model_type; - std::shared_ptr electric( new model_type("electric") ); - // init the toolbox - electric->init(); - electric->printAndSaveInfo(); - - // create a lambda function - auto lambda = [&electric](FeelModels::ModelAlgebraic::DataUpdateLinear & data) { - // build each time, not just for the constant part - bool buildCstPart = data.buildCstPart(); - if( buildCstPart ) - return; - // retrieve matrix and vector already assemble - sparse_matrix_ptrtype& A = data.matrix(); - vector_ptrtype& F = data.rhs(); - // retrieve space and create a test function - auto Xh = electric->spaceElectricPotential(); - auto v = Xh->element(); - - // create a linear form from the vector on the function space - auto f = form1( _test=Xh, _vector=F); - // get an expression from the options (or from another equation) - auto e = expr(soption("myexpr")); - // add the term to the linear form - f += integrate( _range=elements(support(Xh)), _expr=inner(e, id(v)) ); - }; - // add the lambda function to the algebraic factory - electric->algebraicFactory()->addFunctionLinearAssembly(lambda); - // solve the problem - electric->solve(); - electric->exportResults(); - - return 0; -}