diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml new file mode 100644 index 000000000..9b71c1c01 --- /dev/null +++ b/.github/workflows/linux.yml @@ -0,0 +1,119 @@ +on: + #schedule: + # - cron: '30 6 * * *' #0630 UTC kick off full build and regression + #push: + # branches: [ ci/* ] + workflow_dispatch: + inputs: + runTests: + description: "Run tests?" + required: false + type: boolean + default: false + +env: + AWS_REGION: us-east-1 + AWS_SUBNET_ID: subnet-097642e5c610e7583 + AWS_SG_ID: sg-0ca7912782cf1538b + AMI_ID: ami-03a6eb9ec6101d323 + AWS_INSTANCE_TYPE: c5a.16xlarge + OS: linux + KEY_PAIR: win +jobs: + start-runner: + name: Start self-hosted EC2 runner + runs-on: ubuntu-latest + outputs: + label: ${{ steps.start-ec2-runner.outputs.label }} + ec2-instance-id: ${{ steps.start-ec2-runner.outputs.ec2-instance-id }} + steps: + - name: Configure AWS credentials + uses: aws-actions/configure-aws-credentials@v1 + with: + aws-access-key-id: ${{ secrets.AWS_ACCESS_KEY_ID }} + aws-secret-access-key: ${{ secrets.AWS_SECRET_ACCESS_KEY }} + aws-region: ${{ env.AWS_REGION }} + - name: Start EC2 runner + id: start-ec2-runner + uses: tundranerd/ec2-multiple-github-runners@multi-os-multi-runner + with: + mode: start + key-pair: ${{ env.KEY_PAIR }} + os: ${{ env.OS }} + github-token: ${{ secrets.GH_TOKEN}} + ec2-image-id: ${{ env.AMI_ID}} + ec2-instance-type: ${{ env.AWS_INSTANCE_TYPE }} + subnet-id: ${{ env.AWS_SUBNET_ID }} + security-group-id: ${{ env.AWS_SG_ID }} + wait-for-registry-timeout: 2 + aws-resource-tags: > # optional, requires additional permissions + [ + {"Key": "Name", "Value": "ec2-github-runner"}, + {"Key": "GitHubRepository", "Value": "${{ github.repository }}"}, + {"Key": "os", "Value": "${{ env.OS }}"} + ] + build: + name: Clone and build + needs: + - start-runner # required to get output from the start-runner job + runs-on: ${{ needs.start-runner.outputs.label }} # run the job on the newly created runner + steps: + - name: Checkout + uses: actions/checkout@v2 + - name: Build + run: | + ci/linux.sh + - name: Upload Artifacts + uses: actions/upload-artifact@v3 + with: + name: febio3-${{runner.os}}-${{runner.arch}} + path: cmbuild/bin + tests: + if: ${{inputs.runTests}} + name: Run test suite + needs: + - start-runner # required to get output from the start-runner job + - build + runs-on: ${{ needs.start-runner.outputs.label }} # run the job on the newly created runner + steps: + - name: Checkout + uses: actions/checkout@v2 + with: + submodules: 'true' + - uses: actions/download-artifact@v3 + with: + name: febio3-${{runner.os}}-${{runner.arch}} + path: cmbuild/bin + - name: Run test suite + run: | + ci/linux-test.sh + - name: Upload Artifacts + uses: actions/upload-artifact@v3 + with: + name: testsuite-${{runner.os}}-${{runner.arch}}-logs + path: | + TestSuite/Verify3/*.log + TestSuite/Logs*.txt + stop-runner: + name: Stop self-hosted EC2 runner + needs: + - start-runner # required to get output from the start-runner job + - tests # required to wait when the main job is done + runs-on: ubuntu-latest + if: ${{ always() }} # required to stop the runner even if the error happened in the previous jobs + steps: + - name: Configure AWS credentials + uses: aws-actions/configure-aws-credentials@v1 + with: + aws-access-key-id: ${{ secrets.AWS_ACCESS_KEY_ID }} + aws-secret-access-key: ${{ secrets.AWS_SECRET_ACCESS_KEY }} + aws-region: ${{ env.AWS_REGION }} + - name: Stop EC2 runner + uses: tundranerd/ec2-multiple-github-runners@multi-os-multi-runner + with: + mode: stop + os: ${{ env.OS }} + wait-for-deregistry-timeout: 5 + github-token: ${{ secrets.GH_TOKEN}} + label: ${{ needs.start-runner.outputs.label }} + ec2-instance-id: ${{ needs.start-runner.outputs.ec2-instance-id }} diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml new file mode 100644 index 000000000..486e53537 --- /dev/null +++ b/.github/workflows/windows.yml @@ -0,0 +1,118 @@ +on: + # schedule: + # - cron: '00 7 * * *' #0630 UTC kick off full build and regression + # push: + # branches: [ ci/* ] + workflow_dispatch: + inputs: + runTests: + description: "Run tests?" + required: false + type: boolean + default: false +env: + AWS_REGION: us-east-1 + AWS_SUBNET_ID: subnet-097642e5c610e7583 + AWS_SG_ID: sg-0ca7912782cf1538b + AMI_ID: ami-04bb42ae8c0b24120 + AWS_INSTANCE_TYPE: c5a.16xlarge + OS: windows + KEY_PAIR: win +jobs: + start-runner: + name: Start self-hosted EC2 runner + runs-on: ubuntu-latest + outputs: + label: ${{ steps.start-ec2-runner.outputs.label }} + ec2-instance-id: ${{ steps.start-ec2-runner.outputs.ec2-instance-id }} + steps: + - name: Configure AWS credentials + uses: aws-actions/configure-aws-credentials@v1 + with: + aws-access-key-id: ${{ secrets.AWS_ACCESS_KEY_ID }} + aws-secret-access-key: ${{ secrets.AWS_SECRET_ACCESS_KEY }} + aws-region: ${{ env.AWS_REGION }} + - name: Start EC2 runner + id: start-ec2-runner + uses: tundranerd/ec2-multiple-github-runners@multi-os-multi-runner + with: + mode: start + key-pair: ${{ env.KEY_PAIR }} + os: ${{ env.OS }} + github-token: ${{ secrets.GH_TOKEN}} + ec2-image-id: ${{ env.AMI_ID}} + ec2-instance-type: ${{ env.AWS_INSTANCE_TYPE }} + subnet-id: ${{ env.AWS_SUBNET_ID }} + security-group-id: ${{ env.AWS_SG_ID }} + wait-for-registry-timeout: 5 + aws-resource-tags: > # optional, requires additional permissions + [ + {"Key": "Name", "Value": "ec2-github-runner"}, + {"Key": "GitHubRepository", "Value": "${{ github.repository }}"}, + {"Key": "os", "Value": "${{ env.OS }}"} + ] + build: + name: Clone and build + needs: + - start-runner # required to get output from the start-runner job + runs-on: ${{ needs.start-runner.outputs.label }} # run the job on the newly created runner + steps: + - name: Checkout + uses: actions/checkout@v2 + - name: Build + run: | + ci/windows.bat + - name: Upload Artifacts + uses: actions/upload-artifact@v3 + with: + name: febio3-${{runner.os}}-${{runner.arch}} + path: cmbuild/bin + tests: + if: ${{inputs.runTests}} + name: Run test suite + needs: + - start-runner # required to get output from the start-runner job + - build + runs-on: ${{ needs.start-runner.outputs.label }} # run the job on the newly created runner + steps: + - name: Checkout + uses: actions/checkout@v2 + with: + submodules: 'true' + - uses: actions/download-artifact@v3 + with: + name: febio3-${{runner.os}}-${{runner.arch}} + path: cmbuild/bin + - name: Run test suite + run: | + ci/windows-test.bat + - name: Upload Artifacts + uses: actions/upload-artifact@v3 + with: + name: testsuite-${{runner.os}}-${{runner.arch}}-logs + path: | + TestSuite/Verify3/*.log + TestSuite/Logs*.txt + stop-runner: + name: Stop self-hosted EC2 runner + needs: + - start-runner # required to get output from the start-runner job + - tests # required to wait when the main job is done + runs-on: ubuntu-latest + if: ${{ always() }} # required to stop the runner even if the error happened in the previous jobs + steps: + - name: Configure AWS credentials + uses: aws-actions/configure-aws-credentials@v1 + with: + aws-access-key-id: ${{ secrets.AWS_ACCESS_KEY_ID }} + aws-secret-access-key: ${{ secrets.AWS_SECRET_ACCESS_KEY }} + aws-region: ${{ env.AWS_REGION }} + - name: Stop EC2 runner + uses: tundranerd/ec2-multiple-github-runners@multi-os-multi-runner + with: + mode: stop + os: ${{ env.OS }} + wait-for-deregistry-timeout: 5 + github-token: ${{ secrets.GH_TOKEN}} + label: ${{ needs.start-runner.outputs.label }} + ec2-instance-id: ${{ needs.start-runner.outputs.ec2-instance-id }} diff --git a/.gitignore b/.gitignore index 937c65d0a..4a3a02b49 100644 --- a/.gitignore +++ b/.gitignore @@ -33,9 +33,10 @@ bld/ build/ -# Visual Studio 2015/2017 cache/options directory +# Visual Studio cache/options directory .vs/ VS2019/ +VS2022/ # Uncomment if you have tasks that create the project's static files in wwwroot #wwwroot/ diff --git a/Documentation/FEBio3.bib b/Documentation/FEBio3.bib index 65b676045..54f4b8ecf 100644 --- a/Documentation/FEBio3.bib +++ b/Documentation/FEBio3.bib @@ -1,7 +1,7 @@ %% This BibTeX bibliography file was created using BibDesk. %% http://bibdesk.sourceforge.net/ -%% Created for Gerard Ateshian at 2022-04-16 18:28:43 -0400 +%% Created for Gerard Ateshian at 2022-07-20 20:08:16 -0400 %% Saved with string encoding Unicode (UTF-8) @@ -2127,11 +2127,12 @@ @article{BALZANI2012139 bdsk-url-1 = {https://www.sciencedirect.com/science/article/pii/S0045782511003616}, bdsk-url-2 = {https://doi.org/10.1016/j.cma.2011.11.015}} - @article{Birzle2019, - title={A coupled approach for identification of nonlinear and compressible material models for soft tissue based on different experimental setups – exemplified and detailed for Lung Parenchyma}, - volume={94}, - DOI={10.1016/j.jmbbm.2019.02.019}, - journal={Journal of the Mechanical Behavior of Biomedical Materials}, - author={Birzle, Anna M. and Martin, Christian and Uhlig, Stefan and Wall, Wolfgang A.}, - year={2019}, - pages={126–143}} +@article{Birzle2019, + author = {Birzle, Anna M. and Martin, Christian and Uhlig, Stefan and Wall, Wolfgang A.}, + date-modified = {2022-07-20 20:08:15 -0400}, + journal = {Journal of the Mechanical Behavior of Biomedical Materials}, + pages = {126--143}, + title = {A coupled approach for identification of nonlinear and compressible material models for soft tissue based on different experimental setups - exemplified and detailed for Lung Parenchyma}, + volume = {94}, + year = {2019}, + bdsk-url-1 = {https://doi.org/10.1016/j.jmbbm.2019.02.019}} diff --git a/Documentation/FEBio_Theory_Manual.lyx b/Documentation/FEBio_Theory_Manual.lyx index f18a4f3de..9d84ffac2 100644 --- a/Documentation/FEBio_Theory_Manual.lyx +++ b/Documentation/FEBio_Theory_Manual.lyx @@ -179,11 +179,11 @@ status open \series default -Theory Manual Version 3.7 +Theory Manual Version 3.8 \end_layout \begin_layout Date -Last Updated: June 14, 2022 +Last Updated: December 12, 2022 \end_layout \begin_layout Paragraph* @@ -7533,7 +7533,7 @@ When multiple solid species are present, the net solid mass content may , may be evaluated from \begin_inset Formula \begin{equation} -\varphi_{r}^{s}=\varphi_{0}^{s}+\sum\limits _{\sigma}\rho_{r}^{\sigma}/\rho_{T}^{\sigma},\label{eq156} +\varphi_{r}^{s}=\varphi_{0}^{s}+\sum\limits _{\sigma}\rho_{r}^{\sigma}/\rho_{T}^{\sigma},\label{eq:referential-solid-volume-fraction} \end{equation} \end_inset @@ -7551,10 +7551,10 @@ where \end_inset per volume of -\begin_inset Formula $\sigma)$ +\begin_inset Formula $\sigma$ \end_inset - and +) and \begin_inset Formula $\varphi_{0}^{s}$ \end_inset @@ -7591,48 +7591,175 @@ reference "eq154" \end_layout \begin_layout Standard -The various constituents of the solid matrix may be electrically charged. - Let +At the start of an analysis, the formula of eq. +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:referential-solid-volume-fraction" +plural "false" +caps "false" +noprefix "false" + +\end_inset + + takes the form +\begin_inset Formula +\begin{equation} +\varphi_{r}^{s}\left(0\right)=\varphi_{0}^{s}+\sum\limits _{\sigma}\rho_{r}^{\sigma}\left(0\right)/\rho_{T}^{\sigma}\label{eq:initial-referential-solid-fraction} +\end{equation} + +\end_inset + +where +\begin_inset Formula $\rho_{r}^{\sigma}\left(0\right)$ +\end_inset + + represents the initial apparent mass densities of solid constituents. + This initial value of +\begin_inset Formula $\varphi_{r}^{s}$ +\end_inset + + must satisfy +\begin_inset Formula $0\le\varphi_{r}^{s}\left(0\right)\le1$ +\end_inset + +, since +\begin_inset Formula $J=1$ +\end_inset + + at the start of an analysis. +\end_layout + +\begin_layout Standard +The various constituents of the porous solid matrix of a multiphasic mixture + may be electrically charged. + The charge density in the current configuration, normalized by the mixture + volume in the current configuration, is given by +\begin_inset Formula +\begin{equation} +\check{c}^{F}=\frac{z^{F}dn^{F}+\sum_{\sigma}z^{\sigma}dn^{\sigma}}{dV}\,,\label{eq:FCD-mixture-volume} +\end{equation} + +\end_inset + +where +\begin_inset Formula $z^{F}$ +\end_inset + + is the charge number (equivalent charge per mole) and +\begin_inset Formula $dn^{F}$ +\end_inset + + is the elemental number of moles associated with fixed charges present + in the initial volume fraction +\begin_inset Formula $\varphi_{0}^{s}$ +\end_inset + + of the solid matrix. + Similarly, \begin_inset Formula $z^{\sigma}$ \end_inset - be the charge number (equivalent charge per mole) of solid constituent - + is the charge number of evolving solid constituent \begin_inset Formula $\sigma$ \end_inset -, then the net referential fixed charge density of the solid matrix (equivalent - charge per fluid volume in the referential configuration) is given by +. + By convention however, the fixed charge density +\begin_inset Formula $c^{F}$ +\end_inset + + of a multiphasic material is normalized by the fluid volume of the mixture, + thus it can be shown that \begin_inset Formula \begin{equation} -c_{r}^{F}=\frac{1}{1-\varphi_{r}^{s}}\sum\limits _{\sigma}\frac{z^{\sigma}\rho_{r}^{\sigma}}{M^{\sigma}}\,,\label{eq157} +c^{F}=\frac{\check{c}^{F}}{1-\varphi^{s}}=\frac{1}{J-\varphi_{r}^{s}}\left(\frac{z^{F}dn^{F}}{dV_{r}}+\sum_{\sigma}\frac{z^{\sigma}\rho_{r}^{\sigma}}{M^{\sigma}}\right)\,.\label{eq:FCD-fluid-volume} \end{equation} \end_inset -where -\begin_inset Formula $M^{\sigma}$ +In particular, in the reference configuration this formula produces +\begin_inset Formula +\begin{equation} +c^{F}\left(0\right)=\frac{1}{1-\varphi_{r}^{s}\left(0\right)}\left(\frac{z^{F}dn^{F}\left(0\right)}{dV_{r}}+\sum_{\sigma}\frac{z^{\sigma}\rho_{r}^{\sigma}\left(0\right)}{M^{\sigma}}\right)\,.\label{eq:FCD-initial-time} +\end{equation} + \end_inset - is the molar mass of +We define the first term on the right-hand-side of this equation as the + initial value of the user-specified fixed charge density +\begin_inset Formula $c_{0}^{F}$ +\end_inset + + (which is not associated with solid species \begin_inset Formula $\sigma$ \end_inset - (an invariant quantity) and -\begin_inset Formula $1-\varphi_{r}^{s}$ +), +\begin_inset Formula +\begin{equation} +c_{0}^{F}\equiv\frac{1}{1-\varphi_{r}^{s}\left(0\right)}\frac{z^{F}dn^{F}\left(0\right)}{dV_{r}}\,.\label{eq:user-specified-cF0} +\end{equation} + +\end_inset + +Therefore, in the current configuration, we may rewrite +\begin_inset Formula +\begin{equation} +c^{F}=\frac{1-\varphi_{r}^{s}\left(0\right)}{J-\varphi_{r}^{s}}c_{0}^{F}+\frac{1}{J-\varphi_{r}^{s}}\sum_{\sigma}\frac{z^{\sigma}\rho_{r}^{\sigma}}{M^{\sigma}}\,.\label{eq158} +\end{equation} + +\end_inset + +where the user-specified +\begin_inset Formula $c_{0}^{F}$ \end_inset - represents the referential volume fraction of all fluid constituents (solvent + may optionally be associated with a load curve representing the scale factor -\begin_inset Formula $+$ +\begin_inset Formula $dn^{F}/dn^{F}\left(0\right)$ +\end_inset + +, in case the user would like to allow +\begin_inset Formula $c_{0}^{F}$ +\end_inset + + to evolve. + By analogy, we may now define the referential fixed-charge density +\begin_inset Formula $c_{r}^{F}$ +\end_inset + + of the mixture as +\begin_inset Formula +\begin{equation} +c_{r}^{F}\equiv\frac{J-\varphi_{r}^{s}}{1-\varphi_{r}^{s}\left(0\right)}c^{F}=c_{0}^{F}+\frac{1}{1-\varphi_{r}^{s}\left(0\right)}\sum_{\sigma}\frac{z^{\sigma}\rho_{r}^{\sigma}}{M^{\sigma}}\label{eq157} +\end{equation} + +\end_inset + +Here, +\begin_inset Formula $c_{r}^{F}$ \end_inset - solutes) in a saturated mixture. - Based on the kinematics of the continuum, the fixed charge density in the - current configuration is + may evolve with time, but it represents the number of fixed equivalent + charges in the current configuration, per fluid volume in the reference + configuration. + Hence we may also write +\begin_inset Formula +\[ +c^{F}=\frac{1-\varphi_{r}^{s}\left(0\right)}{J-\varphi_{r}^{s}}c_{r}^{F}\,. +\] + +\end_inset + + +\end_layout + +\begin_layout Standard +The molar concentration of a solid-bound molecular constituent, which may + be needed in a reactive process involving solutes, is given by \begin_inset Formula \begin{equation} -c^{F}=\frac{1-\varphi_{r}^{s}}{J-\varphi_{r}^{s}}c_{r}^{F}\,.\label{eq158} +c^{\sigma}=\frac{1}{J-\varphi_{r}^{s}}\frac{\rho_{r}^{\sigma}}{M^{\sigma}}\,.\label{eq:sbm-molar-concentration} \end{equation} \end_inset @@ -37789,7 +37916,7 @@ where is evaluated from \begin_inset CommandInset ref LatexCommand eqref -reference "eq156" +reference "eq:referential-solid-volume-fraction" \end_inset diff --git a/Documentation/FEBio_Theory_Manual.pdf b/Documentation/FEBio_Theory_Manual.pdf index 4d9055c81..fa6c979a0 100644 Binary files a/Documentation/FEBio_Theory_Manual.pdf and b/Documentation/FEBio_Theory_Manual.pdf differ diff --git a/Documentation/FEBio_User_Manual.lyx b/Documentation/FEBio_User_Manual.lyx index c2440170b..5b7e2264d 100644 --- a/Documentation/FEBio_User_Manual.lyx +++ b/Documentation/FEBio_User_Manual.lyx @@ -162,13 +162,13 @@ status open \end_inset -User's Manual Version 3.7 +User's Manual Version 3.8 \end_layout \begin_layout Date \series bold -Last Updated: June 14, 2022 +Last Updated: December 12, 2022 \end_layout \begin_layout Standard @@ -21763,7 +21763,7 @@ name "subsec:Fluid-Resistance" \begin_layout Standard In a fluid or fluid-FSI analysis, this boundary condition prescribes an elastic pressure -\begin_inset Formula $p=R\,Q+p_{0}$ +\begin_inset Formula $p=R\,q+p_{d}$ \end_inset on a surface, where @@ -21771,12 +21771,12 @@ In a fluid or fluid-FSI analysis, this boundary condition prescribes an \end_inset is the flow resistance and -\begin_inset Formula $Q$ +\begin_inset Formula $q$ \end_inset is the volumetric flow rate across the surface (calculated internally); -\begin_inset Formula $p_{0}$ +\begin_inset Formula $p_{d}$ \end_inset is an optional offset pressure. @@ -21804,6 +21804,164 @@ In a fluid or fluid-FSI analysis, this boundary condition prescribes an \end_layout +\begin_layout Subsubsection +Fluid RCR +\begin_inset CommandInset label +LatexCommand label +name "subsec:Fluid-RCR" + +\end_inset + + +\end_layout + +\begin_layout Standard +This surface load models a three-element Windkessel outflow condition on + a surface, as described in +\begin_inset CommandInset citation +LatexCommand citep +key "Vignon-Clementel06" +literal "false" + +\end_inset + +. + It prescribes an elastic pressure +\begin_inset Formula $p$ +\end_inset + + on a surface, which satisfies the ordinary differential equation +\begin_inset Formula +\[ +p+\tau\frac{dp}{dt}=\left(R_{d}+R\right)q+R\tau\frac{dq}{dt}+p_{d}+\tau\frac{dp_{d}}{dt},\quad\tau=R_{d}C +\] + +\end_inset + +where +\begin_inset Formula $R$ +\end_inset + + is the upstream flow resistance, +\begin_inset Formula $R_{d}$ +\end_inset + + is the downstream flow resistance, and +\begin_inset Formula $C$ +\end_inset + + is the downstream flow capacitance, whereas +\begin_inset Formula $q$ +\end_inset + + is the volumetric flow rate across the surface (calculated internally); + +\begin_inset Formula $p_{d}$ +\end_inset + + is an optional downstream offset pressure, which may be associated with + a loadcurve to produce a function of time. + The pressure +\begin_inset Formula $p$ +\end_inset + + obtained by solving this differential equation is converted to a dilatation + and prescribed as an essential boundary condition at the nodes of the facets. + Since the pressure is the solution of an ordinary differential equation, + the user may also specify the initial condition for the pressure. +\end_layout + +\begin_layout LyX-Code + +\end_layout + +\begin_layout LyX-Code + 1 +\end_layout + +\begin_layout LyX-Code + 10 +\end_layout + +\begin_layout LyX-Code + 0.5 +\end_layout + +\begin_layout LyX-Code + 0 +\end_layout + +\begin_layout LyX-Code + 0 +\end_layout + +\begin_layout LyX-Code + +\end_layout + +\begin_layout Standard +By setting +\begin_inset Formula $C=0$ +\end_inset + + we recover the fluid resistance surface load described in Section +\begin_inset space ~ +\end_inset + + +\begin_inset CommandInset ref +LatexCommand ref +reference "subsec:Fluid-Resistance" +plural "false" +caps "false" +noprefix "false" + +\end_inset + +, with flow resistance equal to +\begin_inset Formula $R_{d}+R$ +\end_inset + +. +\end_layout + +\begin_layout Standard +Internally, the ordinary differential equation presented above is solved + numerically using a standard finite difference scheme, +\begin_inset Formula +\[ +p_{n+1}=\left(\frac{R_{d}}{1+\frac{\tau}{\Delta t}}+R\right)q_{n+1}+\frac{\tau}{\Delta t+\tau}\left(p_{n}-p_{d,n}-Rq_{n}\right)+p_{d,n+1} +\] + +\end_inset + +where +\begin_inset Formula $f_{n+1}$ +\end_inset + + is the value of any function +\begin_inset Formula $f\left(t\right)$ +\end_inset + + at the current time step +\begin_inset Formula $t_{n+1}$ +\end_inset + +, +\begin_inset Formula $f_{n}$ +\end_inset + + is its value at the previous time step +\begin_inset Formula $t_{n}$ +\end_inset + +, and +\begin_inset Formula $\Delta t=t_{n+1}-t_{n}$ +\end_inset + + is the time increment. +\end_layout + \begin_layout Subsubsection Fluid-FSI Traction \begin_inset CommandInset label diff --git a/Documentation/FEBio_User_Manual.pdf b/Documentation/FEBio_User_Manual.pdf index c3274507a..11cfb16b2 100644 Binary files a/Documentation/FEBio_User_Manual.pdf and b/Documentation/FEBio_User_Manual.pdf differ diff --git a/Documentation/ReleaseNotes.txt b/Documentation/ReleaseNotes.txt index 0b919b4cc..020638df4 100644 --- a/Documentation/ReleaseNotes.txt +++ b/Documentation/ReleaseNotes.txt @@ -39,6 +39,51 @@ the FEBio User's Manual for a more detailed description of the new features. To report any bugs or request new features, please visit the FEBio forum at https://forums.febio.org/. +=========================================================================== + + R E L E A S E 3.8 12/12/2022 + +=========================================================================== + +* Several issues with the cold restart capability were addressed. +* An issue was caused that could cause a crash when rigid rotational dofs are coupled to a biphasic element. +* Added dump stride parameter, which allows users to specify how often the dump file is written. +* Added "contact_traction" log variable. +* Added names of parent when printing model parameter values. +* Suppor for the MKL-DSS solver was added, which seems to work better for certain types of models than the default pardiso solver. +* Added gap offset option to facet-2-facet sliding. +* Fixed issue with plot output logic when first step does not request output. +* Fixed bug in "element center of mass" plot variable. +* Parameters of "orthotropic elastic" material can now be mapped. +* Added "discrete element signed force" plot variable. +* Added log variables for discrete element force components. +* Fixed bug in compressible Gent material. +* Added "facet area" log variable. +* Added "deformation gradient" plot variable. +* Fixed issue with initializationg of FEParamMat3d parameters. +* Added FEMathExpression that refactored some of the logic previously in FEMathValue. FEMathValueVec3 now also uses FEMathExpression so that it can support maps. +* Added prescribed fiber active contraction models. (Similar to uniaxial, but fiber direction can be defined by user.) +* Added support for mixture stress log variable. +* The "mixture stress" plot variable now also works for uncoupled mixtures. +* Fixed calculation of fixed charge density in multiphasic and triphasic materials (in FEBio3) +* Added error exception in Holmes-Mow permeability to catch pathological case. +* Added "contact area" log variable. +* Updated Theory Manual to fix calculation of fixed charge density in reactive multiphasic materials that use charged solid-bound molecules +* Added description of fluid RCR surface load in User Manual +* Fixed the fluid RCR surface load and made it ready for release. Simplified the fluid resistance to be consistent with fluid RCR. +* Allow the initial referential mass density of solid-bound molecules, m_rho0, to be a FEParamDouble. +* Updated FEMembraneMassActionForward and FEMembraneMassActionReversible to use areal concentration of solid-bound molecules. Updated FEMembraneReactionRateIonChannel to accommodate solid-bound molecules in the calculation of the reaction rate. +* Added exponential damage CDF (useful for reactive viscoelastic materials with exponentially increasing weak bond fractions with increasing strain) +* Added aggressive interpolation type as user-specifiable option in BoomerAMG solver. +* Modified FEFiberPowLinear and FEFiberPowLinearUncoupled to use FEParamDouble material coefficients. +* Updated FEIdealGasIsentropic and FEIdealGasIsothermal to run properly in current version of FEBio. +* Point Source Updates +* Fixed bug in FEContinuousFiberDistribution::Stress +* Fixed issue with exporting plot parameter for mapped parameters of nodal loads. +* Added support for X, Y, Z in generic hyperelastic material. +* Restricted nodal reaction forces plot variable to non-rigid materials. + + =========================================================================== R E L E A S E 3.7 6/15/2022 diff --git a/FEAMR/FEMMGRemesh.cpp b/FEAMR/FEMMGRemesh.cpp index d64cc9895..0a82647b3 100644 --- a/FEAMR/FEMMGRemesh.cpp +++ b/FEAMR/FEMMGRemesh.cpp @@ -59,7 +59,6 @@ class FEMMGRemesh::MMG #endif BEGIN_FECORE_CLASS(FEMMGRemesh, FERefineMesh) - ADD_PARAMETER(m_maxiter, "max_iters"); ADD_PARAMETER(m_hmin, "min_element_size"); ADD_PARAMETER(m_hausd, "hausdorff"); ADD_PARAMETER(m_hgrad, "gradation"); @@ -72,7 +71,6 @@ END_FECORE_CLASS(); FEMMGRemesh::FEMMGRemesh(FEModel* fem) : FERefineMesh(fem) { - m_maxiter = 1; m_maxelem = 0; m_relativeSize = true; m_meshCoarsen = false; diff --git a/FEAMR/FEMMGRemesh.h b/FEAMR/FEMMGRemesh.h index d12ab86bf..2161e5732 100644 --- a/FEAMR/FEMMGRemesh.h +++ b/FEAMR/FEMMGRemesh.h @@ -44,7 +44,6 @@ class FEMMGRemesh : public FERefineMesh FEMeshAdaptorCriterion* GetCriterion() { return m_criterion; } private: - int m_maxiter; bool m_relativeSize; bool m_meshCoarsen; bool m_normalizeData; diff --git a/FEBioFluid/FEFluidRCRBC.cpp b/FEBioFluid/FEFluidRCRBC.cpp index bb9d4d5b8..29d78415f 100644 --- a/FEBioFluid/FEFluidRCRBC.cpp +++ b/FEBioFluid/FEFluidRCRBC.cpp @@ -29,15 +29,15 @@ SOFTWARE.*/ #include "FEFluid.h" #include "FEBioFluid.h" #include +#include //============================================================================= BEGIN_FECORE_CLASS(FEFluidRCRBC, FESurfaceLoad) -ADD_PARAMETER(m_R , "R"); -ADD_PARAMETER(m_Rd , "Rd"); -ADD_PARAMETER(m_p0, "initial_pressure"); -ADD_PARAMETER(m_pd, "pressure_offset"); -ADD_PARAMETER(m_C, "capacitance"); -ADD_PARAMETER(m_Bern, "Bernoulli"); + ADD_PARAMETER(m_R , "R"); + ADD_PARAMETER(m_Rd , "Rd"); + ADD_PARAMETER(m_p0, "initial_pressure"); + ADD_PARAMETER(m_pd, "pressure_offset"); + ADD_PARAMETER(m_C, "capacitance"); END_FECORE_CLASS(); //----------------------------------------------------------------------------- @@ -46,16 +46,10 @@ FEFluidRCRBC::FEFluidRCRBC(FEModel* pfem) : FESurfaceLoad(pfem), m_dofW(pfem) { m_R = 0.0; m_pfluid = nullptr; - m_alpha = 1.0; m_p0 = 0; m_Rd = 0.0; m_pd = 0.0; m_C = 0.0; - m_Bern = false; - - m_stepHist.clear(); - m_timeHist.clear(); - m_flowHist.clear(); m_dofW.AddVariable(FEBioFluid::GetVariableName(FEBioFluid::RELATIVE_FLUID_VELOCITY)); m_dofEF = pfem->GetDOFIndex(FEBioFluid::GetVariableName(FEBioFluid::FLUID_DILATATION), 0); @@ -83,9 +77,10 @@ bool FEFluidRCRBC::Init() m_pfluid = pm->ExtractProperty(); if (m_pfluid == nullptr) return false; - m_stepHist.resize(1,0.0); - m_timeHist.resize(1,0.0); - m_flowHist.resize(1,0.0); + m_pn = m_pp = m_p0; + m_pdn = m_pdp = m_pd; + m_qn = m_qp = 0; + m_tp = 0; return true; } @@ -108,64 +103,30 @@ void FEFluidRCRBC::Activate() //! Evaluate and prescribe the resistance pressure void FEFluidRCRBC::Update() { - // evaluate the flow rate - double Q = FlowRate(); - - int numsteps = GetFEModel()->GetCurrentStep()->m_ntimesteps; - FETimeInfo& tp = GetFEModel()->GetTime(); - - m_flowHist.resize(numsteps + 1, 0); - m_timeHist.resize(numsteps + 1); - m_stepHist.resize(numsteps + 1); - - m_timeHist[numsteps] = tp.currentTime; - m_stepHist[numsteps] = tp.timeIncrement; + // Check if we started a new time, if so, update variables + FETimeInfo& timeInfo = GetFEModel()->GetTime(); + double time = timeInfo.currentTime; + int iter = timeInfo.currentIteration; + double dt = timeInfo.timeIncrement; + if ((time > m_tp) && (iter == 0)) { + m_pp = m_pn; + m_qp = m_qn; + m_pdp = m_pdn; + m_tp = time; + } - m_flowHist[numsteps] = Q; + // evaluate the flow rate at the current time + m_qn = FlowRate(); + m_pdn = m_pd; double tau = m_Rd*m_C; - // calculate the resistance pressure - double pR = 0; - if (m_Bern) - pR = m_R*Q*abs(Q); - else - pR = m_R*Q; - - //calculate initial pressure contribution - double pi = 0.0; - if (tau > 0) - pi = m_p0*exp(-m_timeHist[numsteps]/tau); - - //calculate pressure from capacitor contribution - double pC = 0.0; - if (m_C > 0 && m_Rd > 0) - { - for (int i = 0; i<=numsteps; ++i) - { - double p1 = 0; - double p2 = 0; - if (m_Bern) - { - if (i != 0) - p1 = exp(-(m_timeHist[numsteps]-m_timeHist[i-1])/tau)/m_C*m_flowHist[i-1]*abs(m_flowHist[i-1]); - p2 = exp(-(m_timeHist[numsteps]-m_timeHist[i])/tau)/m_C*m_flowHist[i]*abs(m_flowHist[i]); - } - else - { - if (i != 0) - p1 = exp(-(m_timeHist[numsteps]-m_timeHist[i-1])/tau)/m_C*m_flowHist[i-1]; - p2 = exp(-(m_timeHist[numsteps]-m_timeHist[i])/tau)/m_C*m_flowHist[i]; - } - pC += (p1+p2)/2.0*m_stepHist[i]; - } - } - - double p = pR + pi + m_pd + pC; + // calculate the RCR pressure + m_pn = m_pdn + (m_Rd/(1+tau/dt)+m_R)*m_qn + tau/(dt+tau)*(m_pp - m_pdp - m_R*m_qp); // calculate the dilatation double e = 0; - bool good = m_pfluid->Dilatation(0,p,0,e); + bool good = m_pfluid->Dilatation(0,m_pn,0, e); assert(good); // prescribe this dilatation at the nodes @@ -185,7 +146,7 @@ void FEFluidRCRBC::Update() } //----------------------------------------------------------------------------- -//! evaluate the flow rate across this surface +//! evaluate the flow rate across this surface at current time double FEFluidRCRBC::FlowRate() { double Q = 0; @@ -206,8 +167,8 @@ double FEFluidRCRBC::FlowRate() // nodal coordinates for (int i=0; iGetMesh()->Node(el.m_node[i]); - rt[i] = node.m_rt*m_alpha + node.m_rp*(1-m_alpha); - vt[i] = node.get_vec3d(m_dofW[0], m_dofW[1], m_dofW[2])*m_alphaf + node.get_vec3d_prev(m_dofW[0], m_dofW[1], m_dofW[2])*(1-m_alphaf); + rt[i] = node.m_rt; + vt[i] = node.get_vec3d(m_dofW[0], m_dofW[1], m_dofW[2]); } double* Nr, *Ns; @@ -245,17 +206,6 @@ double FEFluidRCRBC::FlowRate() //! calculate residual void FEFluidRCRBC::LoadVector(FEGlobalVector& R, const FETimeInfo& tp) { - m_alpha = tp.alpha; m_alphaf = tp.alphaf; - - /* - int numsteps = GetFEModel()->GetCurrentStep()->m_ntimesteps; - m_flowHist.resize(numsteps + 1, 0); - m_timeHist.resize(numsteps + 1); - m_stepHist.resize(numsteps + 1); - - m_timeHist[numsteps] = tp.currentTime; - m_stepHist[numsteps] = tp.timeIncrement; - */ } //----------------------------------------------------------------------------- @@ -263,9 +213,5 @@ void FEFluidRCRBC::LoadVector(FEGlobalVector& R, const FETimeInfo& tp) void FEFluidRCRBC::Serialize(DumpStream& ar) { FESurfaceLoad::Serialize(ar); - ar & m_alpha & m_alphaf; - ar & m_pfluid; - ar & m_timeHist; - ar & m_flowHist; - ar & m_stepHist; + ar & m_pn & m_pp & m_qn & m_qp & m_pdn & m_pdp & m_tp; } diff --git a/FEBioFluid/FEFluidRCRBC.h b/FEBioFluid/FEFluidRCRBC.h index ee810840f..99b183308 100644 --- a/FEBioFluid/FEFluidRCRBC.h +++ b/FEBioFluid/FEFluidRCRBC.h @@ -29,8 +29,7 @@ SOFTWARE.*/ #include "FEFluidMaterial.h" //----------------------------------------------------------------------------- -//! FEFluidResistanceBC is a fluid surface that has a normal -//! pressure proportional to the flow rate (resistance). +//! FEFluidRCRBC is a fluid surface load that implements a 3-element Windkessel model //! class FEBIOFLUID_API FEFluidRCRBC : public FESurfaceLoad { @@ -65,15 +64,15 @@ class FEBIOFLUID_API FEFluidRCRBC : public FESurfaceLoad double m_p0; //!< initial fluid pressure double m_C; //!< capacitance double m_pd; //!< downstream pressure - bool m_Bern; //!< Use Bernoulli's Relation (Q*|Q|) - - vector m_stepHist; //!< history of time step size (recorded each step) - vector m_timeHist; //!< history of time at each step - vector m_flowHist; //!< history of flow rate at each step private: - double m_alpha; - double m_alphaf; + double m_pn; //!< fluid pressure at current time point + double m_pp; //!< fluid pressure at previous time point + double m_qn; //!< flow rate at current time point + double m_qp; //!< flow rate at previous time point + double m_pdn; //!< downstream fluid pressure at current time point + double m_pdp; //!< downstream fluid pressure at previous time point + double m_tp; //!< previous time FEFluidMaterial* m_pfluid; //!< pointer to fluid FEDofList m_dofW; diff --git a/FEBioFluid/FEFluidResistanceBC.cpp b/FEBioFluid/FEFluidResistanceBC.cpp index 54b3102db..aab824209 100644 --- a/FEBioFluid/FEFluidResistanceBC.cpp +++ b/FEBioFluid/FEFluidResistanceBC.cpp @@ -42,7 +42,6 @@ FEFluidResistanceBC::FEFluidResistanceBC(FEModel* pfem) : FESurfaceLoad(pfem), m { m_R = 0.0; m_pfluid = nullptr; - m_alpha = 1.0; m_p0 = 0; m_dofW.AddVariable(FEBioFluid::GetVariableName(FEBioFluid::RELATIVE_FLUID_VELOCITY)); @@ -71,6 +70,8 @@ bool FEFluidResistanceBC::Init() m_pfluid = pm->ExtractProperty(); if (m_pfluid == nullptr) return false; + + return true; } @@ -120,7 +121,7 @@ void FEFluidResistanceBC::Update() } //----------------------------------------------------------------------------- -//! evaluate the flow rate across this surface +//! evaluate the flow rate across this surface at the current time double FEFluidResistanceBC::FlowRate() { double Q = 0; @@ -141,8 +142,8 @@ double FEFluidResistanceBC::FlowRate() // nodal coordinates for (int i=0; iGetMesh()->Node(el.m_node[i]); - rt[i] = node.m_rt*m_alpha + node.m_rp*(1-m_alpha); - vt[i] = node.get_vec3d(m_dofW[0], m_dofW[1], m_dofW[2])*m_alphaf + node.get_vec3d_prev(m_dofW[0], m_dofW[1], m_dofW[2])*(1-m_alphaf); + rt[i] = node.m_rt; + vt[i] = node.get_vec3d(m_dofW[0], m_dofW[1], m_dofW[2]); } double* Nr, *Ns; @@ -175,19 +176,3 @@ double FEFluidResistanceBC::FlowRate() return Q; } - -//----------------------------------------------------------------------------- -//! calculate residual -void FEFluidResistanceBC::LoadVector(FEGlobalVector& R, const FETimeInfo& tp) -{ - m_alpha = tp.alpha; m_alphaf = tp.alphaf; -} - -//----------------------------------------------------------------------------- -//! serialization -void FEFluidResistanceBC::Serialize(DumpStream& ar) -{ - FESurfaceLoad::Serialize(ar); - ar & m_alpha & m_alphaf; - ar & m_pfluid; -} diff --git a/FEBioFluid/FEFluidResistanceBC.h b/FEBioFluid/FEFluidResistanceBC.h index a9f48593f..68d1c6b94 100644 --- a/FEBioFluid/FEFluidResistanceBC.h +++ b/FEBioFluid/FEFluidResistanceBC.h @@ -44,7 +44,7 @@ class FEBIOFLUID_API FEFluidResistanceBC : public FESurfaceLoad void StiffnessMatrix(FELinearSystem& LS, const FETimeInfo& tp) override {} //! calculate load vector - void LoadVector(FEGlobalVector& R, const FETimeInfo& tp) override; + void LoadVector(FEGlobalVector& R, const FETimeInfo& tp) override {} //! set the dilatation void Update() override; @@ -58,17 +58,11 @@ class FEBIOFLUID_API FEFluidResistanceBC : public FESurfaceLoad //! activate void Activate() override; - //! serialization - void Serialize(DumpStream& ar) override; - private: double m_R; //!< flow resistance double m_p0; //!< fluid pressure offset private: - double m_alpha; - double m_alphaf; - FEFluidMaterial* m_pfluid; //!< pointer to fluid FEDofList m_dofW; diff --git a/FEBioFluid/FEIdealGasIsentropic.cpp b/FEBioFluid/FEIdealGasIsentropic.cpp index e3454201a..1f7d61dc3 100644 --- a/FEBioFluid/FEIdealGasIsentropic.cpp +++ b/FEBioFluid/FEIdealGasIsentropic.cpp @@ -70,6 +70,14 @@ bool FEIdealGasIsentropic::Init() return true; } +//----------------------------------------------------------------------------- +//! elastic pressure +double FEIdealGasIsentropic::Pressure(FEMaterialPoint& mp) +{ + FEFluidMaterialPoint& fp = *mp.ExtractData(); + return Pressure(fp.m_ef,0); +} + //----------------------------------------------------------------------------- //! elastic pressure from dilatation double FEIdealGasIsentropic::Pressure(const double e, const double T) diff --git a/FEBioFluid/FEIdealGasIsentropic.h b/FEBioFluid/FEIdealGasIsentropic.h index d28f1f40a..7135b5080 100644 --- a/FEBioFluid/FEIdealGasIsentropic.h +++ b/FEBioFluid/FEIdealGasIsentropic.h @@ -42,6 +42,7 @@ class FEBIOFLUID_API FEIdealGasIsentropic : public FEFluid bool Init() override; //! elastic pressure + double Pressure(FEMaterialPoint& mp) override; double Pressure(const double e, const double T = 0) override; //! tangent of elastic pressure with respect to strain J diff --git a/FEBioFluid/FEIdealGasIsothermal.cpp b/FEBioFluid/FEIdealGasIsothermal.cpp index 74358084a..0be2c74c3 100644 --- a/FEBioFluid/FEIdealGasIsothermal.cpp +++ b/FEBioFluid/FEIdealGasIsothermal.cpp @@ -67,6 +67,14 @@ bool FEIdealGasIsothermal::Init() return true; } +//----------------------------------------------------------------------------- +//! elastic pressure +double FEIdealGasIsothermal::Pressure(FEMaterialPoint& mp) +{ + FEFluidMaterialPoint& fp = *mp.ExtractData(); + return Pressure(fp.m_ef,0); +} + //----------------------------------------------------------------------------- //! elastic pressure from dilatation double FEIdealGasIsothermal::Pressure(const double e, const double T) diff --git a/FEBioFluid/FEIdealGasIsothermal.h b/FEBioFluid/FEIdealGasIsothermal.h index b0bd37f4d..f672215f5 100644 --- a/FEBioFluid/FEIdealGasIsothermal.h +++ b/FEBioFluid/FEIdealGasIsothermal.h @@ -42,6 +42,7 @@ class FEBIOFLUID_API FEIdealGasIsothermal : public FEFluid bool Init() override; //! elastic pressure + double Pressure(FEMaterialPoint& mp) override; double Pressure(const double e, const double T = 0) override; //! tangent of elastic pressure with respect to strain J diff --git a/FEBioLib/FEBioModel.cpp b/FEBioLib/FEBioModel.cpp index a02072ebd..841d5883d 100644 --- a/FEBioLib/FEBioModel.cpp +++ b/FEBioLib/FEBioModel.cpp @@ -106,6 +106,7 @@ FEBioModel::FEBioModel() m_logLevel = 1; m_dumpLevel = FE_DUMP_NEVER; + m_dumpStride = 1; // --- I/O-Data --- m_ndebug = 0; @@ -162,6 +163,12 @@ void FEBioModel::SetDumpLevel(int dumpLevel) { m_dumpLevel = dumpLevel; } //! get the dump level int FEBioModel::GetDumpLevel() const { return m_dumpLevel; } +//! Set the dump stride +void FEBioModel::SetDumpStride(int n) { m_dumpStride = n; } + +//! get the dump stride +int FEBioModel::GetDumpStride() const { return m_dumpStride; } + //! Set the log level void FEBioModel::SetLogLevel(int logLevel) { m_logLevel = logLevel; } @@ -383,7 +390,7 @@ void FEBioModel::Write(unsigned int nevent) FEAnalysis* pstep = GetCurrentStep(); // update plot file - if (m_plot) WritePlot(nevent); + WritePlot(nevent); // Dump converged state to the archive DumpData(nevent); @@ -407,6 +414,9 @@ void FEBioModel::WritePlot(unsigned int nevent) // try to open the plot file if ((nevent == CB_INIT) || (nevent == CB_STEP_ACTIVE)) { + // If the first step did not request output, m_plot can still be null + if (m_plot == 0) InitPlotFile(); + if (m_plot->IsValid() == false) { // Add the plot objects @@ -613,13 +623,22 @@ void FEBioModel::DumpData(int nevent) // get the current step FEAnalysis* pstep = GetCurrentStep(); int ndump = GetDumpLevel(); + int stride = GetDumpStride(); if (ndump == FE_DUMP_NEVER) return; bool bdump = false; switch (nevent) { case CB_MAJOR_ITERS: - if (ndump == FE_DUMP_MAJOR_ITRS) bdump = true; + if (ndump == FE_DUMP_MAJOR_ITRS) + { + if (stride <= 1) bdump = true; + else + { + int niter = pstep->m_ntimesteps; + bdump = ((niter % stride) == 0); + } + } if ((ndump == FE_DUMP_MUST_POINTS) && (pstep->m_timeController) && (pstep->m_timeController->m_nmust >= 0)) bdump = true; break; case CB_STEP_SOLVED: if (ndump == FE_DUMP_STEP) bdump = true; break; @@ -1317,6 +1336,35 @@ void FEBioModel::SerializeDataStore(DumpStream& ar) // I N I T I A L I Z A T I O N //============================================================================= +//----------------------------------------------------------------------------- +// Initialize plot file +bool FEBioModel::InitPlotFile() +{ + FEBioPlotFile* pplt = new FEBioPlotFile(this); + m_plot = pplt; + + // set the software string + const char* szver = febio::getVersionString(); + char szbuf[256] = { 0 }; + sprintf(szbuf, "FEBio %s", szver); + pplt->SetSoftwareString(szbuf); + + // see if a valid plot file name is defined. + const std::string& splt = GetPlotFileName(); + if (splt.empty()) + { + // if not, we take the input file name and set the extension to .xplt + char sz[1024] = { 0 }; + strcpy(sz, GetInputFileName().c_str()); + char* ch = strrchr(sz, '.'); + if (ch) *ch = 0; + strcat(sz, ".xplt"); + SetPlotFilename(sz); + } + + return true; +} + //----------------------------------------------------------------------------- //! This function performs one-time-initialization stuff. All the different //! modules are initialized here as well. This routine also performs some @@ -1339,30 +1387,7 @@ bool FEBioModel::Init() FEAnalysis* step = GetCurrentStep(); if (step->GetPlotLevel() != FE_PLOT_NEVER) { - if (m_plot == 0) - { - pplt = new FEBioPlotFile(this); - m_plot = pplt; - - // set the software string - const char* szver = febio::getVersionString(); - char szbuf[256] = { 0 }; - sprintf(szbuf, "FEBio %s", szver); - pplt->SetSoftwareString(szbuf); - } - - // see if a valid plot file name is defined. - const std::string& splt = GetPlotFileName(); - if (splt.empty()) - { - // if not, we take the input file name and set the extension to .xplt - char sz[1024] = {0}; - strcpy(sz, GetInputFileName().c_str()); - char *ch = strrchr(sz, '.'); - if (ch) *ch = 0; - strcat(sz, ".xplt"); - SetPlotFilename(sz); - } + if (m_plot == 0) InitPlotFile(); } // initialize model data diff --git a/FEBioLib/FEBioModel.h b/FEBioLib/FEBioModel.h index 24032d21b..2795aeb7d 100644 --- a/FEBioLib/FEBioModel.h +++ b/FEBioLib/FEBioModel.h @@ -177,6 +177,12 @@ class FEBIOLIB_API FEBioModel : public FEMechModel //! get the dump level int GetDumpLevel() const; + //! Set the dump stride + void SetDumpStride(int n); + + //! get the dump stride + int GetDumpStride() const; + //! Set the log level void SetLogLevel(int logLevel); @@ -205,6 +211,7 @@ class FEBIOLIB_API FEBioModel : public FEMechModel int m_logLevel; //!< output level for log file int m_dumpLevel; //!< level or writing restart file + int m_dumpStride; //!< write dump file every nth iterations private: // accumulative statistics diff --git a/FEBioLib/version.h b/FEBioLib/version.h index 9a552bf30..7a0d5a066 100644 --- a/FEBioLib/version.h +++ b/FEBioLib/version.h @@ -35,7 +35,7 @@ SOFTWARE.*/ // SUBSUBVERSION is incremented when bugs are fixed. #define VERSION 3 -#define SUBVERSION 7 +#define SUBVERSION 8 #define SUBSUBVERSION 0 /////////////////////////////////////////////////////////////////////////////// diff --git a/FEBioMech/FEBioMechData.cpp b/FEBioMech/FEBioMechData.cpp index a9c7cab89..cbc66275b 100644 --- a/FEBioMech/FEBioMechData.cpp +++ b/FEBioMech/FEBioMechData.cpp @@ -231,6 +231,54 @@ double FELogContactPressure::value(FESurfaceElement& el) return Lm; } +//----------------------------------------------------------------------------- +double FELogContactTractionX::value(FESurfaceElement& el) +{ + FEContactSurface* ps = dynamic_cast(el.GetMeshPartition()); + if (ps == nullptr) return 0.0; + + FEContactInterface* pci = ps->GetContactInterface(); assert(pci); + if ((pci == 0) || pci->IsActive()) + { + vec3d tn; + ps->GetSurfaceTraction(el.m_lid, tn); + return tn.x; + } + return 0.0; +} + +//----------------------------------------------------------------------------- +double FELogContactTractionY::value(FESurfaceElement& el) +{ + FEContactSurface* ps = dynamic_cast(el.GetMeshPartition()); + if (ps == nullptr) return 0.0; + + FEContactInterface* pci = ps->GetContactInterface(); assert(pci); + if ((pci == 0) || pci->IsActive()) + { + vec3d tn; + ps->GetSurfaceTraction(el.m_lid, tn); + return tn.y; + } + return 0.0; +} + +//----------------------------------------------------------------------------- +double FELogContactTractionZ::value(FESurfaceElement& el) +{ + FEContactSurface* ps = dynamic_cast(el.GetMeshPartition()); + if (ps == nullptr) return 0.0; + + FEContactInterface* pci = ps->GetContactInterface(); assert(pci); + if ((pci == 0) || pci->IsActive()) + { + vec3d tn; + ps->GetSurfaceTraction(el.m_lid, tn); + return tn.z; + } + return 0.0; +} + //----------------------------------------------------------------------------- double FELogElemPosX::value(FEElement& el) { @@ -1805,6 +1853,63 @@ double FELogDiscreteElementForce::value(FEElement& el) return Fm; } +//----------------------------------------------------------------------------- +double FELogDiscreteElementForceX::value(FEElement& el) +{ + FEDiscreteElasticMaterialPoint* mp = dynamic_cast(el.GetMaterialPoint(0)); + if (mp) return mp->m_Ft.x; + else return 0.0; +} + +//----------------------------------------------------------------------------- +double FELogDiscreteElementForceY::value(FEElement& el) +{ + FEDiscreteElasticMaterialPoint* mp = dynamic_cast(el.GetMaterialPoint(0)); + if (mp) return mp->m_Ft.y; + else return 0.0; +} + +//----------------------------------------------------------------------------- +double FELogDiscreteElementForceZ::value(FEElement& el) +{ + FEDiscreteElasticMaterialPoint* mp = dynamic_cast(el.GetMaterialPoint(0)); + if (mp) return mp->m_Ft.z; + else return 0.0; +} + +//----------------------------------------------------------------------------- +double FELogElementMixtureStress::value(FEElement& el) +{ + if (m_comp < 0) return 0.0; + + double s = 0.0; + for (int n = 0; n < el.GaussPoints(); ++n) + { + FEMaterialPoint& mp = *el.GetMaterialPoint(n); + FEElasticMixtureMaterialPoint* mmp = mp.ExtractData< FEElasticMixtureMaterialPoint>(); + if (mmp) + { + if (m_comp < mmp->Components()) + { + FEElasticMaterialPoint& ep = *mmp->GetPointData(m_comp)->ExtractData(); + + switch (m_metric) + { + case 0: s += ep.m_s.xx(); break; + case 1: s += ep.m_s.xy(); break; + case 2: s += ep.m_s.yy(); break; + case 3: s += ep.m_s.xz(); break; + case 4: s += ep.m_s.yz(); break; + case 5: s += ep.m_s.zz(); break; + } + } + } + } + s /= (double)el.GaussPoints(); + + return s; +} + //----------------------------------------------------------------------------- double FELogRigidBodyR11::value(FERigidBody& rb) { return (rb.GetRotation().RotationMatrix()(0,0)); } double FELogRigidBodyR12::value(FERigidBody& rb) { return (rb.GetRotation().RotationMatrix()(0, 1)); } @@ -1957,3 +2062,20 @@ double FELogVolumePressure::value(FENLConstraint& rc) FEVolumeConstraint* prc = dynamic_cast(&rc); return (prc ? prc->m_s.m_p : 0); } + +//============================================================================= +double FELogContactArea::value(FESurface& surface) +{ + FEContactSurface* pcs = dynamic_cast(&surface); + if (pcs == 0) return 0.0; + + // make sure the corresponding contact interface is active + // (in case the parent was not set, we'll proceed regardless) + FEContactInterface* pci = pcs->GetContactInterface(); assert(pci); + if ((pci == 0) || pci->IsActive()) + { + double area = pcs->GetContactArea(); + return area; + } + return 0.0; +} diff --git a/FEBioMech/FEBioMechData.h b/FEBioMech/FEBioMechData.h index 6fe43f381..22d9a8e9d 100644 --- a/FEBioMech/FEBioMechData.h +++ b/FEBioMech/FEBioMechData.h @@ -33,6 +33,7 @@ SOFTWARE.*/ #include "ObjectDataRecord.h" #include #include +#include //============================================================================= // N O D E D A T A @@ -178,6 +179,30 @@ class FELogContactPressure : public FEFaceLogData double value(FESurfaceElement& el) override; }; +//----------------------------------------------------------------------------- +class FELogContactTractionX : public FEFaceLogData +{ +public: + FELogContactTractionX(FEModel* fem) : FEFaceLogData(fem) {} + double value(FESurfaceElement& el) override; +}; + +//----------------------------------------------------------------------------- +class FELogContactTractionY : public FEFaceLogData +{ +public: + FELogContactTractionY(FEModel* fem) : FEFaceLogData(fem) {} + double value(FESurfaceElement& el) override; +}; + +//----------------------------------------------------------------------------- +class FELogContactTractionZ : public FEFaceLogData +{ +public: + FELogContactTractionZ(FEModel* fem) : FEFaceLogData(fem) {} + double value(FESurfaceElement& el) override; +}; + //============================================================================= // E L E M E N T D A T A //============================================================================= @@ -995,6 +1020,51 @@ class FELogDiscreteElementForce : public FELogElemData double value(FEElement& el); }; +//----------------------------------------------------------------------------- +//! Discrete element force +class FELogDiscreteElementForceX : public FELogElemData +{ +public: + FELogDiscreteElementForceX(FEModel* fem) : FELogElemData(fem) {} + double value(FEElement& el); +}; + +//----------------------------------------------------------------------------- +//! Discrete element force +class FELogDiscreteElementForceY : public FELogElemData +{ +public: + FELogDiscreteElementForceY(FEModel* fem) : FELogElemData(fem) {} + double value(FEElement& el); +}; + +//----------------------------------------------------------------------------- +//! Discrete element force +class FELogDiscreteElementForceZ : public FELogElemData +{ +public: + FELogDiscreteElementForceZ(FEModel* fem) : FELogElemData(fem) {} + double value(FEElement& el); +}; + +//----------------------------------------------------------------------------- +class FELogElementMixtureStress : public FELogElemData +{ +public: + FELogElementMixtureStress(FEModel* fem, int n, int m) : FELogElemData(fem), m_comp(n), m_metric(m) {} + double value(FEElement& el); + +private: + int m_comp; + int m_metric; +}; + +template class FELogElementMixtureStress_T : public FELogElementMixtureStress +{ +public: + FELogElementMixtureStress_T(FEModel* fem) : FELogElementMixtureStress(fem, N, M) {} +}; + //============================================================================= // R I G I D B O D Y D A T A //============================================================================= @@ -1415,3 +1485,13 @@ class FELogVolumePressure : public FELogNLConstraintData double value(FENLConstraint& rc); }; +//============================================================================= +// S U R F A C E D A T A +//============================================================================= + +class FELogContactArea : public FELogSurfaceData +{ +public: + FELogContactArea(FEModel* fem) : FELogSurfaceData(fem) {} + double value(FESurface& surface) override; +}; diff --git a/FEBioMech/FEBioMechModule.cpp b/FEBioMech/FEBioMechModule.cpp index 291e2a6af..0739d2bd2 100644 --- a/FEBioMech/FEBioMechModule.cpp +++ b/FEBioMech/FEBioMechModule.cpp @@ -429,6 +429,8 @@ void FEBioMech::InitModule() REGISTER_FECORE_CLASS(FEPrescribedActiveContractionTransIsoUC, "uncoupled prescribed trans iso active contraction"); REGISTER_FECORE_CLASS(FEPrescribedActiveContractionIsotropic, "prescribed isotropic active contraction"); REGISTER_FECORE_CLASS(FEPrescribedActiveContractionIsotropicUC, "uncoupled prescribed isotropic active contraction"); + REGISTER_FECORE_CLASS(FEPrescribedActiveContractionFiber, "prescribed fiber active contraction"); + REGISTER_FECORE_CLASS(FEPrescribedActiveContractionFiberUC, "uncoupled prescribed fiber active contraction"); REGISTER_FECORE_CLASS(FEActiveFiberStress, "active fiber stress"); REGISTER_FECORE_CLASS(FEActiveFiberStressUC, "uncoupled active fiber stress"); @@ -470,6 +472,7 @@ void FEBioMech::InitModule() REGISTER_FECORE_CLASS(FEDamageCDFGamma, "CDF gamma"); REGISTER_FECORE_CLASS(FEDamageCDFUser, "CDF user"); REGISTER_FECORE_CLASS(FEDamageCDFPower, "CDF power"); + REGISTER_FECORE_CLASS(FEDamageCDFExp, "CDF exponential"); // damage criterion (used by damage and plastic materials) REGISTER_FECORE_CLASS(FEDamageCriterionSimo, "DC Simo"); @@ -665,6 +668,7 @@ void FEBioMech::InitModule() REGISTER_FECORE_CLASS(FEPlotSPRPrincStresses, "SPR principal stress"); REGISTER_FECORE_CLASS(FEPlotNodalStresses, "nodal stress"); REGISTER_FECORE_CLASS(FEPlotShellStrain, "shell strain"); // NOTE: Deprecated + REGISTER_FECORE_CLASS(FEPlotDeformationGradient, "deformation gradient"); REGISTER_FECORE_CLASS(FEPlotLagrangeStrain, "Lagrange strain"); REGISTER_FECORE_CLASS(FEPlotInfStrain, "infinitesimal strain"); REGISTER_FECORE_CLASS(FEPlotSPRLagrangeStrain, "SPR Lagrange strain"); @@ -722,6 +726,7 @@ void FEBioMech::InitModule() REGISTER_FECORE_CLASS(FEPlotDiscreteElementElongation, "discrete element elongation"); REGISTER_FECORE_CLASS(FEPlotDiscreteElementPercentElongation, "discrete element percent elongation"); REGISTER_FECORE_CLASS(FEPlotDiscreteElementForce, "discrete element force"); + REGISTER_FECORE_CLASS(FEPlotDiscreteElementSignedForce, "discrete element signed force"); REGISTER_FECORE_CLASS(FEPlotDiscreteElementStrainEnergy, "discrete element strain energy"); REGISTER_FECORE_CLASS(FEPlotContinuousDamage_D , "continuous damage"); REGISTER_FECORE_CLASS(FEPlotContinuousDamage_D1 , "continuous damage D1"); @@ -767,6 +772,9 @@ void FEBioMech::InitModule() // Derived from FELogFaceData REGISTER_FECORE_CLASS(FELogContactGap , "contact gap"); REGISTER_FECORE_CLASS(FELogContactPressure, "contact pressure"); + REGISTER_FECORE_CLASS(FELogContactTractionX, "contact_traction.x"); + REGISTER_FECORE_CLASS(FELogContactTractionY, "contact_traction.y"); + REGISTER_FECORE_CLASS(FELogContactTractionZ, "contact_traction.z"); //----------------------------------------------------------------------------- // Derived from FELogElemData @@ -896,6 +904,28 @@ void FEBioMech::InitModule() REGISTER_FECORE_CLASS(FELogDiscreteElementStretch , "discrete element stretch"); REGISTER_FECORE_CLASS(FELogDiscreteElementElongation, "discrete element elongation"); REGISTER_FECORE_CLASS(FELogDiscreteElementForce , "discrete element force" ); + REGISTER_FECORE_CLASS(FELogDiscreteElementForceX , "Fde.x"); + REGISTER_FECORE_CLASS(FELogDiscreteElementForceY , "Fde.y"); + REGISTER_FECORE_CLASS(FELogDiscreteElementForceZ , "Fde.z"); + REGISTER_FECORE_CLASS(FELogContactArea, "contact area"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 0, 0, "mixture_stress[0].xx"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 0, 1, "mixture_stress[0].xy"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 0, 2, "mixture_stress[0].yy"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 0, 3, "mixture_stress[0].xz"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 0, 4, "mixture_stress[0].yz"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 0, 5, "mixture_stress[0].zz"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 1, 0, "mixture_stress[1].xx"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 1, 1, "mixture_stress[1].xy"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 1, 2, "mixture_stress[1].yy"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 1, 3, "mixture_stress[1].xz"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 1, 4, "mixture_stress[1].yz"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 1, 5, "mixture_stress[1].zz"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 2, 0, "mixture_stress[2].xx"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 2, 1, "mixture_stress[2].xy"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 2, 2, "mixture_stress[2].yy"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 2, 3, "mixture_stress[2].xz"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 2, 4, "mixture_stress[2].yz"); + REGISTER_FECORE_CLASS_T2(FELogElementMixtureStress_T, 2, 5, "mixture_stress[2].zz"); //----------------------------------------------------------------------------- // Derived from FELogObjectData diff --git a/FEBioMech/FEBioMechPlot.cpp b/FEBioMech/FEBioMechPlot.cpp index 5d9f0dd0c..5deb9aa87 100644 --- a/FEBioMech/FEBioMechPlot.cpp +++ b/FEBioMech/FEBioMechPlot.cpp @@ -106,29 +106,33 @@ bool FEPlotNodeReactionForces::Save(FEMesh& m, FEDataStream& a) FEElasticSolidDomain* dom = dynamic_cast(&m.Domain(i)); if (dom) { - int NE = dom->Elements(); - for (int i = 0; i < NE; ++i) + FERigidMaterial* prm = dynamic_cast(dom->GetMaterial()); + if (prm == nullptr) { - // get the element - FESolidElement& el = dom->Element(i); - - if (el.isActive()) { - // element force vector - vector fe; - vector lm; - - // get the element force vector and initialize it to zero - int ndof = 3 * el.Nodes(); - fe.assign(ndof, 0); - - // calculate internal force vector - dom->ElementInternalForce(el, fe); - - // assemble into F - for (size_t j = 0; j < el.Nodes(); ++j) - { - vec3d rj(fe[3 * j], fe[3 * j + 1], fe[3 * j + 2]); - F[el.m_node[j]] += rj; + int NE = dom->Elements(); + for (int i = 0; i < NE; ++i) + { + // get the element + FESolidElement& el = dom->Element(i); + + if (el.isActive()) { + // element force vector + vector fe; + vector lm; + + // get the element force vector and initialize it to zero + int ndof = 3 * el.Nodes(); + fe.assign(ndof, 0); + + // calculate internal force vector + dom->ElementInternalForce(el, fe); + + // assemble into F + for (size_t j = 0; j < el.Nodes(); ++j) + { + vec3d rj(fe[3 * j], fe[3 * j + 1], fe[3 * j + 2]); + F[el.m_node[j]] += rj; + } } } } @@ -838,7 +842,8 @@ bool FEPlotElementMixtureStress::Save(FEDomain& dom, FEDataStream& a) // make sure this is a mixture FEElasticMixture* pmm = dynamic_cast(pmat); - if (pmm == nullptr) return false; + FEUncoupledElasticMixture* pum = dynamic_cast(pmat); + if ((pmm == nullptr) && (pum == nullptr)) return false; // get the mixture component if (m_comp < 0) return false; @@ -1233,10 +1238,11 @@ bool FEPlotElementCenterOfMass::Save(FEDomain &dom, FEDataStream& a) double m = 0; for (int j = 0; jExtractData()); + FEMaterialPoint& mp = *el.GetMaterialPoint(j); + FEElasticMaterialPoint& pt = *(mp.ExtractData()); double detJ = bd.detJ0(el, j)*gw[j]; - ew += pt.m_rt*(pme->Density(pt)*detJ); - m += pme->Density(pt)*detJ; + ew += pt.m_rt*(pme->Density(mp)*detJ); + m += pme->Density(mp)*detJ; } a << ew / m; @@ -2026,6 +2032,20 @@ bool FEPlotSPRPrincStresses::Save(FEDomain& dom, FEDataStream& a) return true; } +//----------------------------------------------------------------------------- +bool FEPlotDeformationGradient::Save(FEDomain& dom, FEDataStream& a) +{ + FEElasticMaterial* pme = dom.GetMaterial()->ExtractProperty(); + if (pme == nullptr) return false; + + writeAverageElementValue(dom, a, [](const FEMaterialPoint& mp) { + const FEElasticMaterialPoint& pt = *mp.ExtractData(); + return pt.m_F; + }); + + return true; +} + //============================================================================= //! Store the average Lagrangian strain class FELagrangeStrain @@ -3796,6 +3816,35 @@ bool FEPlotDiscreteElementForce::Save(FEDomain& dom, FEDataStream& a) return true; } +bool FEPlotDiscreteElementSignedForce::Save(FEDomain& dom, FEDataStream& a) +{ + FEDiscreteElasticDomain* pdiscreteDomain = dynamic_cast(&dom); + if (pdiscreteDomain == nullptr) return false; + FEDiscreteElasticDomain& discreteDomain = *pdiscreteDomain; + + int NE = discreteDomain.Elements(); + for (int i = 0; i < NE; ++i) + { + FEDiscreteElement& el = discreteDomain.Element(i); + + // get the (one) material point data + FEDiscreteElasticMaterialPoint& mp = dynamic_cast(*el.GetMaterialPoint(0)); + + vec3d ra1 = dom.Node(el.m_lnode[0]).m_rt; + vec3d rb1 = dom.Node(el.m_lnode[1]).m_rt; + vec3d e = rb1 - ra1; e.unit(); + + vec3d F = mp.m_Ft; + + double Fm = F * e; + + // write the force + a << Fm; + } + + return true; +} + bool FEPlotDiscreteElementStrainEnergy::Save(FEDomain& dom, FEDataStream& a) { FEDiscreteElasticDomain* pdiscreteDomain = dynamic_cast(&dom); diff --git a/FEBioMech/FEBioMechPlot.h b/FEBioMech/FEBioMechPlot.h index 604fb764f..c8f009333 100644 --- a/FEBioMech/FEBioMechPlot.h +++ b/FEBioMech/FEBioMechPlot.h @@ -819,6 +819,15 @@ class FEPlotNodalStresses : public FEPlotDomainData bool Save(FEDomain& dom, FEDataStream& a); }; +//----------------------------------------------------------------------------- +//! Deformation gradient +class FEPlotDeformationGradient : public FEPlotDomainData +{ +public: + FEPlotDeformationGradient(FEModel* pfem) : FEPlotDomainData(pfem, PLT_MAT3F, FMT_ITEM) {} + bool Save(FEDomain& dom, FEDataStream& a); +}; + //----------------------------------------------------------------------------- //! Lagrange strains class FEPlotLagrangeStrain : public FEPlotDomainData @@ -1010,6 +1019,15 @@ class FEPlotDiscreteElementForce : public FEPlotDomainData bool Save(FEDomain& dom, FEDataStream& a); }; +//----------------------------------------------------------------------------- +class FEPlotDiscreteElementSignedForce : public FEPlotDomainData +{ +public: + FEPlotDiscreteElementSignedForce(FEModel* fem) : FEPlotDomainData(fem, PLT_FLOAT, FMT_ITEM) {} + bool Save(FEDomain& dom, FEDataStream& a); +}; + + //----------------------------------------------------------------------------- class FEPlotDiscreteElementStrainEnergy : public FEPlotDomainData { diff --git a/FEBioMech/FEContinuousFiberDistribution.cpp b/FEBioMech/FEContinuousFiberDistribution.cpp index 36541da4f..f91200c42 100644 --- a/FEBioMech/FEContinuousFiberDistribution.cpp +++ b/FEBioMech/FEContinuousFiberDistribution.cpp @@ -104,7 +104,7 @@ mat3ds FEContinuousFiberDistribution::Stress(FEMaterialPoint& mp) // calculate the stress double wn = it->m_weight; - s += m_pFmat->FiberStress(pt, fp.FiberPreStretch(n0))*(R*wn); + s += m_pFmat->FiberStress(mp, fp.FiberPreStretch(n0))*(R*wn); } while (it->Next()); } diff --git a/FEBioMech/FEDamageCDF.cpp b/FEBioMech/FEDamageCDF.cpp index a01f49a97..d25745846 100644 --- a/FEBioMech/FEDamageCDF.cpp +++ b/FEBioMech/FEDamageCDF.cpp @@ -367,8 +367,8 @@ double FEDamageCDFUser::pdf(const double X) //----------------------------------------------------------------------------- // define the material parameters BEGIN_FECORE_CLASS(FEDamageCDFPower, FEDamageCDF) - ADD_PARAMETER(m_alpha, FE_RANGE_GREATER_OR_EQUAL(1.0), "alpha" ); - ADD_PARAMETER(m_mu0 , FE_RANGE_GREATER_OR_EQUAL(0.0), "mu0" ); + ADD_PARAMETER(m_alpha, "alpha" ); + ADD_PARAMETER(m_mu0 , "mu0" ); ADD_PARAMETER(m_mu1 , "mu1" ); ADD_PARAMETER(m_s , FE_RANGE_GREATER(0.0) , "scale" ); END_FECORE_CLASS(); @@ -418,6 +418,51 @@ double FEDamageCDFPower::pdf(const double X) return pdf; } +/////////////////////////////////////////////////////////////////////////////// +//----------------------------------------------------------------------------- +// define the material parameters +BEGIN_FECORE_CLASS(FEDamageCDFExp, FEDamageCDF) + ADD_PARAMETER(m_alpha, FE_RANGE_GREATER_OR_EQUAL(1.0), "alpha" ); + ADD_PARAMETER(m_mu0 , FE_RANGE_GREATER_OR_EQUAL(0.0), "mu0" ); + ADD_PARAMETER(m_mu1 , "mu1" ); + ADD_PARAMETER(m_s , FE_RANGE_GREATER(0.0) , "scale" ); +END_FECORE_CLASS(); + +//----------------------------------------------------------------------------- +//! Constructor. +FEDamageCDFExp::FEDamageCDFExp(FEModel* pfem) : FEDamageCDF(pfem) +{ + m_alpha = 2; + m_mu0 = 1; + m_mu1 = 0; + m_s = 1; +} + +//----------------------------------------------------------------------------- +// Exponential cumulative distribution function +double FEDamageCDFExp::cdf(const double X) +{ + double cdf = 0; + + // this CDF only admits positive values + if (X > 0) + cdf = m_mu0*exp(m_mu1*pow(X/m_s,m_alpha)); + + return cdf; +} + +// Exponential probability density function +double FEDamageCDFExp::pdf(const double X) +{ + double pdf = 0; + + // this CDF only admits positive values + if (X > 0) + pdf = m_mu0*m_mu1*m_alpha/m_s*pow(X/m_s,m_alpha-1)*exp(m_mu1*pow(X/m_s,m_alpha)); + + return pdf; +} + /////////////////////////////////////////////////////////////////////////////// //----------------------------------------------------------------------------- // define the material parameters diff --git a/FEBioMech/FEDamageCDF.h b/FEBioMech/FEDamageCDF.h index 0cf05327b..b804938a2 100644 --- a/FEBioMech/FEDamageCDF.h +++ b/FEBioMech/FEDamageCDF.h @@ -244,6 +244,31 @@ class FEDamageCDFPower : public FEDamageCDF DECLARE_FECORE_CLASS(); }; +//----------------------------------------------------------------------------- +// Exponential law cumulative distribution function + +class FEDamageCDFExp : public FEDamageCDF +{ +public: + FEDamageCDFExp(FEModel* pfem); + ~FEDamageCDFExp() {} + + //! cumulative distribution function + double cdf(const double X) override; + + //! probability density function + double pdf(const double X) override; + +public: + double m_alpha; //!< power exponent alpha + double m_mu0; //!< constant coeff + double m_mu1; //!< coeff of power + double m_s; //!< scale factor for argument + + // declare parameter list + DECLARE_FECORE_CLASS(); +}; + //----------------------------------------------------------------------------- // Quadratic polynomial cumulative distribution function diff --git a/FEBioMech/FEExplicitSolidSolver.cpp b/FEBioMech/FEExplicitSolidSolver.cpp index 6807cd150..6dc8419af 100644 --- a/FEBioMech/FEExplicitSolidSolver.cpp +++ b/FEBioMech/FEExplicitSolidSolver.cpp @@ -586,6 +586,18 @@ void FEExplicitSolidSolver::Serialize(DumpStream& ar) { FESolver::Serialize(ar); ar & m_nrhs & m_niter & m_nref & m_ntotref & m_naug & m_neq & m_nreq; + + if (ar.IsShallow()) return; + + ar & m_Ut; + ar & m_Mi; + + if (ar.IsLoading()) + { + m_Fn.assign(m_neq, 0); + m_Fr.assign(m_neq, 0); + m_ui.assign(m_neq, 0); + } } //----------------------------------------------------------------------------- diff --git a/FEBioMech/FEFacet2FacetTied.cpp b/FEBioMech/FEFacet2FacetTied.cpp index 02e7510d2..8ece496b5 100644 --- a/FEBioMech/FEFacet2FacetTied.cpp +++ b/FEBioMech/FEFacet2FacetTied.cpp @@ -42,12 +42,14 @@ BEGIN_FECORE_CLASS(FEFacet2FacetTied, FEContactInterface) ADD_PARAMETER(m_naugmin, "minaug" ); ADD_PARAMETER(m_naugmax, "maxaug" ); ADD_PARAMETER(m_stol , "search_tolerance"); + ADD_PARAMETER(m_gapoff , "gap_offset" ); END_FECORE_CLASS(); //----------------------------------------------------------------------------- FEFacetTiedSurface::Data::Data() { - m_vgap = vec3d(0,0,0); + m_vgap = vec3d(0,0,0); + m_vgap0 = vec3d(0,0,0); m_Lm = vec3d(0,0,0); m_rs = vec2d(0,0); m_pme = (FESurfaceElement*) 0; @@ -56,7 +58,7 @@ FEFacetTiedSurface::Data::Data() void FEFacetTiedSurface::Data::Serialize(DumpStream& ar) { FEContactMaterialPoint::Serialize(ar); - ar & m_vgap & m_Lm & m_rs; + ar & m_vgap & m_vgap0 & m_Lm & m_rs; } //----------------------------------------------------------------------------- @@ -97,6 +99,7 @@ void FEFacetTiedSurface::Serialize(DumpStream &ar) { Data& d = static_cast(*el.GetMaterialPoint(j)); ar << d.m_vgap; + ar << d.m_vgap0; ar << d.m_rs; ar << d.m_Lm; } @@ -112,6 +115,7 @@ void FEFacetTiedSurface::Serialize(DumpStream &ar) { Data& d = static_cast(*el.GetMaterialPoint(j)); ar >> d.m_vgap; + ar >> d.m_vgap0; ar >> d.m_rs; ar >> d.m_Lm; } @@ -130,6 +134,7 @@ void FEFacetTiedSurface::Serialize(DumpStream &ar) { Data& d = static_cast(*el.GetMaterialPoint(j)); ar << d.m_vgap; + ar << d.m_vgap0; ar << d.m_rs; ar << d.m_Lm; } @@ -145,6 +150,7 @@ void FEFacetTiedSurface::Serialize(DumpStream &ar) { Data& d = static_cast(*el.GetMaterialPoint(j)); ar >> d.m_vgap; + ar >> d.m_vgap0; ar >> d.m_rs; ar >> d.m_Lm; } @@ -174,6 +180,7 @@ FEFacet2FacetTied::FEFacet2FacetTied(FEModel* pfem) : FEContactInterface(pfem), m_naugmin = 0; m_naugmax = 10; m_stol = 0.0001; + m_gapoff = false; } //----------------------------------------------------------------------------- @@ -263,6 +270,28 @@ void FEFacet2FacetTied::Activate() // project primary surface onto secondary surface ProjectSurface(m_ss, m_ms); + + if (m_gapoff) + { + // loop over all primary elements + for (int i = 0; i < m_ss.Elements(); ++i) + { + // get the primary element + FESurfaceElement& se = m_ss.Element(i); + + // loop over all its integration points + int nint = se.GaussPoints(); + for (int j = 0; j < nint; ++j) + { + // get integration point data + FEFacetTiedSurface::Data& pt = static_cast(*se.GetMaterialPoint(j)); + + pt.m_vgap0 = pt.m_vgap; + pt.m_vgap = vec3d(0, 0, 0); + pt.m_gap = 0.0; + } + } + } } //----------------------------------------------------------------------------- @@ -310,6 +339,8 @@ void FEFacet2FacetTied::ProjectSurface(FEFacetTiedSurface& ss, FEFacetTiedSurfac // calculate gap pt.m_vgap = x - q; + + pt.m_gap = pt.m_vgap.norm(); } else pt.m_pme = 0; } @@ -365,7 +396,9 @@ void FEFacet2FacetTied::Update() vec3d q = me.eval(y, r, s); // calculate the gap function - pt.m_vgap = rn - q; + pt.m_vgap = (rn - q) - pt.m_vgap0; + + pt.m_gap = pt.m_vgap.norm(); } } } diff --git a/FEBioMech/FEFacet2FacetTied.h b/FEBioMech/FEFacet2FacetTied.h index ca165fa54..00106335b 100644 --- a/FEBioMech/FEFacet2FacetTied.h +++ b/FEBioMech/FEFacet2FacetTied.h @@ -44,7 +44,8 @@ class FEFacetTiedSurface : public FEContactSurface void Serialize(DumpStream& ar); public: - vec3d m_vgap; //!< gap function + vec3d m_vgap; //!< gap function + vec3d m_vgap0; //!< initial gap function vec3d m_Lm; //!< Lagrange multiplier vec2d m_rs; //!< natural coordinates on secondary surface element }; @@ -118,6 +119,7 @@ class FEFacet2FacetTied : public FEContactInterface double m_stol; //!< search tolerance int m_naugmax; //!< maximum nr of augmentations int m_naugmin; //!< minimum nr of augmentations + bool m_gapoff; //!< retain initial gap as offset DECLARE_FECORE_CLASS(); }; diff --git a/FEBioMech/FEFiberPowLinear.cpp b/FEBioMech/FEFiberPowLinear.cpp index 62ba0d455..359b2bc3b 100644 --- a/FEBioMech/FEFiberPowLinear.cpp +++ b/FEBioMech/FEFiberPowLinear.cpp @@ -53,9 +53,12 @@ mat3ds FEFiberPowLinear::FiberStress(FEMaterialPoint& mp, const vec3d& n0) FEElasticMaterialPoint& pt = *mp.ExtractData(); // initialize material constants - double I0 = m_lam0*m_lam0; - double ksi = m_E/4/(m_beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-m_beta); - double b = ksi*pow(I0-1, m_beta-1) + m_E/2/sqrt(I0); + double E = m_E(mp); + double lam0 = m_lam0(mp); + double beta = m_beta(mp); + double I0 = lam0*lam0; + double ksi = E/4/(beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-beta); + double b = ksi*pow(I0-1, beta-1) + E/2/sqrt(I0); // deformation gradient mat3d &F = pt.m_F; @@ -80,8 +83,8 @@ mat3ds FEFiberPowLinear::FiberStress(FEMaterialPoint& mp, const vec3d& n0) // calculate the fiber stress magnitude double sn = (In < I0) ? - 2*In*ksi*pow(In-1, m_beta-1) : - 2*b*In - m_E*sqrt(In); + 2*In*ksi*pow(In-1, beta-1) : + 2*b*In - E*sqrt(In); // calculate the fiber stress s = N*(sn/J); @@ -100,8 +103,11 @@ tens4ds FEFiberPowLinear::FiberTangent(FEMaterialPoint& mp, const vec3d& n0) FEElasticMaterialPoint& pt = *mp.ExtractData(); // initialize material constants - double I0 = m_lam0*m_lam0; - double ksi = m_E/4/(m_beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-m_beta); + double E = m_E(mp); + double lam0 = m_lam0(mp); + double beta = m_beta(mp); + double I0 = lam0*lam0; + double ksi = E/4/(beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-beta); // deformation gradient mat3d &F = pt.m_F; @@ -127,8 +133,8 @@ tens4ds FEFiberPowLinear::FiberTangent(FEMaterialPoint& mp, const vec3d& n0) // calculate modulus double cn = (In < I0) ? - 4*In*In*ksi*(m_beta-1)*pow(In-1, m_beta-2) : - m_E*sqrt(In); + 4*In*In*ksi*(beta-1)*pow(In-1, beta-2) : + E*sqrt(In); // calculate the fiber tangent c = NxN*(cn/J); @@ -149,9 +155,12 @@ double FEFiberPowLinear::FiberStrainEnergyDensity(FEMaterialPoint& mp, const vec FEElasticMaterialPoint& pt = *mp.ExtractData(); // initialize material constants - double I0 = m_lam0*m_lam0; - double ksi = m_E/4/(m_beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-m_beta); - double b = ksi*pow(I0-1, m_beta-1) + m_E/2/sqrt(I0); + double E = m_E(mp); + double lam0 = m_lam0(mp); + double beta = m_beta(mp); + double I0 = lam0*lam0; + double ksi = E/4/(beta-1)*pow(I0, -3./2.)*pow(I0-1, 2-beta); + double b = ksi*pow(I0-1, beta-1) + E/2/sqrt(I0); // loop over all integration points const double eps = 0; @@ -165,8 +174,8 @@ double FEFiberPowLinear::FiberStrainEnergyDensity(FEMaterialPoint& mp, const vec { // calculate strain energy density sed = (In < I0) ? - ksi/m_beta*pow(In-1, m_beta) : - b*(In-I0) - m_E*(sqrt(In)-sqrt(I0)) + ksi/m_beta*pow(I0-1, m_beta); + ksi/beta*pow(In-1, beta) : + b*(In-I0) - E*(sqrt(In)-sqrt(I0)) + ksi/beta*pow(I0-1, beta); } return sed; @@ -200,13 +209,6 @@ bool FEFiberExpPowLinear::Validate() { if (FEElasticFiberMaterial::Validate() == false) return false; - // initialize material constants - m_I0 = m_lam0*m_lam0; - m_ksi = m_E*pow(m_I0-1,2-m_beta)*exp(-m_alpha*pow(m_I0-1,m_beta)) - /(4*pow(m_I0,1.5)*(m_beta-1+m_alpha*m_beta*pow(m_I0-1,m_beta))); - m_b = m_E*(2*m_I0*(m_beta-0.5+m_alpha*m_beta*pow(m_I0-1,m_beta))-1) - /(4*pow(m_I0,1.5)*(m_beta-1+m_alpha*m_beta*pow(m_I0-1,m_beta))); - return true; } @@ -230,6 +232,17 @@ mat3ds FEFiberExpPowLinear::FiberStress(FEMaterialPoint& mp, const vec3d& n0) const double eps = 0; if (In - 1 >= eps) { + // initialize material constants + double E = m_E(mp); + double lam0 = m_lam0(mp); + double alpha = m_alpha(mp); + double beta = m_beta(mp); + double I0 = lam0*lam0; + double ksi = E*pow(I0-1,2-beta)*exp(-alpha*pow(I0-1,beta)) + /(4*pow(I0,1.5)*(beta-1+alpha*beta*pow(I0-1,beta))); + double b = E*(2*I0*(beta-0.5+alpha*beta*pow(I0-1,beta))-1) + /(4*pow(I0,1.5)*(beta-1+alpha*beta*pow(I0-1,beta))); + // get the global spatial fiber direction in current configuration vec3d nt = F*n0 / sqrt(In); @@ -237,9 +250,9 @@ mat3ds FEFiberExpPowLinear::FiberStress(FEMaterialPoint& mp, const vec3d& n0) mat3ds N = dyad(nt); // calculate the fiber stress magnitude - double sn = (In < m_I0) ? - 2 * In*exp(m_alpha*pow(In-1,m_beta))*m_ksi*pow(In - 1, m_beta - 1) : - 2 * m_b*In - m_E*sqrt(In); + double sn = (In < I0) ? + 2 * In*exp(alpha*pow(In-1,beta))*ksi*pow(In - 1, beta - 1) : + 2 * b*In - E*sqrt(In); // calculate the fiber stress s = N*(sn / J); @@ -271,6 +284,17 @@ tens4ds FEFiberExpPowLinear::FiberTangent(FEMaterialPoint& mp, const vec3d& n0) const double eps = m_epsf * std::numeric_limits::epsilon(); if (In >= 1 + eps) { + // initialize material constants + double E = m_E(mp); + double lam0 = m_lam0(mp); + double alpha = m_alpha(mp); + double beta = m_beta(mp); + double I0 = lam0*lam0; + double ksi = E*pow(I0-1,2-beta)*exp(-alpha*pow(I0-1,beta)) + /(4*pow(I0,1.5)*(beta-1+alpha*beta*pow(I0-1,beta))); + double b = E*(2*I0*(beta-0.5+alpha*beta*pow(I0-1,beta))-1) + /(4*pow(I0,1.5)*(beta-1+alpha*beta*pow(I0-1,beta))); + // get the global spatial fiber direction in current configuration vec3d nt = F*n0 / sqrt(In); @@ -279,9 +303,9 @@ tens4ds FEFiberExpPowLinear::FiberTangent(FEMaterialPoint& mp, const vec3d& n0) tens4ds NxN = dyad1s(N); // calculate modulus - double cn = (In < m_I0) ? - 4*m_ksi*In*In*exp(m_alpha*pow(In-1,m_beta))*pow(In-1,m_beta-2)*(m_beta-1+m_alpha*m_beta*pow(In-1,m_beta)) : - m_E*sqrt(In); + double cn = (In < I0) ? + 4*ksi*In*In*exp(alpha*pow(In-1,beta))*pow(In-1,beta-2)*(beta-1+alpha*beta*pow(In-1,beta)) : + E*sqrt(In); // calculate the fiber tangent c = NxN*(cn / J); @@ -312,16 +336,27 @@ double FEFiberExpPowLinear::FiberStrainEnergyDensity(FEMaterialPoint& mp, const const double eps = 0; if (In - 1 >= eps) { + // initialize material constants + double E = m_E(mp); + double lam0 = m_lam0(mp); + double alpha = m_alpha(mp); + double beta = m_beta(mp); + double I0 = lam0*lam0; + double ksi = E*pow(I0-1,2-beta)*exp(-alpha*pow(I0-1,beta)) + /(4*pow(I0,1.5)*(beta-1+alpha*beta*pow(I0-1,beta))); + double b = E*(2*I0*(beta-0.5+alpha*beta*pow(I0-1,beta))-1) + /(4*pow(I0,1.5)*(beta-1+alpha*beta*pow(I0-1,beta))); + // calculate strain energy density - if (m_alpha == 0) { - sed = (In < m_I0) ? - m_ksi / m_beta*pow(In - 1, m_beta) : - m_b*(In - m_I0) - m_E*(sqrt(In) - sqrt(m_I0)) + m_ksi / m_beta*pow(m_I0 - 1, m_beta); + if (alpha == 0) { + sed = (In < I0) ? + ksi / beta*pow(In - 1, beta) : + b*(In - I0) - E*(sqrt(In) - sqrt(I0)) + ksi / beta*pow(I0 - 1, beta); } else { - sed = (In < m_I0) ? - m_ksi / (m_alpha*m_beta)*(exp(m_alpha*pow(In - 1, m_beta))-1) : - m_b*(In - m_I0) - m_E*(sqrt(In) - sqrt(m_I0)) + m_ksi / (m_alpha*m_beta)*(exp(m_alpha*pow(m_I0 - 1, m_beta))-1); + sed = (In < I0) ? + ksi / (alpha*beta)*(exp(alpha*pow(In - 1, beta))-1) : + b*(In - I0) - E*(sqrt(In) - sqrt(I0)) + ksi / (alpha*beta)*(exp(alpha*pow(I0 - 1, beta))-1); } } diff --git a/FEBioMech/FEFiberPowLinear.h b/FEBioMech/FEFiberPowLinear.h index 4a36f8a28..87bc71084 100644 --- a/FEBioMech/FEFiberPowLinear.h +++ b/FEBioMech/FEFiberPowLinear.h @@ -51,9 +51,9 @@ class FEFiberPowLinear : public FEElasticFiberMaterial DECLARE_FECORE_CLASS(); public: - double m_E; // fiber modulus - double m_lam0; // stretch ratio at end of toe region - double m_beta; // power law exponent in toe region + FEParamDouble m_E; // fiber modulus + FEParamDouble m_lam0; // stretch ratio at end of toe region + FEParamDouble m_beta; // power law exponent in toe region double m_epsf; }; @@ -77,17 +77,12 @@ class FEFiberExpPowLinear : public FEElasticFiberMaterial double FiberStrainEnergyDensity(FEMaterialPoint& mp, const vec3d& a0) override; public: - double m_E; // fiber modulus - double m_lam0; // stretch ratio at end of toe region - double m_alpha; // exponential coefficient - double m_beta; // power law exponent in toe region + FEParamDouble m_E; // fiber modulus + FEParamDouble m_lam0; // stretch ratio at end of toe region + FEParamDouble m_alpha; // exponential coefficient + FEParamDouble m_beta; // power law exponent in toe region double m_epsf; -private: - double m_I0; // m_lam0^2 - double m_ksi; // power law coefficient in toe region - double m_b; // coefficient in linear region - // declare the parameter list DECLARE_FECORE_CLASS(); }; diff --git a/FEBioMech/FEFiberPowLinearUncoupled.cpp b/FEBioMech/FEFiberPowLinearUncoupled.cpp index 19803e95a..6d672a6ad 100644 --- a/FEBioMech/FEFiberPowLinearUncoupled.cpp +++ b/FEBioMech/FEFiberPowLinearUncoupled.cpp @@ -60,9 +60,11 @@ mat3ds FEFiberPowLinearUncoupled::DevFiberStress(FEMaterialPoint& mp, const vec3 // initialize material constants double E = m_E(mp); - double I0 = m_lam0*m_lam0; - double ksi = m_E(mp) / 4.0 / (m_beta - 1)*pow(I0, -3.0 / 2.0)*pow(I0 - 1.0, 2.0 - m_beta); - double b = ksi*pow(I0 - 1.0, m_beta - 1.0) + E / 2.0 / sqrt(I0); + double lam0 = m_lam0(mp); + double beta = m_beta(mp); + double I0 = lam0*lam0; + double ksi = E / 4.0 / (beta - 1)*pow(I0, -3.0 / 2.0)*pow(I0 - 1.0, 2.0 - beta); + double b = ksi*pow(I0 - 1.0, beta - 1.0) + E / 2.0 / sqrt(I0); // only take fibers in tension into consideration const double eps = 0; @@ -76,8 +78,8 @@ mat3ds FEFiberPowLinearUncoupled::DevFiberStress(FEMaterialPoint& mp, const vec3 // calculate the fiber stress magnitude double sn = (In < I0) ? - 2*In*ksi*pow(In-1, m_beta-1) : - 2*b*In - m_E(mp)*sqrt(In); + 2*In*ksi*pow(In-1, beta-1) : + 2*b*In - E*sqrt(In); // calculate the fiber stress s = N*(sn/J); @@ -120,20 +122,22 @@ tens4ds FEFiberPowLinearUncoupled::DevFiberTangent(FEMaterialPoint& mp, const ve // initialize material constants double E = m_E(mp); - double I0 = m_lam0 * m_lam0; - double ksi = m_E(mp) / 4.0 / (m_beta - 1) * pow(I0, -3.0 / 2.0) * pow(I0 - 1.0, 2.0 - m_beta); - double b = ksi * pow(I0 - 1.0, m_beta - 1.0) + E / 2.0 / sqrt(I0); + double lam0 = m_lam0(mp); + double beta = m_beta(mp); + double I0 = lam0*lam0; + double ksi = E / 4.0 / (beta - 1)*pow(I0, -3.0 / 2.0)*pow(I0 - 1.0, 2.0 - beta); + double b = ksi*pow(I0 - 1.0, beta - 1.0) + E / 2.0 / sqrt(I0); // calculate the fiber stress magnitude double sn = (In < I0) ? - 2*In*ksi*pow(In-1, m_beta-1) : + 2*In*ksi*pow(In-1, beta-1) : 2*b*In - E*sqrt(In); // calculate the fiber stress s = N*(sn/J); // calculate modulus - double cn = (In < I0) ? 4*In*In*ksi*(m_beta-1)*pow(In-1, m_beta-2) : E*sqrt(In); + double cn = (In < I0) ? 4*In*In*ksi*(beta-1)*pow(In-1, beta-2) : E*sqrt(In); // calculate the fiber tangent c = NxN*(cn/J); @@ -171,14 +175,16 @@ double FEFiberPowLinearUncoupled::DevFiberStrainEnergyDensity(FEMaterialPoint& m { // initialize material constants double E = m_E(mp); - double I0 = m_lam0 * m_lam0; - double ksi = m_E(mp) / 4.0 / (m_beta - 1) * pow(I0, -3.0 / 2.0) * pow(I0 - 1.0, 2.0 - m_beta); - double b = ksi * pow(I0 - 1.0, m_beta - 1.0) + E / 2.0 / sqrt(I0); + double lam0 = m_lam0(mp); + double beta = m_beta(mp); + double I0 = lam0*lam0; + double ksi = E / 4.0 / (beta - 1)*pow(I0, -3.0 / 2.0)*pow(I0 - 1.0, 2.0 - beta); + double b = ksi*pow(I0 - 1.0, beta - 1.0) + E / 2.0 / sqrt(I0); // calculate strain energy density sed = (In < I0) ? - ksi/m_beta*pow(In-1, m_beta) : - b*(In-I0) - E*(sqrt(In) - sqrt(I0)) + ksi/m_beta*pow(I0-1, m_beta); + ksi/beta*pow(In-1, beta) : + b*(In-I0) - E*(sqrt(In) - sqrt(I0)) + ksi/beta*pow(I0-1, beta); } return sed; diff --git a/FEBioMech/FEFiberPowLinearUncoupled.h b/FEBioMech/FEFiberPowLinearUncoupled.h index b425d2f79..5a698ac7b 100644 --- a/FEBioMech/FEFiberPowLinearUncoupled.h +++ b/FEBioMech/FEFiberPowLinearUncoupled.h @@ -40,18 +40,18 @@ class FEFiberPowLinearUncoupled : public FEElasticFiberMaterialUC //! Cauchy stress - virtual mat3ds DevFiberStress(FEMaterialPoint& mp, const vec3d& n0) override; + mat3ds DevFiberStress(FEMaterialPoint& mp, const vec3d& n0) override; // Spatial tangent - virtual tens4ds DevFiberTangent(FEMaterialPoint& mp, const vec3d& n0) override; + tens4ds DevFiberTangent(FEMaterialPoint& mp, const vec3d& n0) override; //! Strain energy density - virtual double DevFiberStrainEnergyDensity(FEMaterialPoint& mp, const vec3d& n0) override; + double DevFiberStrainEnergyDensity(FEMaterialPoint& mp, const vec3d& n0) override; public: - FEParamDouble m_E; // fiber modulus - double m_lam0; // stretch ratio at end of toe region - double m_beta; // power law exponent in toe region + FEParamDouble m_E; // fiber modulus + FEParamDouble m_lam0; // stretch ratio at end of toe region + FEParamDouble m_beta; // power law exponent in toe region // declare the parameter list DECLARE_FECORE_CLASS(); diff --git a/FEBioMech/FEGenericHyperelastic.cpp b/FEBioMech/FEGenericHyperelastic.cpp index c72124b05..720113596 100644 --- a/FEBioMech/FEGenericHyperelastic.cpp +++ b/FEBioMech/FEGenericHyperelastic.cpp @@ -45,7 +45,7 @@ bool FEGenericHyperelastic::Init() bool FEGenericHyperelastic::BuildMathExpressions() { - vector vars = { "I1", "I2", "J" }; + vector vars = { "I1", "I2", "J", "X", "Y", "Z"}; // add all user parameters FEParameterList& pl = GetParameterList(); @@ -122,7 +122,7 @@ mat3ds FEGenericHyperelastic::Stress(FEMaterialPoint& mp) double I1 = B.tr(); double I2 = 0.5*(I1*I1 - B2.tr()); - vector v = { I1, I2, J }; + vector v = { I1, I2, J, mp.m_r0.x, mp.m_r0.y, mp.m_r0.z }; for (int i = 0; i < m_param.size(); ++i) v.push_back(*m_param[i]); double W1 = m_W1.value_s(v); @@ -149,7 +149,7 @@ tens4ds FEGenericHyperelastic::Tangent(FEMaterialPoint& mp) double I1 = B.tr(); double I2 = 0.5*(I1*I1 - B2.tr()); - vector v = { I1, I2, J }; + vector v = { I1, I2, J, mp.m_r0.x, mp.m_r0.y, mp.m_r0.z }; for (int i = 0; i < m_param.size(); ++i) v.push_back(*m_param[i]); double W1 = m_W1.value_s(v); @@ -195,7 +195,7 @@ double FEGenericHyperelastic::StrainEnergyDensity(FEMaterialPoint& mp) double I1 = B.tr(); double I2 = 0.5*(I1*I1 - B2.tr()); - vector v = { I1, I2, J }; + vector v = { I1, I2, J, mp.m_r0.x, mp.m_r0.y, mp.m_r0.z }; for (int i = 0; i < m_param.size(); ++i) v.push_back(*m_param[i]); double W = m_W.value_s(v); diff --git a/FEBioMech/FEGenericHyperelasticUC.cpp b/FEBioMech/FEGenericHyperelasticUC.cpp index 1e5464fbe..ac560279a 100644 --- a/FEBioMech/FEGenericHyperelasticUC.cpp +++ b/FEBioMech/FEGenericHyperelasticUC.cpp @@ -48,7 +48,7 @@ bool FEGenericHyperelasticUC::Init() bool FEGenericHyperelasticUC::BuildMathExpressions() { - vector vars = { "I1", "I2" }; + vector vars = { "I1", "I2", "X", "Y", "Z" }; // add all user parameters FEParameterList& pl = GetParameterList(); @@ -125,7 +125,7 @@ mat3ds FEGenericHyperelasticUC::DevStress(FEMaterialPoint& mp) double I2 = 0.5*(I1*I1 - B2.tr()); // get strain energy derivatives - vector v = { I1, I2 }; + vector v = { I1, I2, mp.m_r0.x, mp.m_r0.y, mp.m_r0.z }; for (int i = 0; i < m_param.size(); ++i) v.push_back(*m_param[i]); double W1 = m_W1.value_s(v); double W2 = m_W2.value_s(v); @@ -158,7 +158,7 @@ tens4ds FEGenericHyperelasticUC::DevTangent(FEMaterialPoint& mp) double I2 = 0.5*(I1*I1 - B2.tr()); // get strain energy derivatives - vector v = { I1, I2 }; + vector v = { I1, I2, mp.m_r0.x, mp.m_r0.y, mp.m_r0.z }; for (int i = 0; i < m_param.size(); ++i) v.push_back(*m_param[i]); double W1 = m_W1.value_s(v); double W2 = m_W2.value_s(v); @@ -213,7 +213,7 @@ double FEGenericHyperelasticUC::DevStrainEnergyDensity(FEMaterialPoint& mp) double I2 = 0.5*(I1*I1 - B2.tr()); // evaluate (deviatoric) strain energy - vector v = { I1, I2 }; + vector v = { I1, I2, mp.m_r0.x, mp.m_r0.y, mp.m_r0.z }; for (int i = 0; i < m_param.size(); ++i) v.push_back(*m_param[i]); double W = m_W.value_s(v); diff --git a/FEBioMech/FEGentMaterial.cpp b/FEBioMech/FEGentMaterial.cpp index bf7882142..42626d610 100644 --- a/FEBioMech/FEGentMaterial.cpp +++ b/FEBioMech/FEGentMaterial.cpp @@ -122,7 +122,7 @@ mat3ds FECompressibleGentMaterial::Stress(FEMaterialPoint& mp) FEElasticMaterialPoint& ep = *mp.ExtractData(); double J = ep.m_J; - mat3ds b = ep.RightCauchyGreen(); + mat3ds b = ep.LeftCauchyGreen(); double I1 = b.tr(); double mu = m_G; @@ -146,7 +146,7 @@ tens4ds FECompressibleGentMaterial::Tangent(FEMaterialPoint& mp) FEElasticMaterialPoint& ep = *mp.ExtractData(); double J = ep.m_J; - mat3ds b = ep.RightCauchyGreen(); + mat3ds b = ep.LeftCauchyGreen(); double I1 = b.tr(); double mu = m_G; diff --git a/FEBioMech/FEMechModel.cpp b/FEBioMech/FEMechModel.cpp index f2dae4221..1882d3444 100644 --- a/FEBioMech/FEMechModel.cpp +++ b/FEBioMech/FEMechModel.cpp @@ -166,8 +166,8 @@ void FEMechModel::Reactivate() //----------------------------------------------------------------------------- bool FEMechModel::Reset() { - if (FEModel::Reset() == false) return false; - return m_prs->Reset(); + if (m_prs->Reset() == false) return false; + return FEModel::Reset(); } //----------------------------------------------------------------------------- diff --git a/FEBioMech/FEOrthoElastic.cpp b/FEBioMech/FEOrthoElastic.cpp index bcf1e7da8..40de8df44 100644 --- a/FEBioMech/FEOrthoElastic.cpp +++ b/FEBioMech/FEOrthoElastic.cpp @@ -33,26 +33,62 @@ SOFTWARE.*/ //----------------------------------------------------------------------------- // define the material parameters BEGIN_FECORE_CLASS(FEOrthoElastic, FEElasticMaterial) - ADD_PARAMETER(E1 , FE_RANGE_GREATER(0.0), "E1"); - ADD_PARAMETER(E2 , FE_RANGE_GREATER(0.0), "E2"); - ADD_PARAMETER(E3 , FE_RANGE_GREATER(0.0), "E3"); - ADD_PARAMETER(G12, FE_RANGE_GREATER_OR_EQUAL(0.0), "G12"); - ADD_PARAMETER(G23, FE_RANGE_GREATER_OR_EQUAL(0.0), "G23"); - ADD_PARAMETER(G31, FE_RANGE_GREATER_OR_EQUAL(0.0), "G31"); - ADD_PARAMETER(v12, "v12"); - ADD_PARAMETER(v23, "v23"); - ADD_PARAMETER(v31, "v31"); + ADD_PARAMETER(m_E1 , FE_RANGE_GREATER(0.0), "E1"); + ADD_PARAMETER(m_E2 , FE_RANGE_GREATER(0.0), "E2"); + ADD_PARAMETER(m_E3 , FE_RANGE_GREATER(0.0), "E3"); + ADD_PARAMETER(m_G12, FE_RANGE_GREATER_OR_EQUAL(0.0), "G12"); + ADD_PARAMETER(m_G23, FE_RANGE_GREATER_OR_EQUAL(0.0), "G23"); + ADD_PARAMETER(m_G31, FE_RANGE_GREATER_OR_EQUAL(0.0), "G31"); + ADD_PARAMETER(m_v12, "v12"); + ADD_PARAMETER(m_v23, "v23"); + ADD_PARAMETER(m_v31, "v31"); END_FECORE_CLASS(); +//----------------------------------------------------------------------------- +FEOrthoElastic::FEOrthoElastic(FEModel* pfem) : FEElasticMaterial(pfem) {} + //----------------------------------------------------------------------------- //! Material initialization. bool FEOrthoElastic::Validate() { if (FEElasticMaterial::Validate() == false) return false; - if (v12 > sqrt(E1/E2)) { feLogError("Invalid value for v12. Let v12 <= sqrt(E1/E2)"); return false; } - if (v23 > sqrt(E2/E3)) { feLogError("Invalid value for v23. Let v23 <= sqrt(E2/E3)"); return false; } - if (v31 > sqrt(E3/E1)) { feLogError("Invalid value for v31. Let v31 <= sqrt(E3/E1)"); return false; } + // do some sanity checks + if (m_v12.isConst() && m_E1.isConst() && m_E2.isConst()) + { + double v12 = m_v12.constValue(); + double E1 = m_E1.constValue(); + double E2 = m_E2.constValue(); + if (v12 > sqrt(E1 / E2)) { feLogError("Invalid value for v12. Let v12 <= sqrt(E1/E2)"); return false; } + } + + if (m_v23.isConst() && m_E2.isConst() && m_E3.isConst()) + { + double v23 = m_v23.constValue(); + double E2 = m_E2.constValue(); + double E3 = m_E3.constValue(); + + if (v23 > sqrt(E2 / E3)) { feLogError("Invalid value for v23. Let v23 <= sqrt(E2/E3)"); return false; } + } + + if (m_v31.isConst() && m_E1.isConst() && m_E3.isConst()) + { + double v31 = m_v31.constValue(); + double E1 = m_E1.constValue(); + double E3 = m_E3.constValue(); + + if (v31 > sqrt(E3 / E1)) { feLogError("Invalid value for v31. Let v31 <= sqrt(E3/E1)"); return false; } + } + + return true; +} + +//----------------------------------------------------------------------------- +bool FEOrthoElastic::EvaluateLameCoefficients(FEMaterialPoint& mp, double lam[3][3], double mu[3]) +{ + double E1 = m_E1 (mp); double E2 = m_E2 (mp); double E3 = m_E3 (mp); + double v12 = m_v12(mp); double v23 = m_v23(mp); double v31 = m_v31(mp); + double G12 = m_G12(mp); double G23 = m_G23(mp); double G31 = m_G31(mp); // Evaluate Lame coefficients mu[0] = G12 + G31 - G23; @@ -61,14 +97,15 @@ bool FEOrthoElastic::Validate() lam[0][0] = 1.0/E1; lam[0][1] = -v12/E1; lam[0][2] = -v31/E3; lam[1][0] = -v12/E1; lam[1][1] = 1.0/E2; lam[1][2] = -v23/E2; lam[2][0] = -v31/E3; lam[2][1] = -v23/E2; lam[2][2] = 1.0/E3; - + // check that compliance matrix is positive definite - mat3ds c(lam[0][0],lam[1][1],lam[2][2],lam[0][1],lam[1][2],lam[0][2]); + mat3ds c(lam[0][0], lam[1][1], lam[2][2], lam[0][1], lam[1][2], lam[0][2]); double l[3]; c.exact_eigen(l); if ((l[0] < 0) || (l[1] < 0) || (l[2] < 0)) { feLogError("Stiffness matrix is not positive definite."); + assert(false); return false; } @@ -118,6 +155,11 @@ mat3ds FEOrthoElastic::Stress(FEMaterialPoint& mp) a[i] = F*a0[i]/sqrt(K[i]); // Evaluate the texture direction in the current configuration A[i] = dyad(a[i]); // Evaluate the texture tensor in the current configuration } + + // evaluate the Lame coefficients + double lam[3][3] = { 0 }; + double mu[3] = { 0 }; + EvaluateLameCoefficients(mp, lam, mu); // Evaluate the stress mat3ds s; @@ -163,6 +205,11 @@ tens4ds FEOrthoElastic::Tangent(FEMaterialPoint& mp) A[i] = dyad(a[i]); // Evaluate the texture tensor in the current configuration } + // evaluate the Lame coefficients + double lam[3][3] = { 0 }; + double mu[3] = { 0 }; + EvaluateLameCoefficients(mp, lam, mu); + tens4ds C(0.0); for (i=0; i<3; i++) { C += mu[i]*K[i]*dyad4s(A[i],b); @@ -199,6 +246,12 @@ double FEOrthoElastic::StrainEnergyDensity(FEMaterialPoint& mp) AE2[i] = A0[i].dotdot(E2); } + // evaluate the Lame coefficients + double lam[3][3] = { 0 }; + double mu[3] = { 0 }; + EvaluateLameCoefficients(mp, lam, mu); + + // calculate strain energy double sed = mu[0]*AE2[0] + mu[1]*AE2[1] + mu[2]*AE2[2] +0.5*(lam[0][0]*AE[0]*AE[0]+lam[1][1]*AE[1]*AE[1]+lam[2][2]*AE[2]*AE[2]) +lam[0][1]*AE[0]*AE[1]+lam[1][2]*AE[1]*AE[2]+lam[2][0]*AE[2]*AE[0]; diff --git a/FEBioMech/FEOrthoElastic.h b/FEBioMech/FEOrthoElastic.h index 569b0d634..b4650780c 100644 --- a/FEBioMech/FEOrthoElastic.h +++ b/FEBioMech/FEOrthoElastic.h @@ -35,14 +35,12 @@ SOFTWARE.*/ class FEOrthoElastic : public FEElasticMaterial { public: - double E1, E2, E3; // Young's moduli - double v12, v23, v31; // Poisson's ratio - double G12, G23, G31; // Shear moduli - double lam[3][3]; // first Lame coefficients - double mu[3]; // second Lame coefficients + FEParamDouble m_E1, m_E2, m_E3; // Young's moduli + FEParamDouble m_v12, m_v23, m_v31; // Poisson's ratio + FEParamDouble m_G12, m_G23, m_G31; // Shear moduli public: - FEOrthoElastic(FEModel* pfem) : FEElasticMaterial(pfem) {} + FEOrthoElastic(FEModel* pfem); //! calculate stress at material point virtual mat3ds Stress(FEMaterialPoint& pt) override; @@ -56,6 +54,10 @@ class FEOrthoElastic : public FEElasticMaterial //! data initialization bool Validate() override; +private: + // Helper function for evaluating Lame coefficients at material point + bool EvaluateLameCoefficients(FEMaterialPoint& pt, double lam[3][3], double mu[3]); + // declare parameter list DECLARE_FECORE_CLASS(); }; diff --git a/FEBioMech/FEPrescribedActiveContractionUniaxial.cpp b/FEBioMech/FEPrescribedActiveContractionUniaxial.cpp index 262f43eb2..c36dc2227 100644 --- a/FEBioMech/FEPrescribedActiveContractionUniaxial.cpp +++ b/FEBioMech/FEPrescribedActiveContractionUniaxial.cpp @@ -96,3 +96,52 @@ tens4ds FEPrescribedActiveContractionUniaxial::Tangent(FEMaterialPoint &mp) return c; } + + +//============================================================================= +// define the material parameters +BEGIN_FECORE_CLASS(FEPrescribedActiveContractionFiber, FEElasticMaterial) + ADD_PARAMETER(m_T0 , "T0" ); + ADD_PARAMETER(m_n0 , "fiber" ); +END_FECORE_CLASS(); + +//----------------------------------------------------------------------------- +FEPrescribedActiveContractionFiber::FEPrescribedActiveContractionFiber(FEModel* pfem) : FEElasticMaterial(pfem) +{ + m_T0 = 0.0; + m_n0 = vec3d(1, 0, 0); +} + +//----------------------------------------------------------------------------- +mat3ds FEPrescribedActiveContractionFiber::Stress(FEMaterialPoint &mp) +{ + FEElasticMaterialPoint& pt = *mp.ExtractData(); + + // deformation gradient + mat3d &F = pt.m_F; + double J = pt.m_J; + + // get the local coordinate systems + mat3d Q = GetLocalCS(mp); + + // evaluate fiber direction in global coordinate system + vec3d n0 = Q*m_n0.unitVector(mp); + + // evaluate the deformed fiber direction + vec3d nt = F*n0; + mat3ds N = dyad(nt); + + // evaluate the active stress + double T0 = m_T0(mp); + mat3ds s = N*(T0/J); + + return s; +} + +//----------------------------------------------------------------------------- +tens4ds FEPrescribedActiveContractionFiber::Tangent(FEMaterialPoint &mp) +{ + tens4ds c; + c.zero(); + return c; +} diff --git a/FEBioMech/FEPrescribedActiveContractionUniaxial.h b/FEBioMech/FEPrescribedActiveContractionUniaxial.h index 75b3c5da6..7f656f30c 100644 --- a/FEBioMech/FEPrescribedActiveContractionUniaxial.h +++ b/FEBioMech/FEPrescribedActiveContractionUniaxial.h @@ -58,3 +58,25 @@ class FEPrescribedActiveContractionUniaxial : public FEElasticMaterial DECLARE_FECORE_CLASS(); }; + +//----------------------------------------------------------------------------- +// This material implements an active contraction model which can be used +// as a component of a solid mixture material. +class FEPrescribedActiveContractionFiber : public FEElasticMaterial +{ +public: + //! constructor + FEPrescribedActiveContractionFiber(FEModel* pfem); + + //! stress + mat3ds Stress(FEMaterialPoint& pt) override; + + //! tangent + tens4ds Tangent(FEMaterialPoint& pt) override; + +private: + FEParamDouble m_T0; // prescribed active stress + FEParamVec3 m_n0; // unit vector along fiber direction (local coordinate system) + + DECLARE_FECORE_CLASS(); +}; \ No newline at end of file diff --git a/FEBioMech/FEPrescribedActiveContractionUniaxialUC.cpp b/FEBioMech/FEPrescribedActiveContractionUniaxialUC.cpp index 9f0521aeb..d70ff1942 100644 --- a/FEBioMech/FEPrescribedActiveContractionUniaxialUC.cpp +++ b/FEBioMech/FEPrescribedActiveContractionUniaxialUC.cpp @@ -93,3 +93,50 @@ tens4ds FEPrescribedActiveContractionUniaxialUC::DevTangent(FEMaterialPoint &mp) { return tens4ds(0.0); } + + +//============================================================================= +// define the material parameters +BEGIN_FECORE_CLASS(FEPrescribedActiveContractionFiberUC, FEUncoupledMaterial) + ADD_PARAMETER(m_T0 , "T0" ); + ADD_PARAMETER(m_n0 , "fiber" ); +END_FECORE_CLASS(); + +//----------------------------------------------------------------------------- +FEPrescribedActiveContractionFiberUC::FEPrescribedActiveContractionFiberUC(FEModel* pfem) : FEUncoupledMaterial(pfem) +{ + m_T0 = 0.0; + m_n0 = vec3d(1, 0, 0); +} + +//----------------------------------------------------------------------------- +mat3ds FEPrescribedActiveContractionFiberUC::DevStress(FEMaterialPoint &mp) +{ + FEElasticMaterialPoint& pt = *mp.ExtractData(); + + // deformation gradient + double J = pt.m_J; + mat3d F = pt.m_F; + + // get the local coordinate systems + mat3d Q = GetLocalCS(mp); + + // evaluate fiber direction in global coordinate system + vec3d n0 = Q*m_n0.unitVector(mp); + + // evaluate the deformed fiber direction + vec3d nt = F*n0; + mat3ds N = dyad(nt); + + // evaluate the active stress + double T0 = m_T0(mp); + mat3ds s = N*(T0/J); + + return s; +} + +//----------------------------------------------------------------------------- +tens4ds FEPrescribedActiveContractionFiberUC::DevTangent(FEMaterialPoint &mp) +{ + return tens4ds(0.0); +} diff --git a/FEBioMech/FEPrescribedActiveContractionUniaxialUC.h b/FEBioMech/FEPrescribedActiveContractionUniaxialUC.h index 96710efb9..c0461a740 100644 --- a/FEBioMech/FEPrescribedActiveContractionUniaxialUC.h +++ b/FEBioMech/FEPrescribedActiveContractionUniaxialUC.h @@ -58,3 +58,25 @@ class FEPrescribedActiveContractionUniaxialUC : public FEUncoupledMaterial DECLARE_FECORE_CLASS(); }; + +//----------------------------------------------------------------------------- +// This material implements an active contraction model which can be used +// as a component of an uncoupled solid mixture material. +class FEPrescribedActiveContractionFiberUC : public FEUncoupledMaterial +{ +public: + //! constructor + FEPrescribedActiveContractionFiberUC(FEModel* pfem); + + //! stress + mat3ds DevStress(FEMaterialPoint& pt) override; + + //! tangent + tens4ds DevTangent(FEMaterialPoint& pt) override; + +private: + FEParamDouble m_T0; // prescribed active stress + FEParamVec3 m_n0; // unit vector along fiber direction (local coordinate system) + + DECLARE_FECORE_CLASS(); +}; diff --git a/FEBioMix/FEBioMix.cpp b/FEBioMix/FEBioMix.cpp index f22f4d316..c59efc36e 100644 --- a/FEBioMix/FEBioMix.cpp +++ b/FEBioMix/FEBioMix.cpp @@ -268,6 +268,7 @@ void FEBioMix::InitModule() REGISTER_FECORE_CLASS(FEPlotActualSoluteConcentration , "solute concentration"); REGISTER_FECORE_CLASS(FEPlotPartitionCoefficient , "partition coefficient"); REGISTER_FECORE_CLASS(FEPlotSoluteFlux , "solute flux" ); + REGISTER_FECORE_CLASS(FEPlotSoluteVolumetricFlux , "solute volumetric flux" ); REGISTER_FECORE_CLASS(FEPlotOsmolarity , "osmolarity"); REGISTER_FECORE_CLASS(FEPlotCurrentDensity , "current density" ); REGISTER_FECORE_CLASS(FEPlotFixedChargeDensity , "fixed charge density"); @@ -379,9 +380,10 @@ void FEBioMix::InitModule() //----------------------------------------------------------------------------- // classes derived from FEPlotData REGISTER_FECORE_CLASS(FEPlotReceptorLigandConcentration , "receptor-ligand concentration" ); - REGISTER_FECORE_CLASS(FEPlotSBMConcentration , "sbm concentration" ); + REGISTER_FECORE_CLASS(FEPlotSBMConcentration , "sbm concentration" ); + REGISTER_FECORE_CLASS(FEPlotSBMArealConcentration , "sbm areal concentration" ); REGISTER_FECORE_CLASS(FEPlotSBMRefAppDensity , "sbm referential apparent density"); - REGISTER_FECORE_CLASS(FEPlotOsmoticCoefficient , "osmotic coefficient" ); + REGISTER_FECORE_CLASS(FEPlotOsmoticCoefficient , "osmotic coefficient" ); //----------------------------------------------------------------------------- // Element log data diff --git a/FEBioMix/FEBioMixData.cpp b/FEBioMix/FEBioMixData.cpp index 7856890f6..6e54b0b82 100644 --- a/FEBioMix/FEBioMixData.cpp +++ b/FEBioMix/FEBioMixData.cpp @@ -301,7 +301,7 @@ double FELogElemPorosity::value(FEElement& el) const FEElasticMaterialPoint* et = (mp->ExtractData()); const FEBiphasicMaterialPoint* pt = (mp->ExtractData()); - double p = (et && pt ? (1 - pt->m_phi0 / et->m_J) : 0.0); + double p = (et && pt ? (1 - pt->m_phi0t / et->m_J) : 0.0); val += p; } return val / (double) nint; diff --git a/FEBioMix/FEBioMixPlot.cpp b/FEBioMix/FEBioMixPlot.cpp index db63d485a..9d2b51596 100644 --- a/FEBioMix/FEBioMixPlot.cpp +++ b/FEBioMix/FEBioMixPlot.cpp @@ -708,6 +708,81 @@ bool FEPlotSoluteFlux::Save(FEDomain &dom, FEDataStream& a) return true; } +//----------------------------------------------------------------------------- +FEPlotSoluteVolumetricFlux::FEPlotSoluteVolumetricFlux(FEModel* pfem) : FEPlotDomainData(pfem, PLT_ARRAY_VEC3F, FMT_ITEM) +{ + DOFS& dofs = pfem->GetDOFS(); + int nsol = dofs.GetVariableSize("concentration"); + SetArraySize(nsol); + + // collect the names + int ndata = pfem->GlobalDataItems(); + vector s; + for (int i = 0; i(pfem->GetGlobalData(i)); + if (ps) + { + s.push_back(ps->GetName()); + m_sol.push_back(ps->GetID()); + } + } + assert(nsol == (int)s.size()); + SetArrayNames(s); +} + +//----------------------------------------------------------------------------- +bool FEPlotSoluteVolumetricFlux::Save(FEDomain &dom, FEDataStream& a) +{ + FESoluteInterface* pm = dynamic_cast(dom.GetMaterial()); + if ((pm == 0) || (pm->Solutes() == 0)) return false; + + // figure out the local solute IDs. This depends on the material + int nsols = (int)m_sol.size(); + vector lid(nsols, -1); + int nsc = 0; + for (int i = 0; i<(int)m_sol.size(); ++i) + { + lid[i] = pm->FindLocalSoluteID(m_sol[i]); + if (lid[i] != -1) nsc++; + } + if (nsc == 0) return false; + + for (int i = 0; i()); + FEFluidSolutesMaterialPoint* spt = (mp.ExtractData()); + FESolutesMaterial::Point* spt2 = (mp.ExtractData()); + FEMultiphasicFSIMaterialPoint* mfpt = (mp.ExtractData()); + + if (pt && (pt->m_ca[nsid] > 0)) ew += pt->m_j[nsid]/pt->m_ca[nsid]; + else if (spt && (spt->m_ca[nsid] > 0)) ew += spt->m_j[nsid]/spt->m_ca[nsid]; + else if (spt2 && (spt2->m_ca[nsid] > 0)) ew += spt2->m_j[nsid]/spt2->m_ca[nsid]; + else if (mfpt && (mfpt->m_ca[nsid] > 0)) ew += mfpt->m_j[nsid]/mfpt->m_ca[nsid]; + } + + ew /= el.GaussPoints(); + + a << ew; + } + } + } + return true; +} + //----------------------------------------------------------------------------- bool FEPlotOsmolarity::Save(FEDomain &dom, FEDataStream& a) { @@ -850,6 +925,78 @@ bool FEPlotSBMConcentration::Save(FEDomain &dom, FEDataStream& a) return true; } +//================================================================================================= +// FEPlotSBMArealConcentration +//================================================================================================= + +//----------------------------------------------------------------------------- +FEPlotSBMArealConcentration::FEPlotSBMArealConcentration(FEModel* pfem) : FEPlotDomainData(pfem, PLT_ARRAY, FMT_ITEM) +{ + // count SBMs + int sbms = 0; + int ndata = pfem->GlobalDataItems(); + vector names; + for (int i=0; i(pfem->GetGlobalData(i)); + if (sbm) + { + names.push_back(sbm->GetName()); + m_sbm.push_back(sbm->GetID()); + sbms++; + } + } + + SetArraySize(sbms); + SetArrayNames(names); +} + +//----------------------------------------------------------------------------- +bool FEPlotSBMArealConcentration::Save(FEDomain &dom, FEDataStream& a) +{ + FEShellDomain* bsd = static_cast(&dom); + if (bsd == nullptr) return false; + + FEMultiphasic* pm = dynamic_cast (dom.GetMaterial()); + if (pm == 0) return false; + + // figure out the local SBM IDs. This depend on the material + int nsbm = (int)m_sbm.size(); + vector lid(nsbm, -1); + for (int i=0; i<(int)m_sbm.size(); ++i) + { + lid[i] = GetLocalSBMID(pm, m_sbm[i]); + } + + int N = dom.Elements(); + for (int i = 0; i()); + + if (pt) ew += pm->SBMArealConcentration(mp, nk); + } + ew /= el.GaussPoints(); + + a << ew; + } + } + } + return true; +} + //----------------------------------------------------------------------------- bool FEPlotElectricPotential::Save(FEDomain &dom, FEDataStream& a) { @@ -921,7 +1068,7 @@ bool FEPlotReferentialSolidVolumeFraction::Save(FEDomain &dom, FEDataStream& a) const FEBiphasicFSIMaterialPoint* bpt = (mp.ExtractData()); double phif0 = 0; if (pt) - phif0 = pt->m_phi0; + phif0 = pt->m_phi0t; else if (bpt) phif0 = bpt->m_phi0; return phif0; @@ -943,7 +1090,7 @@ bool FEPlotPorosity::Save(FEDomain &dom, FEDataStream& a) writeAverageElementValue(dom, a, [](const FEMaterialPoint& mp) { const FEElasticMaterialPoint* et = (mp.ExtractData()); const FEBiphasicMaterialPoint* pt = (mp.ExtractData()); - return (pt ? (1 - pt->m_phi0/et->m_J) : 0.0); + return (pt ? (1 - pt->m_phi0t/et->m_J) : 0.0); }); return true; } @@ -1032,7 +1179,7 @@ bool FEPlotReferentialFixedChargeDensity::Save(FEDomain &dom, FEDataStream& a) const FEBiphasicFSIMaterialPoint* bfspt = mp.ExtractData(); double cf = 0; if (spt) - cf = (ept->m_J - bpt->m_phi0)*spt->m_cF / (1 - bpt->m_phi0); + cf = (ept->m_J - bpt->m_phi0t)*spt->m_cF / (1 - bpt->m_phi0); else if (mfspt) cf = (ept->m_J - bfspt->m_phi0)*mfspt->m_cF / (1 - bfspt->m_phi0); return cf; diff --git a/FEBioMix/FEBioMixPlot.h b/FEBioMix/FEBioMixPlot.h index d98e12fb8..69b5ff6e6 100644 --- a/FEBioMix/FEBioMixPlot.h +++ b/FEBioMix/FEBioMixPlot.h @@ -142,6 +142,18 @@ class FEPlotSoluteFlux : public FEPlotDomainData vector m_sol; }; +//----------------------------------------------------------------------------- +//! Solute volumetric flux +class FEPlotSoluteVolumetricFlux : public FEPlotDomainData +{ +public: + FEPlotSoluteVolumetricFlux(FEModel* pfem); + bool Save(FEDomain& dom, FEDataStream& a); + +protected: + vector m_sol; +}; + //----------------------------------------------------------------------------- //! Osmolarity class FEPlotOsmolarity : public FEPlotDomainData @@ -161,6 +173,16 @@ class FEPlotSBMConcentration : public FEPlotDomainData vector m_sbm; }; +//----------------------------------------------------------------------------- +class FEPlotSBMArealConcentration : public FEPlotDomainData +{ +public: + FEPlotSBMArealConcentration(FEModel* pfem); + bool Save(FEDomain& dom, FEDataStream& a); +protected: + vector m_sbm; +}; + //----------------------------------------------------------------------------- //! Electric potential class FEPlotElectricPotential : public FEPlotDomainData @@ -348,4 +370,3 @@ class FEPlotFluidLoadSupport : public FEPlotSurfaceData FEPlotFluidLoadSupport(FEModel* pfem) : FEPlotSurfaceData(pfem, PLT_FLOAT, FMT_REGION){} bool Save(FESurface& surf, FEDataStream& a); }; - diff --git a/FEBioMix/FEBiphasic.cpp b/FEBioMix/FEBiphasic.cpp index 5a01266bc..f4b3f7a98 100644 --- a/FEBioMix/FEBiphasic.cpp +++ b/FEBioMix/FEBiphasic.cpp @@ -63,7 +63,7 @@ void FEBiphasicMaterialPoint::Serialize(DumpStream& ar) { FEMaterialPoint::Serialize(ar); ar & m_p & m_gradp & m_gradpp; - ar & m_w & m_pa & m_phi0 & m_phi0p & m_phi0hat & m_Jp; + ar & m_w & m_pa & m_phi0 & m_phi0t & m_phi0p & m_phi0hat & m_Jp; ar & m_ss; } @@ -73,7 +73,7 @@ void FEBiphasicMaterialPoint::Init() m_p = m_pa = 0; m_gradp = m_gradpp = vec3d(0,0,0); m_w = vec3d(0,0,0); - m_phi0 = m_phi0p = 0; + m_phi0 = m_phi0t = m_phi0p = 0; m_phi0hat = 0; m_Jp = 1; m_ss.zero(); @@ -145,7 +145,7 @@ double FEBiphasic::Porosity(FEMaterialPoint& pt) double J = et.m_J; // porosity // double phiw = 1 - m_phi0/J; - double phi0 = pet.m_phi0; + double phi0 = pet.m_phi0t; double phiw = 1 - phi0/J; // check for pore collapse // TODO: throw an error if pores collapse diff --git a/FEBioMix/FEBiphasic.h b/FEBioMix/FEBiphasic.h index 8a8a0de0a..5a73209fb 100644 --- a/FEBioMix/FEBiphasic.h +++ b/FEBioMix/FEBiphasic.h @@ -64,7 +64,8 @@ class FEBIOMIX_API FEBiphasicMaterialPoint : public FEMaterialPoint vec3d m_gradpp; //!< gradp at previous time vec3d m_w; //!< fluid flux double m_pa; //!< actual fluid pressure - double m_phi0; //!< referential solid volume fraction at current time + double m_phi0; //!< initial referential solid volume fraction + double m_phi0t; //!< referential solid volume fraction at current time double m_phi0p; //!< referential solid volume fraction at previous time double m_phi0hat; //!< referential solid volume fraction supply at current time double m_Jp; //!< determinant of solid deformation gradient at previous time diff --git a/FEBioMix/FEBiphasicShellDomain.cpp b/FEBioMix/FEBiphasicShellDomain.cpp index 7fe229b7a..f13965876 100644 --- a/FEBioMix/FEBiphasicShellDomain.cpp +++ b/FEBioMix/FEBiphasicShellDomain.cpp @@ -107,7 +107,7 @@ void FEBiphasicShellDomain::PreSolveUpdate(const FETimeInfo& timeInfo) pb.m_p = p; pb.m_gradp = gradient(el, pn, qn, j); pb.m_gradpp = pb.m_gradp; - pb.m_phi0p = pb.m_phi0; + pb.m_phi0p = pb.m_phi0t; mp.Update(timeInfo); } @@ -192,7 +192,7 @@ void FEBiphasicShellDomain::Reset() FEBiphasicMaterialPoint& pt = *(mp.ExtractData()); // initialize referential solid volume fraction - pt.m_phi0 = pmb->m_phi0(mp); + pt.m_phi0 = pt.m_phi0t = pmb->m_phi0(mp); }); } diff --git a/FEBioMix/FEBiphasicSolidDomain.cpp b/FEBioMix/FEBiphasicSolidDomain.cpp index 0f9d75f2a..351bd5142 100644 --- a/FEBioMix/FEBiphasicSolidDomain.cpp +++ b/FEBioMix/FEBiphasicSolidDomain.cpp @@ -121,7 +121,7 @@ void FEBiphasicSolidDomain::PreSolveUpdate(const FETimeInfo& timeInfo) pb.m_p = p; pb.m_gradp = gradient(el, pn, j); pb.m_gradpp = pb.m_gradp; - pb.m_phi0p = pb.m_phi0; + pb.m_phi0p = pb.m_phi0t; mp.Update(timeInfo); } @@ -198,7 +198,7 @@ void FEBiphasicSolidDomain::UnpackLM(FEElement& el, vector& lm) int neln_p = el.ShapeFunctions(degree_p); // allocate lm - lm.resize(neln_d*4); + lm.resize(neln_d*4 + 3*neln_d); // displacement dofs for (int i=0; i& lm) // now the pressure dofs lm[4*i + 3] = id[m_dofP]; + + // rigid rotational dofs + lm[4 * neln_d + 3 * i ] = id[m_dofR[0]]; + lm[4 * neln_d + 3 * i + 1] = id[m_dofR[1]]; + lm[4 * neln_d + 3 * i + 2] = id[m_dofR[2]]; } // substitute interface dofs for solid-shell interfaces @@ -246,7 +251,7 @@ void FEBiphasicSolidDomain::Reset() FEBiphasicMaterialPoint& pt = *(mp.ExtractData()); // initialize referential solid volume fraction - pt.m_phi0 = m_pMat->m_phi0(mp); + pt.m_phi0 = pt.m_phi0t = m_pMat->m_phi0(mp); }); } diff --git a/FEBioMix/FEBiphasicSolute.cpp b/FEBioMix/FEBiphasicSolute.cpp index 4c87e5afb..f5806b906 100644 --- a/FEBioMix/FEBiphasicSolute.cpp +++ b/FEBioMix/FEBiphasicSolute.cpp @@ -128,7 +128,7 @@ double FEBiphasicSolute::Porosity(FEMaterialPoint& pt) double J = et.m_J; // porosity // double phiw = 1 - m_phi0/J; - double phi0 = pet.m_phi0; + double phi0 = pet.m_phi0t; double phiw = 1 - phi0/J; // check for pore collapse // TODO: throw an error if pores collapse diff --git a/FEBioMix/FEBiphasicSoluteShellDomain.cpp b/FEBioMix/FEBiphasicSoluteShellDomain.cpp index dc0db6412..9366a12e2 100644 --- a/FEBioMix/FEBiphasicSoluteShellDomain.cpp +++ b/FEBioMix/FEBiphasicSoluteShellDomain.cpp @@ -152,7 +152,7 @@ void FEBiphasicSoluteShellDomain::InitMaterialPoints() pt.m_pa = m_pMat->Pressure(mp); // initialize referential solid volume fraction - pt.m_phi0 = m_pMat->m_phi0(mp); + pt.m_phi0t = m_pMat->m_phi0(mp); // calculate stress pm.m_s = m_pMat->Stress(mp); @@ -217,7 +217,7 @@ void FEBiphasicSoluteShellDomain::Reset() FESolutesMaterialPoint& ps = *(mp.ExtractData()); // initialize referential solid volume fraction - pt.m_phi0 = m_pMat->m_phi0(mp); + pt.m_phi0 = pt.m_phi0t = m_pMat->m_phi0(mp); // initialize multiphasic solutes ps.m_nsol = nsol; @@ -285,7 +285,7 @@ void FEBiphasicSoluteShellDomain::PreSolveUpdate(const FETimeInfo& timeInfo) pb.m_p = p; pb.m_gradp = gradient(el, pn, qn, j); - pb.m_phi0p = pb.m_phi0; + pb.m_phi0p = pb.m_phi0t; ps.m_c[0] = c; ps.m_gradc[0] = gradient(el, cn, dn, j); @@ -1325,7 +1325,7 @@ void FEBiphasicSoluteShellDomain::UpdateElementStress(int iel) ppt.m_phi0hat = 0; // ppt.m_phi0hat = pmb->GetSolid()->MolarMass()/pmb->GetSolid()->Density()*pmb->GetSolute()->m_pSupp->SolidSupply(mp); - ppt.m_phi0 = ppt.m_phi0p + ppt.m_phi0hat*dt; + ppt.m_phi0t = ppt.m_phi0p + ppt.m_phi0hat*dt; } } diff --git a/FEBioMix/FEBiphasicSoluteSolidDomain.cpp b/FEBioMix/FEBiphasicSoluteSolidDomain.cpp index d5ab6d1eb..a1aeaaca1 100644 --- a/FEBioMix/FEBiphasicSoluteSolidDomain.cpp +++ b/FEBioMix/FEBiphasicSoluteSolidDomain.cpp @@ -154,7 +154,7 @@ void FEBiphasicSoluteSolidDomain::InitMaterialPoints() pt.m_pa = m_pMat->Pressure(mp); // initialize referential solid volume fraction - pt.m_phi0 = m_pMat->m_phi0(mp); + pt.m_phi0t = m_pMat->m_phi0(mp); // calculate stress pm.m_s = m_pMat->Stress(mp); @@ -233,7 +233,7 @@ void FEBiphasicSoluteSolidDomain::Reset() FESolutesMaterialPoint& ps = *(mp.ExtractData()); // initialize referential solid volume fraction - pt.m_phi0 = m_pMat->m_phi0(mp); + pt.m_phi0 = pt.m_phi0t = m_pMat->m_phi0(mp); // initialize multiphasic solutes ps.m_nsol = nsol; @@ -304,7 +304,7 @@ void FEBiphasicSoluteSolidDomain::PreSolveUpdate(const FETimeInfo& timeInfo) pb.m_p = p; pb.m_gradp = gradient(el, pn, j); - pb.m_phi0p = pb.m_phi0; + pb.m_phi0p = pb.m_phi0t; ps.m_c[0] = c; ps.m_gradc[0] = gradient(el, ct, j); @@ -1146,7 +1146,7 @@ void FEBiphasicSoluteSolidDomain::UpdateElementStress(int iel) ppt.m_phi0hat = 0; // ppt.m_phi0hat = pmb->GetSolid()->MolarMass()/pmb->GetSolid()->Density()*pmb->GetSolute()->m_pSupp->SolidSupply(mp); - ppt.m_phi0 = ppt.m_phi0p + ppt.m_phi0hat*dt; + ppt.m_phi0t = ppt.m_phi0p + ppt.m_phi0hat*dt; } } diff --git a/FEBioMix/FEBiphasicSolver.cpp b/FEBioMix/FEBiphasicSolver.cpp index c0a1c2c6e..6fe069162 100644 --- a/FEBioMix/FEBiphasicSolver.cpp +++ b/FEBioMix/FEBiphasicSolver.cpp @@ -176,13 +176,13 @@ bool FEBiphasicSolver::Quasin() double normE1; // energy norm double normD; // displacement norm double normd; // displacement increment norm - double normRi; // initial residual norm - double normEi; // initial energy norm - double normEm; // max energy norm - double normDi; // initial displacement norm + double normRi = 0; // initial residual norm + double normEi = 0; // initial energy norm + double normEm = 0; // max energy norm + double normDi = 0; // initial displacement norm // poro convergence norms data - double normPi; // initial pressure norm + double normPi = 0; // initial pressure norm double normP; // current pressure norm double normp; // incremement pressure norm diff --git a/FEBioMix/FECarterHayes.cpp b/FEBioMix/FECarterHayes.cpp index a597e2943..f5b59b406 100644 --- a/FEBioMix/FECarterHayes.cpp +++ b/FEBioMix/FECarterHayes.cpp @@ -47,7 +47,7 @@ bool FECarterHayes::Init() if (FEElasticMaterial::Init() == false) return false; // get the parent material which must be a multiphasic material - FEMultiphasic* pMP = dynamic_cast (GetParent()); + FEMultiphasic* pMP = GetAncestor()->ExtractProperty(); if (pMP == 0) { feLogError("Parent material must be multiphasic"); return false; diff --git a/FEBioMix/FEChemicalReaction.cpp b/FEBioMix/FEChemicalReaction.cpp index d2e8c1348..6665bfa06 100644 --- a/FEBioMix/FEChemicalReaction.cpp +++ b/FEBioMix/FEChemicalReaction.cpp @@ -63,13 +63,13 @@ FEChemicalReaction::FEChemicalReaction(FEModel* pfem) : FEReaction(pfem) //----------------------------------------------------------------------------- bool FEChemicalReaction::Init() { - // initialize base class - FEReaction::Init(); - // set the parents for the reaction rates if (m_pFwd) m_pFwd->m_pReact = this; if (m_pRev) m_pRev->m_pReact = this; + // initialize base class + FEReaction::Init(); + // initialize the reaction coefficients int isol, isbm, itot; diff --git a/FEBioMix/FEConcentrationIndependentReaction.cpp b/FEBioMix/FEConcentrationIndependentReaction.cpp index 4e1699d1d..d62a4c7d8 100644 --- a/FEBioMix/FEConcentrationIndependentReaction.cpp +++ b/FEBioMix/FEConcentrationIndependentReaction.cpp @@ -66,7 +66,7 @@ mat3ds FEConcentrationIndependentReaction::Tangent_ReactionSupply_Strain(FEMater const int nsol = m_nsol; const int nsbm = (int)m_v.size() - nsol; double J = ept.m_J; - double phi0 = bpt.m_phi0; + double phi0 = bpt.m_phi0t; double kF = m_pFwd->ReactionRate(pt); mat3ds dkFde = m_pFwd->Tangent_ReactionRate_Strain(pt); diff --git a/FEBioMix/FEDiffAlbroIso.cpp b/FEBioMix/FEDiffAlbroIso.cpp index 7b960e360..4d2efa03d 100644 --- a/FEBioMix/FEDiffAlbroIso.cpp +++ b/FEBioMix/FEDiffAlbroIso.cpp @@ -157,7 +157,7 @@ mat3ds FEDiffAlbroIso::Diffusivity(FEMaterialPoint& mp) // solid volume fraction in reference configuration double phi0 = 0; if (ppt) - phi0 = ppt->m_phi0; + phi0 = ppt->m_phi0t; else if (bfpt) phi0 = bfpt->m_phi0; // porosity in current configuration @@ -200,7 +200,7 @@ tens4dmm FEDiffAlbroIso::Tangent_Diffusivity_Strain(FEMaterialPoint &mp) // solid volume fraction in reference configuration double phi0 = 0; if (ppt) - phi0 = ppt->m_phi0; + phi0 = ppt->m_phi0t; else if (bfpt) phi0 = bfpt->m_phi0; // porosity in current configuration @@ -256,7 +256,7 @@ mat3ds FEDiffAlbroIso::Tangent_Diffusivity_Concentration(FEMaterialPoint &mp, co // solid volume fraction in reference configuration double phi0 = 0; if (ppt) - phi0 = ppt->m_phi0; + phi0 = ppt->m_phi0t; else if (bfpt) phi0 = bfpt->m_phi0; // porosity in current configuration diff --git a/FEBioMix/FEDiffRefIso.cpp b/FEBioMix/FEDiffRefIso.cpp index 946f0868e..333d55b81 100644 --- a/FEBioMix/FEDiffRefIso.cpp +++ b/FEBioMix/FEDiffRefIso.cpp @@ -84,7 +84,7 @@ mat3ds FEDiffRefIso::Diffusivity(FEMaterialPoint& mp) // solid volume fraction in reference configuration double phi0 = 0; if (ppt) - phi0 = ppt->m_phi0; + phi0 = ppt->m_phi0t; else if (bfpt) phi0 = bfpt->m_phi0; @@ -119,7 +119,7 @@ tens4dmm FEDiffRefIso::Tangent_Diffusivity_Strain(FEMaterialPoint &mp) // solid volume fraction in reference configuration double phi0 = 0; if (ppt) - phi0 = ppt->m_phi0; + phi0 = ppt->m_phi0t; else if (bfpt) phi0 = bfpt->m_phi0; diff --git a/FEBioMix/FEMassActionForward.cpp b/FEBioMix/FEMassActionForward.cpp index bd5abc8f4..ea7dafd1c 100644 --- a/FEBioMix/FEMassActionForward.cpp +++ b/FEBioMix/FEMassActionForward.cpp @@ -141,7 +141,7 @@ mat3ds FEMassActionForward::Tangent_ReactionSupply_Strain(FEMaterialPoint& pt) FESolutesMaterialPoint& spt = *pt.ExtractData(); double J = ept.m_J; - double phi0 = bpt.m_phi0; + double phi0 = bpt.m_phi0t; double kF = m_pFwd->ReactionRate(pt); mat3ds dkFde = m_pFwd->Tangent_ReactionRate_Strain(pt); diff --git a/FEBioMix/FEMassActionForwardEffective.cpp b/FEBioMix/FEMassActionForwardEffective.cpp index 6ed507636..d4607913b 100644 --- a/FEBioMix/FEMassActionForwardEffective.cpp +++ b/FEBioMix/FEMassActionForwardEffective.cpp @@ -115,7 +115,7 @@ mat3ds FEMassActionForwardEffective::Tangent_ReactionSupply_Strain(FEMaterialPoi else { double J = ept.m_J; - double phi0 = bpt.m_phi0; + double phi0 = bpt.m_phi0t; double kF = m_pFwd->ReactionRate(pt); mat3ds dkFde = m_pFwd->Tangent_ReactionRate_Strain(pt); diff --git a/FEBioMix/FEMassActionReversible.cpp b/FEBioMix/FEMassActionReversible.cpp index ea1514e27..abf5892d4 100644 --- a/FEBioMix/FEMassActionReversible.cpp +++ b/FEBioMix/FEMassActionReversible.cpp @@ -219,7 +219,7 @@ mat3ds FEMassActionReversible::Tangent_ReactionSupply_Strain(FEMaterialPoint& pt else { J = ept.m_J; - phi0 = bpt.m_phi0; + phi0 = bpt.m_phi0t; } diff --git a/FEBioMix/FEMassActionReversibleEffective.cpp b/FEBioMix/FEMassActionReversibleEffective.cpp index b8936d471..0d4d01405 100644 --- a/FEBioMix/FEMassActionReversibleEffective.cpp +++ b/FEBioMix/FEMassActionReversibleEffective.cpp @@ -217,7 +217,7 @@ mat3ds FEMassActionReversibleEffective::Tangent_ReactionSupply_Strain(FEMaterial } else { - phi0 = bpt.m_phi0; + phi0 = bpt.m_phi0t; } // forward reaction diff --git a/FEBioMix/FEMembraneMassActionForward.cpp b/FEBioMix/FEMembraneMassActionForward.cpp index ec2f344d6..cd6e393f1 100644 --- a/FEBioMix/FEMembraneMassActionForward.cpp +++ b/FEBioMix/FEMembraneMassActionForward.cpp @@ -28,7 +28,7 @@ SOFTWARE.*/ #include "stdafx.h" #include "FEMembraneMassActionForward.h" -#include "FECore/DOFS.h" +#include //----------------------------------------------------------------------------- //! molar supply at material point @@ -57,7 +57,7 @@ double FEMembraneMassActionForward::ReactionSupply(FEMaterialPoint& pt) for (int i=0; i 0) { - double c = m_pMP->SBMConcentration(pt, i); + double c = m_pMP->SBMArealConcentration(pt, i); zhat *= pow(c, vR); } } diff --git a/FEBioMix/FEMembraneMassActionReversible.cpp b/FEBioMix/FEMembraneMassActionReversible.cpp index 25e0a3771..6724a6710 100644 --- a/FEBioMix/FEMembraneMassActionReversible.cpp +++ b/FEBioMix/FEMembraneMassActionReversible.cpp @@ -55,7 +55,7 @@ double FEMembraneMassActionReversible::FwdReactionSupply(FEMaterialPoint& pt) for (int i=0; i 0) { - double c = m_pMP->SBMConcentration(pt, i); + double c = m_pMP->SBMArealConcentration(pt, i); zhat *= pow(c, vR); } } @@ -110,7 +110,7 @@ double FEMembraneMassActionReversible::RevReactionSupply(FEMaterialPoint& pt) for (int i=0; i 0) { - double c = m_pMP->SBMConcentration(pt, i); + double c = m_pMP->SBMArealConcentration(pt, i); zhat *= pow(c, vP); } } diff --git a/FEBioMix/FEMembraneReactionRateIonChannel.cpp b/FEBioMix/FEMembraneReactionRateIonChannel.cpp index 9613b0d1a..bf259bb25 100644 --- a/FEBioMix/FEMembraneReactionRateIonChannel.cpp +++ b/FEBioMix/FEMembraneReactionRateIonChannel.cpp @@ -33,12 +33,14 @@ SOFTWARE.*/ BEGIN_FECORE_CLASS(FEMembraneReactionRateIonChannel, FEMembraneReactionRate) ADD_PARAMETER(m_g, FE_RANGE_GREATER_OR_EQUAL(0.0), "g"); ADD_PARAMETER(m_sol, "sol"); + ADD_PARAMETER(m_sbm, "sbm"); END_FECORE_CLASS(); bool FEMembraneReactionRateIonChannel::Init() { - // reset m_sol to be zero-based + // reset m_sol and m_sbm to be zero-based m_sol -= 1; + m_sbm -= 1; // membrane reaction rate is child of membrane reaction FEMembraneReaction* m_MRp = dynamic_cast(GetParent()); @@ -55,14 +57,28 @@ double FEMembraneReactionRateIonChannel::ReactionRate(FEMaterialPoint& pt) double R = GetFEModel()->GetGlobalConstant("R"); double T = GetFEModel()->GetGlobalConstant("T"); double Fc = GetFEModel()->GetGlobalConstant("Fc"); + FEMultiphasic* pMP = m_pReact->m_pMP; + assert(pMP); + double ksi = pMP->SBMArealConcentration(pt, m_sbm); double k = 0; if ((ci > 0) && (ce > 0)) - k = (ci != ce) ? R*T*m_g/pow(Fc*m_z,2)*log(ci/ce)/(ci-ce) : R*T*m_g/pow(Fc*m_z,2)/ce; + k = (ci != ce) ? R*T*m_g/ksi/pow(Fc*m_z,2)*log(ci/ce)/(ci-ce) : R*T*m_g/pow(Fc*m_z,2)/ce; return k; } +//! tangent of reaction rate with strain at material point +double FEMembraneReactionRateIonChannel::Tangent_ReactionRate_Strain(FEMaterialPoint& pt) +{ + // get the areal strain + FEElasticMaterialPoint& pe = *(pt.ExtractData()); + FEShellElement*sel = dynamic_cast(pt.m_elem); + assert(sel); + double Jg = pe.m_J*sel->Evaluate(sel->m_h0, pt.m_index)/sel->Evaluate(sel->m_ht, pt.m_index); + return ReactionRate(pt)/Jg; +} + double FEMembraneReactionRateIonChannel::Tangent_ReactionRate_Ci(FEMaterialPoint& pt, const int isol) { if (isol != m_sol) return 0; @@ -73,10 +89,13 @@ double FEMembraneReactionRateIonChannel::Tangent_ReactionRate_Ci(FEMaterialPoint double R = GetFEModel()->GetGlobalConstant("R"); double T = GetFEModel()->GetGlobalConstant("T"); double Fc = GetFEModel()->GetGlobalConstant("Fc"); - + FEMultiphasic* pMP = m_pReact->m_pMP; + assert(pMP); + double ksi = pMP->SBMArealConcentration(pt, m_sbm); + double dkdc = 0; if ((ci > 0) && (ce > 0)) - dkdc = (ci != ce) ? R*T/pow(Fc*m_z,2)*m_g*(ci*(1-log(ci/ce))-ce)/pow(ci-ce,2)/ci : -R*T*m_g/pow(ci*Fc*m_z,2)/2; + dkdc = (ci != ce) ? R*T/pow(Fc*m_z,2)*m_g/ksi*(ci*(1-log(ci/ce))-ce)/pow(ci-ce,2)/ci : -R*T*m_g/pow(ci*Fc*m_z,2)/2; return dkdc; } @@ -91,10 +110,13 @@ double FEMembraneReactionRateIonChannel::Tangent_ReactionRate_Ce(FEMaterialPoint double R = GetFEModel()->GetGlobalConstant("R"); double T = GetFEModel()->GetGlobalConstant("T"); double Fc = GetFEModel()->GetGlobalConstant("Fc"); - + FEMultiphasic* pMP = m_pReact->m_pMP; + assert(pMP); + double ksi = pMP->SBMArealConcentration(pt, m_sbm); + double dkdc = 0; if ((ci > 0) && (ce > 0)) - dkdc = (ci != ce) ? R*T*m_g/pow(Fc*m_z,2)*(ce*(1+log(ci/ce))-ci)/pow(ci-ce,2)/ce : -R*T*m_g/pow(ci*Fc*m_z,2)/2; + dkdc = (ci != ce) ? R*T*m_g/ksi/pow(Fc*m_z,2)*(ce*(1+log(ci/ce))-ci)/pow(ci-ce,2)/ce : -R*T*m_g/pow(ci*Fc*m_z,2)/2; return dkdc; } diff --git a/FEBioMix/FEMembraneReactionRateIonChannel.h b/FEBioMix/FEMembraneReactionRateIonChannel.h index 9eeb7cea8..8285aa359 100644 --- a/FEBioMix/FEMembraneReactionRateIonChannel.h +++ b/FEBioMix/FEMembraneReactionRateIonChannel.h @@ -33,7 +33,7 @@ class FEBIOMIX_API FEMembraneReactionRateIonChannel : public FEMembraneReactionR { public: //! constructor - FEMembraneReactionRateIonChannel(FEModel* pfem) : FEMembraneReactionRate(pfem) { m_g = 0; m_sol = -1; m_z = 0; } + FEMembraneReactionRateIonChannel(FEModel* pfem) : FEMembraneReactionRate(pfem) { m_g = 0; m_sol = -1; m_z = 0; m_sbm = -1;} // initialization bool Init() override; @@ -42,7 +42,7 @@ class FEBIOMIX_API FEMembraneReactionRateIonChannel : public FEMembraneReactionR double ReactionRate(FEMaterialPoint& pt) override; //! tangent of reaction rate with strain at material point - double Tangent_ReactionRate_Strain(FEMaterialPoint& pt) override { return 0; } + double Tangent_ReactionRate_Strain(FEMaterialPoint& pt) override; //! tangent of reaction rate with effective fluid pressure at material point double Tangent_ReactionRate_Pressure(FEMaterialPoint& pt) override {return 0; } @@ -55,8 +55,9 @@ class FEBIOMIX_API FEMembraneReactionRateIonChannel : public FEMembraneReactionR double Tangent_ReactionRate_Ci(FEMaterialPoint& pt, const int isol) override; public: - int m_sol; //!< solute id (1-based) + int m_sol; //!< solute id (1-based) for ion int m_z; //!< charge number of channel ion + int m_sbm; //!< sbm id (1-based) for channel protein double m_g; //!< channel conductance diff --git a/FEBioMix/FEMichaelisMenten.cpp b/FEBioMix/FEMichaelisMenten.cpp index 5906c97e4..5b6f6dcee 100644 --- a/FEBioMix/FEMichaelisMenten.cpp +++ b/FEBioMix/FEMichaelisMenten.cpp @@ -131,7 +131,7 @@ mat3ds FEMichaelisMenten::Tangent_ReactionSupply_Strain(FEMaterialPoint& pt) double J = ept.m_J; double phi0 = 0.0; if (m_pMP) - phi0 = bpt.m_phi0; + phi0 = bpt.m_phi0t; else if (m_pSM || m_pMF || m_pFS) phi0 = bfpt.m_phi0; dcdJ = -c/(J-phi0); diff --git a/FEBioMix/FEMultiphasic.cpp b/FEBioMix/FEMultiphasic.cpp index 605d05268..f35bd236b 100644 --- a/FEBioMix/FEMultiphasic.cpp +++ b/FEBioMix/FEMultiphasic.cpp @@ -422,7 +422,7 @@ double FEMultiphasic::SolidReferentialApparentDensity(FEMaterialPoint& pt) // evaluate referential apparent density of base solid double density = m_pSolid->Density(pt); - double rhosr = pet.m_phi0*density; + double rhosr = pet.m_phi0t*density; // add contribution from solid-bound molecules for (int isbm=0; isbm<(int)spt.m_sbmr.size(); ++isbm) @@ -453,7 +453,7 @@ double FEMultiphasic::Porosity(FEMaterialPoint& pt) FEBiphasicMaterialPoint& bt = *pt.ExtractData(); // solid referential volume fraction - double phisr = bt.m_phi0; + double phisr = bt.m_phi0t; // relative volume double J = et.m_J; @@ -476,7 +476,8 @@ double FEMultiphasic::FixedChargeDensity(FEMaterialPoint& pt) // relative volume double J = et.m_J; - double phi0 = bt.m_phi0; + double phi0 = bt.m_phi0; + double phir = bt.m_phi0t; double ce = 0; // add contribution from charged solid-bound molecules @@ -484,7 +485,7 @@ double FEMultiphasic::FixedChargeDensity(FEMaterialPoint& pt) ce += SBMChargeNumber(isbm)*spt.m_sbmr[isbm]/SBMMolarMass(isbm); double cFr = m_cFr(pt); - double cF = (cFr*(1-phi0)+ce)/(J-phi0); + double cF = (cFr*(1-phi0)+ce)/(J-phir); return cF; } @@ -625,7 +626,7 @@ void FEMultiphasic::PartitionCoefficientFunctions(FEMaterialPoint& mp, vector(); FEBiphasicMaterialPoint& bpt = *pt.ExtractData(); FESolutesMaterialPoint& spt = *pt.ExtractData(); - return spt.m_sbmr[sbm]/(ept.m_J-bpt.m_phi0)/SBMMolarMass(sbm); + return spt.m_sbmr[sbm]/(ept.m_J-bpt.m_phi0t)/SBMMolarMass(sbm); } //! SBM areal concentration (mole per shell area) -- should only be called from shell domains diff --git a/FEBioMix/FEMultiphasicDomain.cpp b/FEBioMix/FEMultiphasicDomain.cpp index 50a79d71b..2e9dfcd5e 100644 --- a/FEBioMix/FEMultiphasicDomain.cpp +++ b/FEBioMix/FEMultiphasicDomain.cpp @@ -39,4 +39,6 @@ FEMultiphasicDomain::FEMultiphasicDomain(FEModel* pfem) : FEElasticDomain(pfem) m_dofC = pfem->GetDOFIndex("concentration", 0); m_dofD = pfem->GetDOFIndex("shell concentration", 0); + + m_breset = false; } diff --git a/FEBioMix/FEMultiphasicDomain.h b/FEBioMix/FEMultiphasicDomain.h index d12f401b5..022cfc8ce 100644 --- a/FEBioMix/FEMultiphasicDomain.h +++ b/FEBioMix/FEMultiphasicDomain.h @@ -73,4 +73,7 @@ class FEBIOMIX_API FEMultiphasicDomain : public FEElasticDomain int m_dofVX; int m_dofVY; int m_dofVZ; + +protected: + bool m_breset; //!< flag for when Reset() is called }; diff --git a/FEBioMix/FEMultiphasicMultigeneration.cpp b/FEBioMix/FEMultiphasicMultigeneration.cpp index 5277d24be..05c327cf0 100644 --- a/FEBioMix/FEMultiphasicMultigeneration.cpp +++ b/FEBioMix/FEMultiphasicMultigeneration.cpp @@ -86,7 +86,7 @@ void FEMultiphasicMultigeneration::UpdateSolidBoundMolecules(FEMaterialPoint& mp FESolutesMaterialPoint& spt = *(mp.ExtractData()); FEMultigenSBMMaterialPoint& mpt = *(mp.ExtractData()); - double phi0 = ppt.m_phi0; + double phi0 = ppt.m_phi0t; int nsbm = SBMs(); int nsol = Solutes(); int ngen = mpt.m_ngen; diff --git a/FEBioMix/FEMultiphasicShellDomain.cpp b/FEBioMix/FEMultiphasicShellDomain.cpp index 4b81616dd..f7dc6b03c 100644 --- a/FEBioMix/FEMultiphasicShellDomain.cpp +++ b/FEBioMix/FEMultiphasicShellDomain.cpp @@ -191,9 +191,6 @@ bool FEMultiphasicShellDomain::Init() const int nsbm = m_pMat->SBMs(); const int nsol = m_pMat->Solutes(); vector sbmr(nsbm, 0); - for (int i = 0; iGetSBM(i)->m_rho0; - } for (int i = 0; i<(int)m_Elem.size(); ++i) { @@ -211,12 +208,17 @@ bool FEMultiphasicShellDomain::Init() FEMaterialPoint& mp = *el.GetMaterialPoint(n); FEBiphasicMaterialPoint& pb = *(mp.ExtractData()); FESolutesMaterialPoint& ps = *(mp.ExtractData()); + double h = el.Evaluate(el.m_ht, n); // shell thickness + + // extract the initial apparent densities of the solid-bound molecules + for (int i = 0; iGetSBM(i)->m_rho0(mp); ps.m_sbmr = sbmr; ps.m_sbmrp = sbmr; ps.m_sbmrhat.assign(nsbm, 0); ps.m_sbmrhatp.assign(nsbm, 0); - pb.m_phi0 = m_pMat->SolidReferentialVolumeFraction(mp); + pb.m_phi0t = m_pMat->SolidReferentialVolumeFraction(mp); ps.m_cF = m_pMat->FixedChargeDensity(mp); // evaluate reaction rates at initial time @@ -228,7 +230,7 @@ bool FEMultiphasicShellDomain::Init() // multiphasic material point data FEElasticMaterialPoint& pt = *(mp.ExtractData()); - double phi0 = pb.m_phi0; + double phi0 = pb.m_phi0t; for (int isbm=0; isbm()); - double phi0 = pb.m_phi0; for (int isbm=0; isbmGetMembraneReaction(k)->ReactionSupply(mp); double v = m_pMat->GetMembraneReaction(k)->m_v[nsol+isbm]; // remember to convert from molar supply to referential mass supply - ps.m_sbmrhat[isbm] += (pt.m_J-phi0)*m_pMat->SBMMolarMass(isbm)*v*zetahat; + ps.m_sbmrhat[isbm] += pt.m_J/h*m_pMat->SBMMolarMass(isbm)*v*zetahat; } } } @@ -376,7 +377,7 @@ void FEMultiphasicShellDomain::InitMaterialPoints() pt.m_pa = m_pMat->Pressure(mp); // initialize referential solid volume fraction - pt.m_phi0 = m_pMat->SolidReferentialVolumeFraction(mp); + pt.m_phi0t = m_pMat->SolidReferentialVolumeFraction(mp); // calculate FCD, current and stress ps.m_cF = m_pMat->FixedChargeDensity(mp); @@ -394,12 +395,7 @@ void FEMultiphasicShellDomain::Reset() const int nsol = m_pMat->Solutes(); const int nsbm = m_pMat->SBMs(); - - // extract the initial concentrations of the solid-bound molecules vector sbmr(nsbm,0); - for (int i=0; iGetSBM(i)->m_rho0; - } for (int i=0; i<(int) m_Elem.size(); ++i) { @@ -417,7 +413,11 @@ void FEMultiphasicShellDomain::Reset() FESolutesMaterialPoint& ps = *(mp.ExtractData()); // initialize referential solid volume fraction - pt.m_phi0 = m_pMat->m_phi0(mp); + pt.m_phi0 = pt.m_phi0t = m_pMat->m_phi0(mp); + + // extract the initial apparent densities of the solid-bound molecules + for (int i = 0; iGetSBM(i)->m_rho0(mp); // initialize multiphasic solutes ps.m_nsol = nsol; @@ -454,6 +454,7 @@ void FEMultiphasicShellDomain::Reset() } } } + m_breset = true; } //----------------------------------------------------------------------------- @@ -496,7 +497,7 @@ void FEMultiphasicShellDomain::PreSolveUpdate(const FETimeInfo& timeInfo) pt.m_Jp = pe.m_J; // reset referential solid volume fraction at previous time - pt.m_phi0p = pt.m_phi0; + pt.m_phi0p = pt.m_phi0t; // reset referential actual solute concentration at previous time for (int k=0; kSolutes(); ++k) { @@ -1036,7 +1037,7 @@ bool FEMultiphasicShellDomain::ElementMultiphasicStiffness(FEShellElement& el, m // evaluate the porosity and its derivative double phiw = m_pMat->Porosity(mp); - double phi0 = ppt.m_phi0; + double phi0 = ppt.m_phi0t; double phis = 1. - phiw; double dpdJ = phis/J; @@ -2377,8 +2378,9 @@ void FEMultiphasicShellDomain::UpdateElementStress(int iel, const FETimeInfo& tp pmb->UpdateSolidBoundMolecules(mp); // evaluate referential solid volume fraction - ppt.m_phi0 = pmb->SolidReferentialVolumeFraction(mp); - + ppt.m_phi0t = pmb->SolidReferentialVolumeFraction(mp); + if (m_breset) ppt.m_phi0 = ppt.m_phi0t; + // evaluate fluid pressure at gauss-point ppt.m_p = evaluate(el, pn, qn, n); @@ -2420,6 +2422,7 @@ void FEMultiphasicShellDomain::UpdateElementStress(int iel, const FETimeInfo& tp pmb->GetMembraneReaction(j)->UpdateElementData(mp); } + if (m_breset) m_breset = false; } //----------------------------------------------------------------------------- diff --git a/FEBioMix/FEMultiphasicSolidDomain.cpp b/FEBioMix/FEMultiphasicSolidDomain.cpp index 8edbe0d81..417a8773a 100644 --- a/FEBioMix/FEMultiphasicSolidDomain.cpp +++ b/FEBioMix/FEMultiphasicSolidDomain.cpp @@ -136,9 +136,6 @@ bool FEMultiphasicSolidDomain::Init() const int nsbm = m_pMat->SBMs(); const int nsol = m_pMat->Solutes(); vector sbmr(nsbm, 0); - for (int i = 0; iGetSBM(i)->m_rho0; - } for (int i = 0; i<(int)m_Elem.size(); ++i) { @@ -155,11 +152,13 @@ bool FEMultiphasicSolidDomain::Init() FEBiphasicMaterialPoint& pb = *(mp.ExtractData()); FESolutesMaterialPoint& ps = *(mp.ExtractData()); + for (int i = 0; iGetSBM(i)->m_rho0(mp); ps.m_sbmr = sbmr; ps.m_sbmrp.assign(nsbm, 0); ps.m_sbmrhat.assign(nsbm, 0); ps.m_sbmrhatp.assign(nsbm, 0); - pb.m_phi0 = m_pMat->SolidReferentialVolumeFraction(mp); + pb.m_phi0t = m_pMat->SolidReferentialVolumeFraction(mp); ps.m_cF = m_pMat->FixedChargeDensity(mp); // evaluate reaction rates at initial time @@ -171,7 +170,7 @@ bool FEMultiphasicSolidDomain::Init() // multiphasic material point data FEElasticMaterialPoint& pt = *(mp.ExtractData()); - double phi0 = pb.m_phi0; + double phi0 = pb.m_phi0t; for (int isbm=0; isbmPressure(mp); // initialize referential solid volume fraction - pt.m_phi0 = m_pMat->SolidReferentialVolumeFraction(mp); + pt.m_phi0t = m_pMat->SolidReferentialVolumeFraction(mp); // calculate FCD, current and stress ps.m_cF = m_pMat->FixedChargeDensity(mp); @@ -368,9 +367,6 @@ void FEMultiphasicSolidDomain::Reset() // extract the initial concentrations of the solid-bound molecules vector sbmr(nsbm,0); - for (int i=0; iGetSBM(i)->m_rho0; - } for (int i=0; i<(int) m_Elem.size(); ++i) { @@ -388,7 +384,11 @@ void FEMultiphasicSolidDomain::Reset() FESolutesMaterialPoint& ps = *(mp.ExtractData()); // initialize referential solid volume fraction - pt.m_phi0 = m_pMat->m_phi0(mp); + pt.m_phi0 = pt.m_phi0t = m_pMat->m_phi0(mp); + + // initialize sbm apparent densities + for (int i = 0; iGetSBM(i)->m_rho0(mp); // initialize multiphasic solutes ps.m_nsol = nsol; @@ -413,6 +413,8 @@ void FEMultiphasicSolidDomain::Reset() m_pMat->GetReaction(j)->ResetElementData(mp); } } + + m_breset = true; } //----------------------------------------------------------------------------- @@ -454,7 +456,7 @@ void FEMultiphasicSolidDomain::PreSolveUpdate(const FETimeInfo& timeInfo) pt.m_Jp = pe.m_J; // reset referential solid volume fraction at previous time - pt.m_phi0p = pt.m_phi0; + pt.m_phi0p = pt.m_phi0t; // reset referential actual solute concentration at previous time for (int j=0; jSolutes(); ++j) { @@ -924,7 +926,7 @@ bool FEMultiphasicSolidDomain::ElementMultiphasicStiffness(FESolidElement& el, m // evaluate the porosity and its derivative double phiw = m_pMat->Porosity(mp); - double phi0 = ppt.m_phi0; + double phi0 = ppt.m_phi0t; double phis = 1. - phiw; double dpdJ = phis/J; @@ -1669,7 +1671,8 @@ void FEMultiphasicSolidDomain::UpdateElementStress(int iel, double dt) pmb->UpdateSolidBoundMolecules(mp); // evaluate referential solid volume fraction - ppt.m_phi0 = pmb->SolidReferentialVolumeFraction(mp); + ppt.m_phi0t = pmb->SolidReferentialVolumeFraction(mp); + if (m_breset) ppt.m_phi0 = ppt.m_phi0t; // evaluate fluid pressure at gauss-point ppt.m_p = el.Evaluate(pn, n); @@ -1715,4 +1718,6 @@ void FEMultiphasicSolidDomain::UpdateElementStress(int iel, double dt) pmb->GetReaction(j)->UpdateElementData(mp); } + + if (m_breset) m_breset = false; } diff --git a/FEBioMix/FEMultiphasicStandard.cpp b/FEBioMix/FEMultiphasicStandard.cpp index abd95ebe3..d151bfa20 100644 --- a/FEBioMix/FEMultiphasicStandard.cpp +++ b/FEBioMix/FEMultiphasicStandard.cpp @@ -57,8 +57,10 @@ void FEMultiphasicStandard::UpdateSolidBoundMolecules(FEMaterialPoint& mp) FEElasticMaterialPoint& pt = *(mp.ExtractData()); FEBiphasicMaterialPoint& ppt = *(mp.ExtractData()); FESolutesMaterialPoint& spt = *(mp.ExtractData()); - - double phi0 = ppt.m_phi0; + FEShellElement* sel = dynamic_cast(mp.m_elem); + double h = (sel) ? sel->Evaluate(sel->m_ht, mp.m_index) : 0; // shell thickness + + double phi0 = ppt.m_phi0t; int nsbm = SBMs(); int nsol = Solutes(); // create a temporary container for spt.m_sbmr so that this variable remains @@ -77,7 +79,7 @@ void FEMultiphasicStandard::UpdateSolidBoundMolecules(FEMaterialPoint& mp) double zetahat = GetMembraneReaction(k)->ReactionSupply(mp); double v = GetMembraneReaction(k)->m_v[nsol+isbm]; // remember to convert from molar supply to referential mass supply - spt.m_sbmrhat[isbm] += (pt.m_J-phi0)*SBMMolarMass(isbm)*v*zetahat; + spt.m_sbmrhat[isbm] += pt.m_J/h*SBMMolarMass(isbm)*v*zetahat; } // perform the time integration (midpoint rule) sbmr[isbm] = spt.m_sbmrp[isbm] + dt*(spt.m_sbmrhat[isbm]+spt.m_sbmrhatp[isbm])/2; diff --git a/FEBioMix/FEOsmCoefManning.cpp b/FEBioMix/FEOsmCoefManning.cpp index cb608c628..86c317817 100644 --- a/FEBioMix/FEOsmCoefManning.cpp +++ b/FEBioMix/FEOsmCoefManning.cpp @@ -219,7 +219,7 @@ double FEOsmCoefManning::Tangent_OsmoticCoefficient_Strain_Manning(FEMaterialPoi double phisr = 0.0; if (m_pMP) - phisr = bpt.m_phi0; + phisr = bpt.m_phi0t; if (m_pMF) phisr = bfpt.m_phi0; diff --git a/FEBioMix/FEPermExpIso.cpp b/FEBioMix/FEPermExpIso.cpp index 45c2599ed..0488e137c 100644 --- a/FEBioMix/FEPermExpIso.cpp +++ b/FEBioMix/FEPermExpIso.cpp @@ -56,7 +56,7 @@ mat3ds FEPermExpIso::Permeability(FEMaterialPoint& mp) // referential solid volume fraction double phi0 = 0.0; if(pt) - phi0 = pt->m_phi0; + phi0 = pt->m_phi0t; else if (bpt) phi0 = bpt->m_phi0; @@ -80,7 +80,7 @@ tens4dmm FEPermExpIso::Tangent_Permeability_Strain(FEMaterialPoint &mp) // referential solid volume fraction double phi0 = 0.0; if(pt) - phi0 = pt->m_phi0; + phi0 = pt->m_phi0t; else if (bpt) phi0 = bpt->m_phi0; diff --git a/FEBioMix/FEPermHolmesMow.cpp b/FEBioMix/FEPermHolmesMow.cpp index b07536470..856abca87 100644 --- a/FEBioMix/FEPermHolmesMow.cpp +++ b/FEBioMix/FEPermHolmesMow.cpp @@ -28,6 +28,7 @@ SOFTWARE.*/ #include "stdafx.h" #include "FEPermHolmesMow.h" +#include // define the material parameters BEGIN_FECORE_CLASS(FEPermHolmesMow, FEHydraulicPermeability) @@ -57,10 +58,13 @@ mat3ds FEPermHolmesMow::Permeability(FEMaterialPoint& mp) // referential solid volume fraction also check if bfsi double phi0 = 0.0; if (pt) - phi0 = pt->m_phi0; + phi0 = pt->m_phi0t; else if (bpt) phi0 = bpt->m_phi0; + // check for potential error + if (J <= phi0) feLogError("The Holmes-Mow permeability calculation failed!\nThe volume ratio (J=%g) dropped below its theoretical minimum phi0=%g.",J,phi0); + // --- strain-dependent isotropic permeability --- return mat3dd(m_perm*pow((J-phi0)/(1.0-phi0),m_alpha)*exp(m_M*(J*J-1.0)/2.0)); @@ -79,10 +83,13 @@ tens4dmm FEPermHolmesMow::Tangent_Permeability_Strain(FEMaterialPoint &mp) // referential solid volume fraction double phi0 = 0.0; if (pt) - phi0 = pt->m_phi0; + phi0 = pt->m_phi0t; else if (bpt) phi0 = bpt->m_phi0; + // check for potential error + if (J <= phi0) feLogError("The Holmes-Mow permeability calculation failed!\nThe volume ratio (J=%g) dropped below its theoretical minimum phi0=%g.",J,phi0); + mat3dd I(1); // Identity double k0 = m_perm*pow((J-phi0)/(1.0-phi0),m_alpha)*exp(m_M*(J*J-1.0)/2.0); diff --git a/FEBioMix/FEPermRefIso.cpp b/FEBioMix/FEPermRefIso.cpp index e33719d6b..138555ba9 100644 --- a/FEBioMix/FEPermRefIso.cpp +++ b/FEBioMix/FEPermRefIso.cpp @@ -79,7 +79,7 @@ mat3ds FEPermRefIso::Permeability(FEMaterialPoint& mp) // relative volume double J = et.m_J; // referential solid volume fraction - double phi0 = pt.m_phi0; + double phi0 = pt.m_phi0t; // --- strain-dependent permeability --- @@ -108,7 +108,7 @@ tens4dmm FEPermRefIso::Tangent_Permeability_Strain(FEMaterialPoint &mp) // relative volume double J = et.m_J; // referential solid volume fraction - double phi0 = pt.m_phi0; + double phi0 = pt.m_phi0t; double f = pow((J-phi0)/(1-phi0),m_alpha)*exp(m_M*(J*J-1.0)/2.0); double k0 = m_perm0*f; diff --git a/FEBioMix/FEPermRefOrtho.cpp b/FEBioMix/FEPermRefOrtho.cpp index 4ebc38fd5..52e554f42 100644 --- a/FEBioMix/FEPermRefOrtho.cpp +++ b/FEBioMix/FEPermRefOrtho.cpp @@ -79,7 +79,7 @@ mat3ds FEPermRefOrtho::Permeability(FEMaterialPoint& mp) // relative volume double J = et.m_J; // referential solid volume fraction - double phi0 = pt.m_phi0; + double phi0 = pt.m_phi0t; // get the local coordinate systems mat3d Q = GetLocalCS(mp); @@ -135,7 +135,7 @@ tens4dmm FEPermRefOrtho::Tangent_Permeability_Strain(FEMaterialPoint &mp) // relative volume double J = et.m_J; // referential solid volume fraction - double phi0 = pt.m_phi0; + double phi0 = pt.m_phi0t; // get local coordinates mat3d Q = GetLocalCS(mp); diff --git a/FEBioMix/FEPermRefTransIso.cpp b/FEBioMix/FEPermRefTransIso.cpp index e1e5faad5..533fb6a26 100644 --- a/FEBioMix/FEPermRefTransIso.cpp +++ b/FEBioMix/FEPermRefTransIso.cpp @@ -78,7 +78,7 @@ mat3ds FEPermRefTransIso::Permeability(FEMaterialPoint& mp) // relative volume double J = et.m_J; // referential solid volume fraction - double phi0 = pt.m_phi0; + double phi0 = pt.m_phi0t; // get the local coordinate systems mat3d Q = GetLocalCS(mp); @@ -127,7 +127,7 @@ tens4dmm FEPermRefTransIso::Tangent_Permeability_Strain(FEMaterialPoint &mp) // relative volume double J = et.m_J; // referential solid volume fraction - double phi0 = pt.m_phi0; + double phi0 = pt.m_phi0t; // get the local coordinate systems mat3d Q = GetLocalCS(mp); diff --git a/FEBioMix/FEPorousNeoHookean.cpp b/FEBioMix/FEPorousNeoHookean.cpp index 676f1800c..a7e22d3d9 100644 --- a/FEBioMix/FEPorousNeoHookean.cpp +++ b/FEBioMix/FEPorousNeoHookean.cpp @@ -128,5 +128,5 @@ double FEPorousNeoHookean::ReferentialSolidVolumeFraction(FEMaterialPoint& mp) if (m_phisr < 1) return m_phisr; FEBiphasicMaterialPoint& pt = *mp.ExtractData(); - return pt.m_phi0; + return pt.m_phi0t; } diff --git a/FEBioMix/FESBMPointSource.cpp b/FEBioMix/FESBMPointSource.cpp index 99427253d..1d81bb5fc 100644 --- a/FEBioMix/FESBMPointSource.cpp +++ b/FEBioMix/FESBMPointSource.cpp @@ -26,42 +26,120 @@ SOFTWARE.*/ #include "stdafx.h" #include "FESBMPointSource.h" #include "FEMultiphasic.h" +#include "FESolute.h" #include +#include +#include #include +#include #include +#include +#include +#include BEGIN_FECORE_CLASS(FESBMPointSource, FEBodyLoad) - ADD_PARAMETER(m_sbm, "sbm"); + ADD_PARAMETER(m_sbmId, "sbm"); + ADD_PARAMETER(m_rate, "rate"); ADD_PARAMETER(m_pos.x, "x"); ADD_PARAMETER(m_pos.y, "y"); ADD_PARAMETER(m_pos.z, "z"); - ADD_PARAMETER(m_val, "value"); ADD_PARAMETER(m_weighVolume, "weigh_volume"); END_FECORE_CLASS(); FESBMPointSource::FESBMPointSource(FEModel* fem) : FEBodyLoad(fem), m_search(&fem->GetMesh()) { - static bool bfirst = true; - m_sbm = -1; - m_pos = vec3d(0,0,0); - m_val = 0.0; - m_reset = bfirst; - m_doReset = true; + //static bool bfirst = true; + m_sbmId = -1; + m_pos = vec3d(0, 0, 0); + m_rate = 0.0; + //m_reset = bfirst; + //m_doReset = true; m_weighVolume = true; - bfirst = false; + //bfirst = false; } +vec3d FESBMPointSource::GetPosition() const +{ + return m_pos; +} + +void FESBMPointSource::SetPosition(const vec3d& pos) +{ + m_pos = pos; +} + +int FESBMPointSource::GetSBMID() const +{ + return m_sbmId; +} + + +void FESBMPointSource::SetSBMID(int sbmID) +{ + m_sbmId = sbmID; +} + +double FESBMPointSource::GetRate() const +{ + return m_rate; +} + +void FESBMPointSource::SetRate(double rate) +{ + m_rate = rate; +} + +double FESBMPointSource::GetdC() const +{ + return m_dC; +} + +double FESBMPointSource::GetdCp() const +{ + return m_dCp; +} + +void FESBMPointSource::SetdC(double dC) +{ + m_dC = dC; +} + + +void FESBMPointSource::SetRadius(double radius) +{ + m_radius = radius; + m_Vc = (4.0 / 3.0) * PI * pow(m_radius, 3.0); +} + +void FESBMPointSource::SetAccumulateFlag(bool b) +{ + m_accumulate = b; +} + +//void FESBMPointSource::SetWeighVolume(bool b) +//{ +// m_weighVolume = b; +//} + +//void FESBMPointSource::SetResetFlag(bool b) +//{ +// m_doReset = b; +//} + + + bool FESBMPointSource::Init() { - if (m_sbm == -1) return false; + if (m_sbmId == -1) return false; if (m_search.Init() == false) return false; + m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); return FEBodyLoad::Init(); } // allow species to accumulate at the point source void FESBMPointSource::Accumulate(double dc) { - double rt[3] = { 0, 0, 0 }; - m_el = dynamic_cast(m_search.FindElement(m_pos, rt)); + // find the element in which the point lies + m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); if (m_el == nullptr) return; // make sure this element is part of a multiphasic domain @@ -69,12 +147,13 @@ void FESBMPointSource::Accumulate(double dc) { FEMultiphasic* mat = dynamic_cast(dom->GetMaterial()); if (mat == nullptr) return; + // Make sure the material has the correct solute int sbmid = -1; int sbms = mat->SBMs(); for (int j = 0; j < sbms; ++j) { int sbmj = mat->GetSBM(j)->GetSBMID(); - if (sbmj == m_sbm) + if (sbmj == m_sbmId) { sbmid = j; break; @@ -82,18 +161,17 @@ void FESBMPointSource::Accumulate(double dc) { } if (sbmid == -1) return; - m_val = dc + m_val; // prevent negative concentrations + m_rate = dc + m_rate; m_accumulate = true; } void FESBMPointSource::Update() { - if (m_reset && m_doReset) ResetSBM(); - + //if (m_reset && m_doReset) ResetSBM(); if (m_accumulate) { + m_dCp = m_dC; // find the element in which the point lies - double rt[3] = { 0, 0, 0 }; - m_el = dynamic_cast(m_search.FindElement(m_pos, rt)); + m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); if (m_el == nullptr) return; // make sure this element is part of a multiphasic domain @@ -107,167 +185,258 @@ void FESBMPointSource::Update() // we prescribe the element average to the integration points const int nint = m_el->GaussPoints(); - double val = (m_weighVolume ? m_val / Ve : m_val); - // Make sure the material has the correct sbm int sbmid = -1; int sbms = mat->SBMs(); for (int j = 0; j < sbms; ++j) { int sbmj = mat->GetSBM(j)->GetSBMID(); - if (sbmj == m_sbm) + if (sbmj == m_sbmId) { sbmid = j; break; } } if (sbmid == -1) return; - std::vector possible_ints = FindIntInRadius(); - if (possible_ints.size() == 0) { + + double m_dt = mesh->GetFEModel()->GetCurrentStep()->m_dt; + std::vector possible_ints; + double total_elem = 0; + FindIntInRadius(possible_ints, total_elem); + double total_change = 0.0; + // if the cell is not projecting on top of an integration point (i.e. too small) + // then project it's concentration to the nodes via the shape functions + //SL TODO: Currently assigns just to the current element. Should instead find the ~8 closest integration points + // regardless of element. + if (possible_ints.size() == 0) + { + FindNodesInRadius(possible_ints, total_elem); + } + if (possible_ints.size() == 0) + { // set the concentration of all the integration points double H[FEElement::MAX_NODES]; - double m_q[3]; - m_q[0] = m_q[1] = m_q[2] = 0.0; - m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); - if (m_el == nullptr) return; m_el->shape_fnc(H, m_q[0], m_q[1], m_q[2]); - for (int i = 0; i < nint; ++i) { + double H[FEElement::MAX_NODES]; + m_el->shape_fnc(H, m_q[0], m_q[1], m_q[2]); FEMaterialPoint& mp = *m_el->GetMaterialPoint(i); FESolutesMaterialPoint& pd = *(mp.ExtractData()); - pd.m_sbmr[sbmid] = std::max(0.0, nint * H[i] * val + pd.m_sbmrp[sbmid]); + double new_r = H[i] * m_rate * nint / Ve + pd.m_sbmr[sbmid]; + pd.m_sbmr[sbmid] = new_r; pd.m_sbmrp[sbmid] = pd.m_sbmr[sbmid]; + total_change += m_rate / nint; } } + // else evenly distribute it among the integration points that the cell is on top of. + //SL TODO: currently may lead to a little bias when neighboring elements are smaller resulting in + // higher density of integration points else { int nint_in = possible_ints.size(); - double m_q[3]; - m_q[0] = m_q[1] = m_q[2] = 0.0; - m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); - if (m_el == nullptr) return; for (auto iter = possible_ints.begin(); iter != possible_ints.end(); ++iter) { FEMaterialPoint& mp = **iter; FESolutesMaterialPoint& pd = *(mp.ExtractData()); - pd.m_sbmr[sbmid] = std::max(0.0, nint * val / nint_in + pd.m_sbmrp[sbmid]); + // scale by H so that the total integral over each element is consistent. + double H = double(nint) / double(nint_in); + double new_r = H * m_rate / Ve + pd.m_sbmr[sbmid]; + pd.m_sbmr[sbmid] = new_r; pd.m_sbmrp[sbmid] = pd.m_sbmr[sbmid]; + total_change += m_rate / (nint_in); } - } - - // multiply by # integration points to prevent smoothing/distilling concentration + } + m_dC = -total_change / m_Vc; m_accumulate = false; // don't double count a point source - m_val = 0; } } -void FESBMPointSource::SetPosition(const vec3d& pos) -{ - m_pos = pos; -} - -vec3d FESBMPointSource::GetPosition() const -{ - return m_pos; -} - -void FESBMPointSource::SetSBMID(int id) -{ - m_sbm = id; -} - -int FESBMPointSource::GetSBMID() const -{ - return m_sbm; -} - -void FESBMPointSource::SetValue(double val) -{ - m_val = val; -} - -void FESBMPointSource::SetRadius(double radius) -{ - m_radius = radius; -} - -double FESBMPointSource::GetValue() const -{ - return m_val; -} - -void FESBMPointSource::SetWeighVolume(bool b) -{ - m_weighVolume = b; -} - -void FESBMPointSource::SetResetFlag(bool b) -{ - m_doReset = b; -} +void FESBMPointSource::FindIntInRadius(std::vector &possible_ints, double &total_elem) { + + // get element and set up buffers + m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); + std::unordered_set visited; + std::set next; + //std::vector possible_ints; + visited.reserve(1000); + possible_ints.reserve(500); + + //we will need to check the current element first + next.insert(m_el); + // create the element adjacency list. + FEDomain* dom = dynamic_cast(m_el->GetMeshPartition()); + FEMultiphasic* mat = dynamic_cast(dom->GetMaterial()); + if (mat == nullptr) return; + // calculate the element volume + auto mesh = dom->GetMesh(); + // create the element-element list + FEElemElemList EEL; + EEL.Create(mesh); + + //while there are still elements to evaluate + while (next.size()) { + // get the element to be evaluated + FESolidElement* cur = *next.begin(); + // remove the current element from the next buffer and add it to the visited buffer + next.erase(next.begin()); + visited.insert(cur); + // get the current element bounds + std::vector cur_element_bounds; + // add integration points within the radius + bool int_flag = false; + for (int i = 0; i < cur->GaussPoints(); i++) + { + FEMaterialPoint* mp = cur->GetMaterialPoint(i); + vec3d disp = mp->m_r0 - m_pos; + if (disp.norm() <= m_radius) { + possible_ints.push_back(mp); + int_flag = true; + } + } + if (int_flag) + { + total_elem++; + } -void FESBMPointSource::SetAccumulateFlag(bool b) -{ - m_accumulate = b; + // Add neighboring element to the next buffer as long as they haven't been visited. + // get the global ID of the current element + int cur_id = cur->GetID()-1; + // for each neighboring element + for (int i = 0; i < EEL.NeighborSize(); i++) + { + if (EEL.Neighbor(cur_id, i)) + { + // if that element has not been visited yet add it to the next list + if (!visited.count(dynamic_cast(EEL.Neighbor(cur_id, i)))) + { + next.insert(dynamic_cast(EEL.Neighbor(cur_id, i))); + } + } + } + } } -void FESBMPointSource::ResetSBM() -{ - FEModel& fem = *GetFEModel(); - FEMesh& mesh = fem.GetMesh(); +void FESBMPointSource::FindNodesInRadius(std::vector& possible_ints, double& total_elem) { - for (int i = 0; i < mesh.Domains(); ++i) - { - FEDomain& dom = mesh.Domain(i); - int NE = dom.Elements(); + // get element and set up buffers + m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); + std::unordered_set visited; + std::set next; + //std::vector possible_ints; + visited.reserve(1000); + std::vector possible_nodes; + possible_nodes.reserve(500); - FEMultiphasic* mat = dynamic_cast(dom.GetMaterial()); - if (mat) + //we will need to check the current element first + next.insert(m_el); + // create the element adjacency list. + FEDomain* dom = dynamic_cast(m_el->GetMeshPartition()); + FEMultiphasic* mat = dynamic_cast(dom->GetMaterial()); + if (mat == nullptr) return; + // calculate the element volume + auto mesh = dom->GetMesh(); + // create the element-element list + FEElemElemList EEL; + EEL.Create(mesh); + + //while there are still elements to evaluate + while (next.size()) { + // get the element to be evaluated + FESolidElement* cur = *next.begin(); + // remove the current element from the next buffer and add it to the visited buffer + next.erase(next.begin()); + visited.insert(cur); + // get the current element bounds + std::vector cur_element_bounds; + // add integration points within the radius + bool int_flag = false; + for (int i = 0; i < cur->Nodes(); i++) { - // Make sure the material has the correct sbm - int sbmid = -1; - int sbms = mat->SBMs(); - for (int j = 0; jGetSBM(j)->GetSBMID(); - if (sbmj == m_sbm) + FENode* mn = &(mesh->Node(cur->m_node[i])); + vec3d disp = mn->m_rt - m_pos; + FEMaterialPoint* closest; + if (disp.norm() <= m_radius) { + // find the closest mp in the element + double min_d = std::numeric_limits::max(); + for (int j = 0; j < cur->GaussPoints(); j++) { - sbmid = j; - break; + double disp2 = (mn->m_rt - cur->GetMaterialPoint(j)->m_rt).norm(); + if (disp2 < min_d) + { + min_d = disp2; + closest = cur->GetMaterialPoint(j); + } } + possible_ints.push_back(closest); + int_flag = true; } + } + if (int_flag) + { + total_elem++; + } - if (sbmid != -1) + // Add neighboring element to the next buffer as long as they haven't been visited. + // get the global ID of the current element + int cur_id = cur->GetID() - 1; + // for each neighboring element + for (int i = 0; i < EEL.NeighborSize(); i++) + { + if (EEL.Neighbor(cur_id, i)) { - for (int j = 0; j < NE; ++j) + // if that element has not been visited yet add it to the next list + if (!visited.count(dynamic_cast(EEL.Neighbor(cur_id, i)))) { - FEElement& el = dom.ElementRef(j); - - // set the concentration of all the integration points - int nint = el.GaussPoints(); - for (int k = 0; k < nint; ++k) - { - FEMaterialPoint* mp = el.GetMaterialPoint(k); - FESolutesMaterialPoint& pd = *(mp->ExtractData()); - pd.m_sbmr[sbmid] = 0.0; - pd.m_sbmrp[sbmid] = 0.0; - } + next.insert(dynamic_cast(EEL.Neighbor(cur_id, i))); } } } } } -std::vector FESBMPointSource::FindIntInRadius() { - // determine if the radius exceeds the boundaries of the element - int nint = m_el->GaussPoints(); - std::vector possible_ints; - for (int i = 0; i < nint; i++) { - FEMaterialPoint* mp = m_el->GetMaterialPoint(i); - vec3d disp = mp->m_r0 - m_pos; - if (disp.norm() <= m_radius) { - possible_ints.push_back(mp); - } - } - return possible_ints; -} \ No newline at end of file +//void FESBMPointSource::ResetSBM() +//{ +// FEModel& fem = *GetFEModel(); +// FEMesh& mesh = fem.GetMesh(); +// +// for (int i = 0; i < mesh.Domains(); ++i) +// { +// FEDomain& dom = mesh.Domain(i); +// int NE = dom.Elements(); +// +// FEMultiphasic* mat = dynamic_cast(dom.GetMaterial()); +// if (mat) +// { +// // Make sure the material has the correct sbm +// int sbmid = -1; +// int sbms = mat->SBMs(); +// for (int j = 0; jGetSBM(j)->GetSBMID(); +// if (sbmj == m_sbmId) +// { +// sbmid = j; +// break; +// } +// } +// +// if (sbmid != -1) +// { +// for (int j = 0; j < NE; ++j) +// { +// FEElement& el = dom.ElementRef(j); +// +// // set the concentration of all the integration points +// int nint = el.GaussPoints(); +// for (int k = 0; k < nint; ++k) +// { +// FEMaterialPoint* mp = el.GetMaterialPoint(k); +// FESolutesMaterialPoint& pd = *(mp->ExtractData()); +// pd.m_sbmr[sbmid] = 0.0; +// pd.m_sbmrp[sbmid] = 0.0; +// } +// } +// } +// } +// } +//} \ No newline at end of file diff --git a/FEBioMix/FESBMPointSource.h b/FEBioMix/FESBMPointSource.h index 49c2f991a..dbaa136e5 100644 --- a/FEBioMix/FESBMPointSource.h +++ b/FEBioMix/FESBMPointSource.h @@ -26,6 +26,7 @@ SOFTWARE.*/ #pragma once #include #include +#include class FESolidElement; @@ -48,11 +49,17 @@ class FESBMPointSource : public FEBodyLoad int GetSBMID() const; - void SetValue(double val); + void SetRate(double rate); void SetRadius(double radius); - double GetValue() const; + double GetRate() const; + + double GetdC() const; + + double GetdCp() const; + + void SetdC(double dC); void SetWeighVolume(bool b); @@ -60,24 +67,32 @@ class FESBMPointSource : public FEBodyLoad void SetAccumulateFlag(bool b); - std::vector FindIntInRadius(); + //std::vector FindIntInRadius(); + void FindIntInRadius(std::vector &possible_ints, double &total_elem); + + //! return all the elements in the given radius + void FindNodesInRadius(std::vector& possible_ints, double& total_elem); private: - void ResetSBM(); + //void ResetSBM(); private: - int m_sbm; // The SBM ID that defins the cell's "concentration" + int m_sbmId; // The SBM ID that defins the cell's "concentration" vec3d m_pos; // the position (in reference coordinates) - double m_val; // density value at point source + double m_rate; // density value at point source double m_radius; + double m_Vc; bool m_reset; bool m_doReset; bool m_weighVolume; bool m_accumulate; // accumulate species flag for the update + double m_dC = 0.0; // total change of a species + double m_dCp = 0.0; private: FEOctreeSearch m_search; FESolidElement* m_el; + double m_q[3]; DECLARE_FECORE_CLASS(); }; diff --git a/FEBioMix/FESolubManning.cpp b/FEBioMix/FESolubManning.cpp index 0925a5407..2a39f74b2 100644 --- a/FEBioMix/FESolubManning.cpp +++ b/FEBioMix/FESolubManning.cpp @@ -266,7 +266,7 @@ double FESolubManning::Tangent_Solubility_Strain_Manning(FEMaterialPoint &mp) double J = pt.m_J; double phisr = 0.0; if (m_pMP) - phisr = bpt.m_phi0; + phisr = bpt.m_phi0t; if (m_pMF) phisr = bfpt.m_phi0; diff --git a/FEBioMix/FESolute.h b/FEBioMix/FESolute.h index 171e38624..b54bc7b7b 100644 --- a/FEBioMix/FESolute.h +++ b/FEBioMix/FESolute.h @@ -276,7 +276,7 @@ class FESolidBoundMolecule : public FEMaterial double m_rhoT; //!< true SBM density double m_M; //!< SBM molar mass int m_z; //!< charge number of SBM - double m_rho0; //!< initial referential (apparent) density of SBM + FEParamDouble m_rho0; //!< initial referential (apparent) density of SBM double m_rhomin; //!< minimum referential (apparent) density of SBM double m_rhomax; //!< maximum referential (apparent) density of SBM diff --git a/FEBioMix/FESolutePointSource.cpp b/FEBioMix/FESolutePointSource.cpp index 47336866c..b43d64c94 100644 --- a/FEBioMix/FESolutePointSource.cpp +++ b/FEBioMix/FESolutePointSource.cpp @@ -29,8 +29,13 @@ SOFTWARE.*/ #include "FEMultiphasic.h" #include #include +#include #include #include +#include +#include +#include +#include BEGIN_FECORE_CLASS(FESolutePointSource, FEBodyLoad) ADD_PARAMETER(m_soluteId, "solute"); @@ -48,14 +53,19 @@ FESolutePointSource::FESolutePointSource(FEModel* fem) : FEBodyLoad(fem), m_sear m_rate = 0.0; } -void FESolutePointSource::SetPosition(const vec3d& v) +vec3d FESolutePointSource::GetPosition() const { - m_pos = v; + return m_pos; } -vec3d FESolutePointSource::GetPosition() const +void FESolutePointSource::SetPosition(const vec3d& pos) { - return m_pos; + m_pos = pos; +} + +int FESolutePointSource::GetSoluteID() const +{ + return m_soluteId; } void FESolutePointSource::SetSoluteID(int soluteID) @@ -63,9 +73,9 @@ void FESolutePointSource::SetSoluteID(int soluteID) m_soluteId = soluteID; } -int FESolutePointSource::GetSoluteID() const +double FESolutePointSource::GetRate() const { - return m_soluteId; + return m_rate; } void FESolutePointSource::SetRate(double rate) @@ -73,16 +83,41 @@ void FESolutePointSource::SetRate(double rate) m_rate = rate; } +double FESolutePointSource::GetdC() const +{ + return m_dC; +} + +double FESolutePointSource::GetdCp() const +{ + return m_dCp; +} + +void FESolutePointSource::SetdC(double dC) +{ + m_dC = dC; +} + +void FESolutePointSource::SetdCp(double dCp) +{ + m_dCp = dCp; +} + void FESolutePointSource::SetRadius(double radius) { m_radius = radius; + m_Vc = (4.0 / 3.0) * PI * pow(m_radius, 3.0); } -double FESolutePointSource::GetRate() const -{ - return m_rate; +void FESolutePointSource::SetAccumulateFlag(bool b) { + m_accumulate = b; +} + +void FESolutePointSource::SetAccumulateCAFlag(bool b) { + m_accumulate_ca = b; } + bool FESolutePointSource::Init() { // see if the solute exists @@ -107,12 +142,14 @@ bool FESolutePointSource::Init() // get the degree of freedom of the concentration m_dofC = fem->GetDOFIndex("concentration", m_soluteId - 1); + m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); + return FEBodyLoad::Init(); } +// allow species to accumulate at the point source void FESolutePointSource::Accumulate(double dc) { // find the element in which the point lies - m_q[0] = m_q[1] = m_q[2] = 0.0; m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); if (m_el == nullptr) return; @@ -121,8 +158,6 @@ void FESolutePointSource::Accumulate(double dc) { FEMultiphasic* mat = dynamic_cast(dom->GetMaterial()); if (mat == nullptr) return; - const int nint = m_el->GaussPoints(); - // Make sure the material has the correct solute int solid = -1; int sols = mat->Solutes(); @@ -139,78 +174,117 @@ void FESolutePointSource::Accumulate(double dc) { m_rate = dc + m_rate; m_accumulate = true; - m_accumulate_ca = true; - } void FESolutePointSource::Update() { - //if (m_accumulate_ca) { - // // find the element in which the point lies - // m_q[0] = m_q[1] = m_q[2] = 0.0; - // m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); - // if (m_el == nullptr) return; - - // // make sure this element is part of a multiphasic domain - // FEDomain* dom = dynamic_cast(m_el->GetMeshPartition()); - // FEMultiphasic* mat = dynamic_cast(dom->GetMaterial()); - // if (mat == nullptr) return; - - // // calculate the element volume - // FEMesh* mesh = dom->GetMesh(); - // double Ve = mesh->ElementVolume(*m_el); - - // // we prescribe the element average to the integration points - // const int nint = m_el->GaussPoints(); - // double val = m_rate / Ve; - - // // Make sure the material has the correct solute - // int solid = -1; - // int sols = mat->Solutes(); - // for (int j = 0; j < sols; ++j) - // { - // int solj = mat->GetSolute(j)->GetSoluteID(); - // if (solj == m_soluteId) - // { - // solid = j; - // break; - // } - // } - // if (solid == -1) return; - - // // set the concentration of all the integration points - // double H[FEElement::MAX_NODES]; - // double m_q[3]; - // m_q[0] = m_q[1] = m_q[2] = 0.0; - // FESolidElement* m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); - // if (m_el == nullptr) return; - // m_el->shape_fnc(H, m_q[0], m_q[1], m_q[2]); - - // for (int i = 0; i < nint; ++i) - // { - // FEMaterialPoint* mp = m_el->GetMaterialPoint(i); - // FESolutesMaterialPoint& pd = *(mp->ExtractData()); - // // if this point source has not yet been added to the integration points do it then turn off the flag - // if (m_accumulate_ca) { - // pd.m_ca[solid] = std::max(0.0, pd.m_crp[solid] + val); // prevent negative concentrations - // pd.m_crp[solid] = pd.m_ca[solid]; - // } - // else { - // pd.m_crp[solid] = pd.m_crp[solid]; - // pd.m_ca[solid] = pd.m_ca[solid]; - // } - - // m_accumulate_ca = false; - // m_rate = 0; - - // } - //} + if (m_accumulate) { + // find the element in which the point lies + m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); + if (m_el == nullptr) return; + + // make sure this element is part of a multiphasic domain + FESolidDomain* dom = dynamic_cast(m_el->GetMeshPartition()); + FEMultiphasic* mat = dynamic_cast(dom->GetMaterial()); + if (mat == nullptr) return; + + // calculate the element volume + FEMesh* mesh = dom->GetMesh(); + double Ve = mesh->ElementVolume(*m_el); + + std::cout << "rate is " << m_rate << endl; + + int solid = -1; + int sols = mat->Solutes(); + for (int j = 0; j < sols; ++j) { + int solj = mat->GetSolute(j)->GetSoluteID(); + if (solj == m_soluteId) + { + solid = j; + break; + } + } + if (solid == -1) return; + + // evaluate the concentration at this point + int neln = m_el->Nodes(); + double c[FEElement::MAX_NODES]; + for (int i = 0; i < neln; ++i) c[i] = mesh->Node(m_el->m_node[i]).get(m_dofC); + double cx = m_el->evaluate(c, m_q[0], m_q[1], m_q[2]); + double dt = mesh->GetFEModel()->GetCurrentStep()->m_dt; + double cxx = cx * Ve + dt * m_rate; + + // assemble the element load vector + vector fe(neln, 0.0); + vector lm(neln, -1); + std::vector possible_nodes; + double total_elem = 0.0; + FindNodesInRadius(possible_nodes, total_elem); + double total_change = 0.0; + //SL: Currently this just assigns within the current element. We will instead want this to assign to the closest integration points. + if (possible_nodes.size() == 0) { + // set the concentration of all nodes via the shape functions for each integration point + double H[FEElement::MAX_NODES]; + m_el->shape_fnc(H, m_q[0], m_q[1], m_q[2]); + int nint = m_el->GaussPoints(); double* w = m_el->GaussWeights(); + for (int n = 0; n < nint; ++n) + { + double* H_int = m_el->H(n); + FEMaterialPoint* mp = m_el->GetMaterialPoint(n); + double m_J = mp->m_J0; + // loop over all nodes + for (int j = 0; j < neln; ++j) + { + // only allow internalization if the species won't go negative + if (cxx > 0.0) { + total_change += m_rate * H[n] * H_int[j] * dt * w[n]; + } + } + } + } + else { + // set the concentration of all nodes via the shape functions for each integration point + double n_elem = possible_nodes.size(); + double H[FEElement::MAX_NODES]; + for (auto iter = possible_nodes.begin(); iter != possible_nodes.end(); iter++) + { + FESolidElement* c_el = dynamic_cast(*iter); + double r[3]; + dom->ProjectToElement(*c_el, m_pos, r); + vec3d r3 = ClampNatC(r); + r[0] = r3.x; r[1] = r3.y; r[2] = r3.z; + int nint = c_el->GaussPoints(); + double* w = c_el->GaussWeights(); + int neln = c_el->Nodes(); + // for this element get the shape functions and constants + c_el->shape_fnc(H, r[0], r[1], r[2]); + // for each integration point in this element project to the nodes. + for (int n = 0; n < nint; ++n) + { + double* H_int = c_el->H(n); + FEMaterialPoint* mp = c_el->GetMaterialPoint(n); + double m_J = mp->m_J0; + // loop over all nodes + for (int j = 0; j < neln; ++j) + { + // only allow internalization if the species won't go negative + if (cxx > 0.0) { + total_change += m_rate * H[n] * H_int[j] * dt * w[n] / n_elem; + } + } + } + } + } + m_dC = -total_change / (m_Vc); + m_accumulate = false; + } } //! Evaluate force vector void FESolutePointSource::LoadVector(FEGlobalVector& R, const FETimeInfo& tp) { // get the domain in which this element resides + m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); FESolidDomain* dom = dynamic_cast(m_el->GetMeshPartition()); FEMesh* mesh = dom->GetMesh(); @@ -227,19 +301,17 @@ void FESolutePointSource::LoadVector(FEGlobalVector& R, const FETimeInfo& tp) for (int i = 0; i < neln; ++i) c[i] = mesh->Node(m_el->m_node[i]).get(m_dofC); double cx = m_el->evaluate(c, m_q[0], m_q[1], m_q[2]); double Ve = mesh->ElementVolume(*m_el); + double cxx = cx * Ve + dt * m_rate; double v_rate = m_rate / Ve; - //if (m_rate > 0) { - // v_rate = m_rate / Ve; - //} - //else { - // v_rate = m_rate * Ve; - //} // assemble the element load vector vector fe(neln, 0.0); vector lm(neln, -1); - std::vector possible_ints = FindIntInRadius(); - if (possible_ints.size() == 0) { + std::vector possible_nodes; + double total_elem = 0; + FindNodesInRadius(possible_nodes, total_elem); + //SL: Currently this just assigns within the current element. We will instead want this to assign to the closest integration points. + if (possible_nodes.size() == 0) { int nint = m_el->GaussPoints(); double* w = m_el->GaussWeights(); for (int n = 0; n < nint; ++n) @@ -248,60 +320,43 @@ void FESolutePointSource::LoadVector(FEGlobalVector& R, const FETimeInfo& tp) FEMaterialPoint* mp = m_el->GetMaterialPoint(n); double m_J = mp->m_J0; // loop over all nodes - for (int j = 0; j < neln; ++j) - { - // get the value of the integrand for this node - // -ve needed since external forces are subtracted from internal forces - fe[j] += -v_rate * m_J * H[n] * H_int[j] * dt * w[n] * neln; - } + for (int j = 0; j < neln; ++j) fe[j] += -v_rate * m_J * H[n] * H_int[j] * dt * w[n] * neln; } - ////if (m_accumulate) { m_rate = 0; m_accumulate = false; } - //// get the LM vector for (int i = 0; i < neln; ++i) lm[i] = mesh->Node(m_el->m_node[i]).m_ID[m_dofC]; + R.Assemble(lm, fe); } - else { - int nint = m_el->GaussPoints(); - double* w = m_el->GaussWeights(); - for (int n = 0; n < nint; ++n) + else + { + double n_elem = possible_nodes.size(); + //vec3d global_pos = GetGlobalPos(vec3d(m_q[0],m_q[1],m_q[2]), m_el); + for (auto iter = possible_nodes.begin(); iter != possible_nodes.end(); iter++) { - double* H_int = m_el->H(n); - FEMaterialPoint* mp = m_el->GetMaterialPoint(n); - double m_J = mp->m_J0; - // loop over all nodes - for (int j = 0; j < neln; ++j) + FESolidElement* c_el = dynamic_cast(*iter); + //std::cout << "element " << c_el->GetID() << endl; + double r[3]; + dom->ProjectToElement(*c_el, m_pos, r); + vec3d r3 = ClampNatC(r); + r[0] = r3.x; r[1] = r3.y; r[2] = r3.z; + //std::cout << "r is " << r[0] << ", " << r[1] << ", " << r[2] << endl; + int nint = c_el->GaussPoints(); + double* w = c_el->GaussWeights(); + int neln = c_el->Nodes(); + vector fe(neln, 0.0); + vector lm(neln, -1); + for (int n = 0; n < nint; ++n) { - // get the value of the integrand for this node - // -ve needed since external forces are subtracted from internal forces - fe[j] += -v_rate * m_J * H[n] * H_int[j] * dt * w[n] * neln; + double* H_int = c_el->H(n); + c_el->shape_fnc(H, r[0], r[1], r[2]); + FEMaterialPoint* mp = c_el->GetMaterialPoint(n); + double m_J = mp->m_J0; + for (int j = 0; j < neln; ++j) fe[j] += -v_rate * m_J * H[n] * H_int[j] * dt * w[n] * neln / n_elem; } + for (int i = 0; i < neln; ++i) lm[i] = mesh->Node(c_el->m_node[i]).m_ID[m_dofC]; + R.Assemble(lm, fe); } - - ////if (m_accumulate) { m_rate = 0; m_accumulate = false; } - - //// get the LM vector - for (int i = 0; i < neln; ++i) lm[i] = mesh->Node(m_el->m_node[i]).m_ID[m_dofC]; } - // assemble into global vector - - //std::vector dist_2_pt; - //for (int i = 0; i < neln; ++i) - //{ - // FENode* m_n = &mesh->Node(m_el->m_node[i]); - // vec3d r = m_n->m_rt; - // vec3d disp_i = this->m_pos - r; - // dist_2_pt.push_back(disp_i.norm()); - //} - //int closest = std::min_element(dist_2_pt.begin(), dist_2_pt.end()) - dist_2_pt.begin(); - //vector lm(1, -1); - //lm[0] = mesh->Node(m_el->m_node[closest]).m_ID[m_dofC]; - //FEMaterialPoint* mp = m_el->GetMaterialPoint(closest); - //double m_J = mp->m_J0; - //vector fe(1, 0.0); - //fe[0] = -v_rate * dt * m_J * sqrt(2); - - R.Assemble(lm, fe); } //! evaluate stiffness matrix @@ -312,6 +367,7 @@ void FESolutePointSource::StiffnessMatrix(FELinearSystem& S, const FETimeInfo& t double dt = tp.timeIncrement; // get the domain in which this element resides + m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); FESolidDomain* dom = dynamic_cast(m_el->GetMeshPartition()); FEMesh* mesh = dom->GetMesh(); @@ -328,47 +384,143 @@ void FESolutePointSource::StiffnessMatrix(FELinearSystem& S, const FETimeInfo& t double* w = m_el->GaussWeights(); double Ve = mesh->ElementVolume(*m_el); double v_rate = m_rate / Ve; - for (int k = 0; k < nint; k++) { - FEMaterialPoint* mp = m_el->GetMaterialPoint(k); - double m_J = mp->m_J0; - for (int i = 0; i < neln; ++i) - for (int j = 0; j < neln; ++j) - { - ke[i][j] = H[i] * m_J * H[j] * v_rate * w[k] * dt; - } - } - // get the LM vector - vector lm(neln, -1); - for (int i = 0; i < neln; ++i) lm[i] = mesh->Node(m_el->m_node[i]).m_ID[m_dofC]; - // get the nodes - vector nodes(neln, -1); - for (int i = 0; i < neln; ++i) nodes[i] = m_el->m_node[i]; + double c[FEElement::MAX_NODES]; + for (int i = 0; i < neln; ++i) c[i] = mesh->Node(m_el->m_node[i]).get(m_dofC); - // assemble into global matrix - ke.SetIndices(lm); - ke.SetNodes(nodes); - S.Assemble(ke); -} + std::vector possible_nodes; + double total_elem = 0; + FindNodesInRadius(possible_nodes, total_elem); + if (possible_nodes.size() == 0) + { + for (int k = 0; k < nint; k++) { + FEMaterialPoint* mp = m_el->GetMaterialPoint(k); + double m_J = mp->m_J0; + for (int i = 0; i < neln; ++i) + for (int j = 0; j < neln; ++j) + { + ke[i][j] = H[i] * m_J * H[j] * v_rate * w[k] * dt; + } + } + // get the LM vector + vector lm(neln, -1); + for (int i = 0; i < neln; ++i) lm[i] = mesh->Node(m_el->m_node[i]).m_ID[m_dofC]; -void FESolutePointSource::SetAccumulateFlag(bool b) { - m_accumulate = b; + // get the nodes + vector nodes(neln, -1); + for (int i = 0; i < neln; ++i) nodes[i] = m_el->m_node[i]; + // assemble into global matrix + ke.SetIndices(lm); + ke.SetNodes(nodes); + S.Assemble(ke); + } + else + { + for (auto iter = possible_nodes.begin(); iter != possible_nodes.end(); ++iter) + { + FESolidElement* c_el = dynamic_cast(*iter); + double r[3]; + dom->ProjectToElement(*c_el, m_pos, r); + vec3d r3 = ClampNatC(r); + r[0] = r3.x; r[1] = r3.y; r[2] = r3.z; + double* w = c_el->GaussWeights(); + c_el->shape_fnc(H, r[0], r[1], r[2]); + nint = c_el->GaussPoints(); + neln = c_el->Nodes(); + FEElementMatrix ke(neln, neln); ke.zero(); + for (int k = 0; k < nint; k++) { + FEMaterialPoint* mp = c_el->GetMaterialPoint(k); + double m_J = mp->m_J0; + for (int i = 0; i < neln; ++i) + for (int j = 0; j < neln; ++j) + { + ke[i][j] = H[i] * m_J * H[j] * v_rate * w[k] * dt; + } + } + // get the LM vector + vector lm(neln, -1); + for (int i = 0; i < neln; ++i) lm[i] = mesh->Node(c_el->m_node[i]).m_ID[m_dofC]; + + // get the nodes + vector nodes(neln, -1); + for (int i = 0; i < neln; ++i) nodes[i] = c_el->m_node[i]; + // assemble into global matrix + ke.SetIndices(lm); + ke.SetNodes(nodes); + S.Assemble(ke); + } + } } -void FESolutePointSource::SetAccumulateCAFlag(bool b) { - m_accumulate_ca = b; -} +void FESolutePointSource::FindNodesInRadius(std::vector& possible_nodes, double& total_elem) { -std::vector FESolutePointSource::FindIntInRadius() { - // determine if the radius exceeds the boundaries of the element - int nint = m_el->GaussPoints(); - std::vector possible_ints; - for (int i = 0; i < nint; i++) { - FEMaterialPoint* mp = m_el->GetMaterialPoint(i); - vec3d disp = mp->m_r0 - m_pos; - if (disp.norm() <= m_radius) { - possible_ints.push_back(mp); + // get element and set up buffers + m_el = dynamic_cast(m_search.FindElement(m_pos, m_q)); + std::unordered_set visited; + std::set next; + //std::vector possible_ints; + visited.reserve(1000); + possible_nodes.reserve(500); + + //we will need to check the current element first + next.insert(m_el); + // create the element adjacency list. + FEDomain* dom = dynamic_cast(m_el->GetMeshPartition()); + FEMultiphasic* mat = dynamic_cast(dom->GetMaterial()); + if (mat == nullptr) return; + // calculate the element volume + auto mesh = dom->GetMesh(); + // create the element-element list + FEElemElemList EEL; + EEL.Create(mesh); + + //while there are still elements to evaluate + while (next.size()) { + // get the element to be evaluated + FESolidElement* cur = *next.begin(); + // remove the current element from the next buffer and add it to the visited buffer + next.erase(next.begin()); + visited.insert(cur); + // get the current element bounds + std::vector cur_element_bounds; + // add integration points within the radius + bool int_flag = false; + for (int i = 0; i < cur->Nodes(); i++) + { + FENode* mn = &(mesh->Node(cur->m_node[i])); + vec3d disp = mn->m_rt - m_pos; + if (disp.norm() <= m_radius) { + possible_nodes.push_back(cur); + int_flag = true; + } + } + if (int_flag) + { + total_elem++; + } + + // Add neighboring element to the next buffer as long as they haven't been visited. + // get the global ID of the current element + int cur_id = cur->GetID() - 1; + // for each neighboring element + for (int i = 0; i < EEL.NeighborSize(); i++) + { + if (EEL.Neighbor(cur_id, i)) + { + // if that element has not been visited yet add it to the next list + if (!visited.count(dynamic_cast(EEL.Neighbor(cur_id, i)))) + { + next.insert(dynamic_cast(EEL.Neighbor(cur_id, i))); + } + } } } - return possible_ints; +} + +vec3d FESolutePointSource::ClampNatC(double r[3]) +{ + return vec3d(std::max(std::min(1.0, r[0]), -1.0), + std::max(std::min(1.0, r[1]), -1.0), + std::max(std::min(1.0, r[2]), -1.0) + ); } \ No newline at end of file diff --git a/FEBioMix/FESolutePointSource.h b/FEBioMix/FESolutePointSource.h index 3b4765f17..5463aaaa7 100644 --- a/FEBioMix/FESolutePointSource.h +++ b/FEBioMix/FESolutePointSource.h @@ -26,6 +26,7 @@ SOFTWARE.*/ #pragma once #include #include +#include class FESolidElement; @@ -54,6 +55,14 @@ class FESolutePointSource : public FEBodyLoad double GetRate() const; + double GetdC() const; + + double GetdCp() const; + + void SetdC(double dC); + + void SetdCp(double dCp); + void SetAccumulateFlag(bool b); void SetAccumulateCAFlag(bool b); @@ -65,15 +74,20 @@ class FESolutePointSource : public FEBodyLoad void StiffnessMatrix(FELinearSystem& S, const FETimeInfo& tp) override; //! return all the elements in the given radius - std::vector FindIntInRadius(); + void FindNodesInRadius(std::vector& possible_nodes, double& total_elem); + + vec3d ClampNatC(double r[3]); private: int m_soluteId; //!< solute ID double m_rate; //!< production rate vec3d m_pos; //!< position of source - bool m_accumulate; //!< accumulate flag + bool m_accumulate = false; //!< accumulate flag bool m_accumulate_ca; //! < accumulate actual concentration flag double m_radius; + double m_Vc; + double m_dC = 0.0; + double m_dCp = 0.0; private: FEOctreeSearch m_search; diff --git a/FEBioMix/FESupplyBinding.cpp b/FEBioMix/FESupplyBinding.cpp index 9627deef4..d39afd8ec 100644 --- a/FEBioMix/FESupplyBinding.cpp +++ b/FEBioMix/FESupplyBinding.cpp @@ -80,7 +80,7 @@ double FESupplyBinding::ReceptorLigandSupply(FEMaterialPoint& mp) double J = et.m_J; double ca = spt.m_ca[0]; - double phi0 = ppt.m_phi0; + double phi0 = ppt.m_phi0t; double cr = (J-phi0)*ca; double crc = spt.m_sbmr[0]; double crchat = m_kf*cr*(m_crt-crc) - m_kr*crc; @@ -104,7 +104,7 @@ double FESupplyBinding::ReceptorLigandConcentrationSS(FEMaterialPoint& mp) double J = et.m_J; double ca = spt.m_ca[0]; - double phi0 = ppt.m_phi0; + double phi0 = ppt.m_phi0t; double cr = (J-phi0)*ca; double Kd = m_kr/m_kf; // dissociation constant double crc = m_crt*cr/(Kd+cr); diff --git a/FEBioMix/FESupplyMichaelisMenten.cpp b/FEBioMix/FESupplyMichaelisMenten.cpp index 327dabf60..23010c056 100644 --- a/FEBioMix/FESupplyMichaelisMenten.cpp +++ b/FEBioMix/FESupplyMichaelisMenten.cpp @@ -56,7 +56,7 @@ double FESupplyMichaelisMenten::Supply(FEMaterialPoint& mp) double J = et.m_J; double ca = spt.m_ca[0]; - double phi0 = ppt.m_phi0; + double phi0 = ppt.m_phi0t; double cr = (J-phi0)*ca; double crhat = -m_Vmax*cr/(m_Km+cr); @@ -80,7 +80,7 @@ double FESupplyMichaelisMenten::Tangent_Supply_Concentration(FEMaterialPoint &mp double J = et.m_J; double ca = spt.m_ca[0]; - double phi0 = ppt.m_phi0; + double phi0 = ppt.m_phi0t; double cr = (J-phi0)*ca; double dcrhatdcr = -m_Vmax*m_Km/SQR(m_Km+cr); diff --git a/FEBioMix/FESupplySynthesisBinding.cpp b/FEBioMix/FESupplySynthesisBinding.cpp index 294c52763..93c26782b 100644 --- a/FEBioMix/FESupplySynthesisBinding.cpp +++ b/FEBioMix/FESupplySynthesisBinding.cpp @@ -81,7 +81,7 @@ double FESupplySynthesisBinding::ReceptorLigandSupply(FEMaterialPoint& mp) double J = et.m_J; double ca = spt.m_ca[0]; - double phi0 = ppt.m_phi0; + double phi0 = ppt.m_phi0t; double cr = (J-phi0)*ca; double crc = spt.m_sbmr[0]; double crchat = m_kf*cr*(m_crt-crc) - m_kr*crc; @@ -105,7 +105,7 @@ double FESupplySynthesisBinding::ReceptorLigandConcentrationSS(FEMaterialPoint& double J = et.m_J; double ca = spt.m_ca[0]; - double phi0 = ppt.m_phi0; + double phi0 = ppt.m_phi0t; double cr = (J-phi0)*ca; double Kd = m_kr/m_kf; // dissociation constant double crc = m_crt*cr/(Kd+cr); diff --git a/FEBioMix/FETriphasic.cpp b/FEBioMix/FETriphasic.cpp index 9e36b5102..e2cb966c2 100644 --- a/FEBioMix/FETriphasic.cpp +++ b/FEBioMix/FETriphasic.cpp @@ -151,7 +151,7 @@ double FETriphasic::Porosity(FEMaterialPoint& pt) double J = et.m_J; // porosity // double phiw = 1 - m_phi0/J; - double phi0 = pet.m_phi0; + double phi0 = pet.m_phi0t; double phiw = 1 - phi0/J; // check for pore collapse // TODO: throw an error if pores collapse @@ -169,8 +169,9 @@ double FETriphasic::FixedChargeDensity(FEMaterialPoint& pt) // relative volume double J = et.m_J; - double phi0 = pet.m_phi0; - double cF = m_cFr(pt)*(1-phi0)/(J-phi0); + double phi0 = pet.m_phi0; + double phir = pet.m_phi0t; + double cF = m_cFr(pt)*(1-phi0)/(J-phir); return cF; } @@ -276,7 +277,7 @@ tens4ds FETriphasic::Tangent(FEMaterialPoint& mp) // relative volume and solid volume fraction double J = ept.m_J; - double phi0 = ppt.m_phi0; + double phi0 = ppt.m_phi0t; // get the charge density and its derivatives double cF = FixedChargeDensity(mp); @@ -543,7 +544,7 @@ void FETriphasic::PartitionCoefficientFunctions(FEMaterialPoint& mp, vector()); // initialize referential solid volume fraction - pt.m_phi0 = pmb->m_phi0(mp); + pt.m_phi0t = pmb->m_phi0(mp); // initialize effective fluid pressure, its gradient, and fluid flux pt.m_p = el.Evaluate(p0, n); @@ -221,7 +221,7 @@ void FETriphasicDomain::Reset() FESolutesMaterialPoint& ps = *(mp.ExtractData()); // initialize referential solid volume fraction - pt.m_phi0 = pmb->m_phi0(mp); + pt.m_phi0 = pt.m_phi0t = pmb->m_phi0(mp); // initialize multiphasic solutes ps.m_nsol = nsol; @@ -275,7 +275,7 @@ void FETriphasicDomain::PreSolveUpdate(const FETimeInfo& timeInfo) pe.m_J = defgrad(el, pe.m_F, j); // reset referential solid volume fraction at previous time - pt.m_phi0p = pt.m_phi0; + pt.m_phi0p = pt.m_phi0t; // reset determinant of solid deformation gradient at previous time pt.m_Jp = pe.m_J; diff --git a/FEBioXML/FEBioOutputSection.cpp b/FEBioXML/FEBioOutputSection.cpp index dafccd88b..6a75c2c38 100644 --- a/FEBioXML/FEBioOutputSection.cpp +++ b/FEBioXML/FEBioOutputSection.cpp @@ -226,7 +226,17 @@ void FEBioOutputSection::ParseLogfile(XMLTag &tag) const char* sz = tag.AttributeValue("surface"); FESurface* surf = mesh.FindSurface(sz); - if (surf == nullptr) throw XMLReader::InvalidAttributeValue(tag, "surface", sz); + if (surf == nullptr) + { + FEFacetSet* pfs = mesh.FindFacetSet(sz); + if (pfs == nullptr) throw XMLReader::InvalidAttributeValue(tag, "surface", sz); + + surf = new FESurface(&fem); + surf->Create(*pfs); + surf->SetName(sz); + surf->Init(); + mesh.AddSurface(surf); + } std::vector items; string_to_int_vector(tag.szvalue(), items); diff --git a/FEBioXML/FileImport.cpp b/FEBioXML/FileImport.cpp index 20f669949..a7d66ef32 100644 --- a/FEBioXML/FileImport.cpp +++ b/FEBioXML/FileImport.cpp @@ -266,7 +266,7 @@ int enumValue(const char* val, const char* szenum) // get the string's length. // there could be a comma, so correct for that. - int L = strlen(val); + size_t L = strlen(val); const char* c = strchr(val, ','); if (c) L = c - val; @@ -719,10 +719,6 @@ bool FEFileSection::ReadParameter(XMLTag& tag, FEParameterList& pl, const char* // assign the valuator to the parameter p.setValuator(val); - - // do the initialization. - // TODO: Is this a good place to do this? - if (val->Init() == false) throw XMLReader::InvalidTag(tag); } break; case FE_PARAM_MAT3D_MAPPED: @@ -757,10 +753,6 @@ bool FEFileSection::ReadParameter(XMLTag& tag, FEParameterList& pl, const char* // assign the valuator to the parameter p.setValuator(val); - - // do the initialization. - // TODO: Is this a good place to do this? - if (val->Init() == false) throw XMLReader::InvalidTag(tag); } break; case FE_PARAM_MAT3DS_MAPPED: diff --git a/FECore/Callback.h b/FECore/Callback.h index 159576347..00d5ae122 100644 --- a/FECore/Callback.h +++ b/FECore/Callback.h @@ -51,7 +51,9 @@ class FEModel; #define CB_RESET 0x00000800 //!< Called after FEModel::Reset #define CB_MODEL_UPDATE 0x00001000 //!< Called at the end of FEModel::Update #define CB_TIMESTEP_SOLVED 0x00002000 //!< Called at FEAnalysis::SolveTimeStep after the solver returns. -#define CB_USER1 0x00010000 //!< can be used by users +#define CB_SERIALIZE_SAVE 0x00004000 //!< Called at the end of FEModel::Serialize when saving +#define CB_SERIALIZE_LOAD 0x00008000 //!< Called at the end of FEModel::Serialize when loading +#define CB_USER1 0x01000000 //!< can be used by users typedef unsigned int FECORE_CB_WHEN; typedef bool(*FECORE_CB_FNC)(FEModel*, unsigned int, void*); diff --git a/FECore/FEConstValueVec3.cpp b/FECore/FEConstValueVec3.cpp index a86237dcd..c2f8d90a8 100644 --- a/FECore/FEConstValueVec3.cpp +++ b/FECore/FEConstValueVec3.cpp @@ -24,6 +24,7 @@ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.*/ #include "stdafx.h" +#include "FEModel.h" #include "FEConstValueVec3.h" #include "FEMaterialPoint.h" #include "FEMeshPartition.h" @@ -69,29 +70,34 @@ bool FEMathValueVec3::Init() //--------------------------------------------------------------------------------------- bool FEMathValueVec3::create(const std::string& sx, const std::string& sy, const std::string& sz) { - for (int i = 0; i < 3; ++i) + FECoreBase* pc = nullptr; + if (pc == nullptr) { - m_math[i].AddVariable("X"); - m_math[i].AddVariable("Y"); - m_math[i].AddVariable("Z"); + // try to find the owner of this parameter + // First, we need the model parameter + FEModelParam* param = GetModelParam(); + if (param == nullptr) return false; + + // we'll need the model for this + FEModel* fem = GetFEModel(); + if (fem == nullptr) return false; + + // Now try to find the owner of this parameter + pc = fem->FindParameterOwner(param); } - bool b; - b = m_math[0].Create(sx); assert(b); - b = m_math[1].Create(sy); assert(b); - b = m_math[2].Create(sz); assert(b); + + if (m_math[0].Init(sx, pc) == false) return false; + if (m_math[1].Init(sy, pc) == false) return false; + if (m_math[2].Init(sz, pc) == false) return false; return true; } vec3d FEMathValueVec3::operator()(const FEMaterialPoint& pt) { - std::vector var(3); - var[0] = pt.m_r0.x; - var[1] = pt.m_r0.y; - var[2] = pt.m_r0.z; - double vx = m_math[0].value_s(var); - double vy = m_math[1].value_s(var); - double vz = m_math[2].value_s(var); + double vx = m_math[0].value(GetFEModel(), pt); + double vy = m_math[1].value(GetFEModel(), pt); + double vz = m_math[2].value(GetFEModel(), pt); return vec3d(vx, vy, vz); } diff --git a/FECore/FEConstValueVec3.h b/FECore/FEConstValueVec3.h index ad55009ed..83194bb3d 100644 --- a/FECore/FEConstValueVec3.h +++ b/FECore/FEConstValueVec3.h @@ -69,7 +69,7 @@ class FECORE_API FEMathValueVec3 : public FEVec3dValuator private: std::string m_expr; - MSimpleExpression m_math[3]; + FEMathExpression m_math[3]; DECLARE_FECORE_CLASS(); }; diff --git a/FECore/FECore.cpp b/FECore/FECore.cpp index e9009f151..57172305d 100644 --- a/FECore/FECore.cpp +++ b/FECore/FECore.cpp @@ -168,4 +168,5 @@ REGISTER_FECORE_CLASS(FELogDomainVolume, "volume"); REGISTER_FECORE_CLASS(FELogAvgDomainData, "avg"); REGISTER_FECORE_CLASS(FELogPctDomainData, "pct"); REGISTER_FECORE_CLASS(FELogSolutionNorm, "solution_norm"); +REGISTER_FECORE_CLASS(FELogFaceArea , "facet area"); } diff --git a/FECore/FECoreBase.cpp b/FECore/FECoreBase.cpp index a270159ba..4a7583264 100644 --- a/FECore/FECoreBase.cpp +++ b/FECore/FECoreBase.cpp @@ -300,7 +300,32 @@ bool FECoreBase::Init() } } } + else if (pi.type() == FE_PARAM_VEC3D_MAPPED) + { + for (int j = 0; j < pi.dim(); ++j) + { + FEParamVec3& pd = pi.value(j); + if (pd.Init() == false) + { + feLogError("Failed to initialize parameter %s", pi.name()); + return false; + } + } + } + else if (pi.type() == FE_PARAM_MAT3D_MAPPED) + { + for (int j = 0; j < pi.dim(); ++j) + { + FEParamMat3d& pd = pi.value(j); + if (pd.Init() == false) + { + feLogError("Failed to initialize parameter %s", pi.name()); + return false; + } + } + } } + // check the parameter ranges if (Validate() == false) return false; diff --git a/FECore/FECorePlot.cpp b/FECore/FECorePlot.cpp index 6c888e2ea..f1f0ae6ab 100644 --- a/FECore/FECorePlot.cpp +++ b/FECore/FECorePlot.cpp @@ -32,6 +32,7 @@ SOFTWARE.*/ #include "FESolidDomain.h" #include "FEModelParam.h" #include "FEBodyLoad.h" +#include "FENodalLoad.h" #include "FEPlotData.h" #include "FESurface.h" #include "FEMaterialPointProperty.h" @@ -83,8 +84,9 @@ bool FEPlotParameter::SetFilter(const char* sz) if (itemList == 0) { // for some classes, the item list can be empty - if (dynamic_cast(pc)) { SetRegionType(FE_REGION_DOMAIN); m_dom = &(dynamic_cast(pc))->GetDomainList(); } + if (dynamic_cast(pc)) { SetRegionType(FE_REGION_DOMAIN); m_dom = &(dynamic_cast(pc))->GetDomainList(); } // else if (dynamic_cast(pc)) { SetRegionType(FE_REGION_SURFACE); m_surf = dynamic_cast(pc)->GetSurface().GetFacetSet(); } + else if (dynamic_cast(pc)) SetRegionType(FE_REGION_NODE); else return false; } else @@ -439,10 +441,18 @@ bool FEPlotParameter::Save(FEMesh& mesh, FEDataStream& a) { FEParamDouble& map = m_param.value(); FENodeSet* nset = dynamic_cast(map.GetItemList()); - if (nset == 0) return false; - - // write the nodal values - writeNodalValues(*nset, a, map); + if (nset) + writeNodalValues(*nset, a, map); + else + { + writeNodalValues(mesh, a, [&](const FENode& node) { + FEMaterialPoint mp; + mp.m_r0 = node.m_r0; + mp.m_index = -1; + double v = map(mp); + return v; + }); + } return true; } diff --git a/FECore/FEDataMap.cpp b/FECore/FEDataMap.cpp index c0414b08c..937d670fc 100644 --- a/FECore/FEDataMap.cpp +++ b/FECore/FEDataMap.cpp @@ -29,6 +29,7 @@ SOFTWARE.*/ #include "stdafx.h" #include "FEDataMap.h" #include "fecore_type.h" +#include "DumpStream.h" //----------------------------------------------------------------------------- FEDataMap::FEDataMap(FEDataMapType mapType, FEDataType dataType) : FEDataArray(mapType, dataType) {} @@ -44,3 +45,13 @@ void FEDataMap::SetName(const std::string& name) //----------------------------------------------------------------------------- const std::string& FEDataMap::GetName() const { return m_name; } + +//----------------------------------------------------------------------------- +void FEDataMap::Serialize(DumpStream& ar) +{ + FEDataArray::Serialize(ar); + if (ar.IsShallow() == false) + { + ar & m_name; + } +} diff --git a/FECore/FEDataMap.h b/FECore/FEDataMap.h index 3647eeeee..c4a994ba2 100644 --- a/FECore/FEDataMap.h +++ b/FECore/FEDataMap.h @@ -59,6 +59,9 @@ class FECORE_API FEDataMap : public FEDataArray // return the item list associated with this map virtual FEItemList* GetItemList() = 0; +public: + void Serialize(DumpStream& ar) override; + protected: std::string m_name; // name of data map TODO: Move to base class? }; diff --git a/FECore/FEElemElemList.h b/FECore/FEElemElemList.h index 39f13df34..4672f86f2 100644 --- a/FECore/FEElemElemList.h +++ b/FECore/FEElemElemList.h @@ -59,6 +59,9 @@ class FECORE_API FEElemElemList //! Find the j-th neighbor element of element n int NeighborIndex(int n, int j) { return m_peli[m_ref[n] + j]; } + //! Return the size of the neighbor vector + int NeighborSize() { return m_pel.size()/m_ref.size(); } + protected: //! Initialization void Init(); diff --git a/FECore/FELogElementVolume.cpp b/FECore/FELogElementVolume.cpp index 05a09e456..74035d76f 100644 --- a/FECore/FELogElementVolume.cpp +++ b/FECore/FELogElementVolume.cpp @@ -25,6 +25,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.*/ #include "FELogElementVolume.h" #include "FEModel.h" +#include "FESurface.h" FELogElementVolume::FELogElementVolume(FEModel* fem) : FELogElemData(fem) { @@ -36,3 +37,12 @@ double FELogElementVolume::value(FEElement& el) FEMesh& mesh = GetFEModel()->GetMesh(); return mesh.CurrentElementVolume(el); } + +FELogFaceArea::FELogFaceArea(FEModel* fem) : FEFaceLogData(fem) {} + +double FELogFaceArea::value(FESurfaceElement& el) +{ + FESurface* surface = dynamic_cast(el.GetMeshPartition()); + if (surface == nullptr) return 0.0; + return surface->CurrentFaceArea(el); +} diff --git a/FECore/FELogElementVolume.h b/FECore/FELogElementVolume.h index 973f19574..5c8619616 100644 --- a/FECore/FELogElementVolume.h +++ b/FECore/FELogElementVolume.h @@ -25,6 +25,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.*/ #pragma once #include "ElementDataRecord.h" +#include "FaceDataRecord.h" class FELogElementVolume : public FELogElemData { @@ -32,3 +33,10 @@ class FELogElementVolume : public FELogElemData FELogElementVolume(FEModel* fem); double value(FEElement& el) override; }; + +class FELogFaceArea : public FEFaceLogData +{ +public: + FELogFaceArea(FEModel* fem); + double value(FESurfaceElement& el) override; +}; diff --git a/FECore/FEModel.cpp b/FECore/FEModel.cpp index 053bee4a1..bd8cbb385 100644 --- a/FECore/FEModel.cpp +++ b/FECore/FEModel.cpp @@ -1161,6 +1161,14 @@ bool FEModel::RCI_Init() return true; } +bool FEModel::RCI_Restart() +{ + FEAnalysis* step = GetCurrentStep(); + if (step == nullptr) return false; + + return step->InitSolver(); +} + bool FEModel::RCI_Advance() { // get the current step @@ -1646,38 +1654,44 @@ bool FEModel::EvaluateLoadParameters() { double s = GetLoadController(nlc)->Value(); FEParam* p = pi.param; + FECoreBase* parent = dynamic_cast(p->parent()); + if (m_imp->m_printParams) + { + if (parent && (parent->GetName().empty() == false)) + { + const char* pname = parent->GetName().c_str(); + feLog("Setting parameter \"%s.%s\" to : ", pname, p->name()); + } + else + feLog("Setting parameter \"%s\" to : ", p->name()); + }; switch (p->type()) { case FE_PARAM_INT: { p->value() = (int)s; - if (m_imp->m_printParams) - feLog("Setting parameter \"%s\" to : %d\n", p->name(), p->value()); + if (m_imp->m_printParams) feLog("%d\n", p->value()); } break; case FE_PARAM_DOUBLE: { p->value() = pi.m_scl*s; - if (m_imp->m_printParams) - feLog("Setting parameter \"%s\" to : %lg\n", p->name(), p->value()); + if (m_imp->m_printParams) feLog("%lg\n", p->value()); } break; case FE_PARAM_BOOL: { p->value() = (s > 0 ? true : false); - if (m_imp->m_printParams) - feLog("Setting parameter \"%s\" to : %s\n", p->name(), (p->value() ? "true" : "false")); + if (m_imp->m_printParams) feLog("%s\n", (p->value() ? "true" : "false")); } break; case FE_PARAM_VEC3D: { vec3d& v = p->value(); p->value() = pi.m_vscl*s; - if (m_imp->m_printParams) - feLog("Setting parameter \"%s\" to : %lg, %lg, %lg\n", p->name(), v.x, v.y, v.z); + if (m_imp->m_printParams) feLog("%lg, %lg, %lg\n", v.x, v.y, v.z); } break; case FE_PARAM_DOUBLE_MAPPED: { p->value().SetScaleFactor(s * pi.m_scl); - if (m_imp->m_printParams) - feLog("Setting parameter \"%s\" to : %lg\n", p->name(), p->value().GetScaleFactor()); + if (m_imp->m_printParams) feLog("%lg\n", p->value().GetScaleFactor()); } break; case FE_PARAM_VEC3D_MAPPED : p->value().SetScaleFactor(s* pi.m_scl); break; @@ -1733,6 +1747,10 @@ bool FEModel::DoCallback(unsigned int nevent) { throw; } + catch (DoRunningRestart) + { + throw; + } catch (std::exception e) { throw; @@ -2202,7 +2220,7 @@ void FEModel::Serialize(DumpStream& ar) TRACK_TIME(TimerID::Timer_Update); m_imp->Serialize(ar); - if (ar.IsShallow()) return; + DoCallback(ar.IsSaving() ? CB_SERIALIZE_SAVE : CB_SERIALIZE_LOAD); } //----------------------------------------------------------------------------- diff --git a/FECore/FEModel.h b/FECore/FEModel.h index 7f9d9a877..e54ea89db 100644 --- a/FECore/FEModel.h +++ b/FECore/FEModel.h @@ -139,6 +139,7 @@ class FECORE_API FEModel : public FECoreBase, public CallbackHandler public: // reverse control solver interface bool RCI_Init(); + bool RCI_Restart(); bool RCI_Rewind(); bool RCI_Advance(); bool RCI_Finish(); diff --git a/FECore/FEModelParam.cpp b/FECore/FEModelParam.cpp index 6ba01d4de..eddb11ff9 100644 --- a/FECore/FEModelParam.cpp +++ b/FECore/FEModelParam.cpp @@ -48,6 +48,10 @@ FEModelParam::~FEModelParam() void FEModelParam::Serialize(DumpStream& ar) { ar & m_scl; + if (ar.IsShallow() == false) + { + ar & m_dom; + } } //--------------------------------------------------------------------------------------- @@ -134,6 +138,11 @@ FEParamVec3::FEParamVec3(const FEParamVec3& p) m_dom = p.m_dom; } +bool FEParamVec3::Init() +{ + return (m_val ? m_val->Init() : true); +} + // set the value void FEParamVec3::operator = (const vec3d& v) { @@ -197,6 +206,11 @@ void FEParamMat3d::operator = (const FEParamMat3d& p) // m_dom = p.m_dom; } +bool FEParamMat3d::Init() +{ + return (m_val ? m_val->Init() : true); +} + // set the valuator void FEParamMat3d::setValuator(FEMat3dValuator* val) { diff --git a/FECore/FEModelParam.h b/FECore/FEModelParam.h index 2cd7f2627..80f2d97f2 100644 --- a/FECore/FEModelParam.h +++ b/FECore/FEModelParam.h @@ -109,6 +109,8 @@ class FECORE_API FEParamVec3 : public FEModelParam FEParamVec3(const FEParamVec3& p); + bool Init(); + // set the value void operator = (const vec3d& v); void operator = (const FEParamVec3& p); @@ -149,6 +151,8 @@ class FECORE_API FEParamMat3d : public FEModelParam void operator = (const mat3d& v); void operator = (const FEParamMat3d& v); + bool Init(); + // set the valuator void setValuator(FEMat3dValuator* val); diff --git a/FECore/FENodeDataMap.cpp b/FECore/FENodeDataMap.cpp index c14aa91c7..6df1f56fa 100644 --- a/FECore/FENodeDataMap.cpp +++ b/FECore/FENodeDataMap.cpp @@ -135,3 +135,23 @@ FEItemList* FENodeDataMap::GetItemList() { return const_cast(m_nodeSet); } + +void FENodeDataMap::Serialize(DumpStream& ar) +{ + FEDataMap::Serialize(ar); + if (ar.IsShallow() == false) + { + if (ar.IsSaving()) + { + // We have to cast the const away before serializing + FENodeSet* ns = const_cast(m_nodeSet); + ar << ns; + } + else + { + FENodeSet* ns; + ar >> ns; + m_nodeSet = ns; + } + } +} diff --git a/FECore/FENodeDataMap.h b/FECore/FENodeDataMap.h index db42fc893..5e21ebdaf 100644 --- a/FECore/FENodeDataMap.h +++ b/FECore/FENodeDataMap.h @@ -44,6 +44,8 @@ class FECORE_API FENodeDataMap : public FEDataMap // return the item list associated with this map FEItemList* GetItemList() override; + void Serialize(DumpStream& ar) override; + public: void setValue(int n, double v) override; void setValue(int n, const vec2d& v) override; diff --git a/FECore/FENormalProjection.cpp b/FECore/FENormalProjection.cpp index ecce2e712..21ce5c44f 100644 --- a/FECore/FENormalProjection.cpp +++ b/FECore/FENormalProjection.cpp @@ -54,7 +54,7 @@ FESurfaceElement* FENormalProjection::Project(vec3d r, vec3d n, double rs[2]) { // let's find all the candidate surface elements setselist; - m_OT.FindCandidateSurfaceElements(r, n, selist); + m_OT.FindCandidateSurfaceElements(r, n, selist, m_rad); // now that we found candidate surface elements, lets see if we can find // those that intersect the ray, then pick the closest intersection @@ -98,7 +98,7 @@ FESurfaceElement* FENormalProjection::Project2(vec3d r, vec3d n, double rs[2]) { // let's find all the candidate surface elements setselist; - m_OT.FindCandidateSurfaceElements(r, n, selist); + m_OT.FindCandidateSurfaceElements(r, n, selist, m_rad); // now that we found candidate surface elements, lets see if we can find // those that intersect the ray, then pick the closest intersection @@ -141,7 +141,7 @@ FESurfaceElement* FENormalProjection::Project3(const vec3d& r, const vec3d& n, d { // let's find all the candidate surface elements setselist; - m_OT.FindCandidateSurfaceElements(r, n, selist); + m_OT.FindCandidateSurfaceElements(r, n, selist, m_rad); double g, gmax = -1e99, r2[2] = {rs[0], rs[1]}; int imin = -1; diff --git a/FECore/FEOctree.cpp b/FECore/FEOctree.cpp index d27b3a571..ad085eee9 100644 --- a/FECore/FEOctree.cpp +++ b/FECore/FEOctree.cpp @@ -188,15 +188,19 @@ bool OTnode::RayIntersectsNode(vec3d p, vec3d n) //----------------------------------------------------------------------------- // Find intersected octree leaves and return a set of their surface elements - -void OTnode::FindIntersectedLeaves(vec3d p, vec3d n, set& sel) +void OTnode::FindIntersectedLeaves(vec3d p, vec3d n, set& sel, double srad) { - if (RayIntersectsNode(p, n)) { + // Check if octree node is within search radius from p. + bool bNodeWithinSRad = ( (cmin.x - srad <= p.x) && (cmax.x + srad >= p.x) && + (cmin.y - srad <= p.y) && (cmax.y + srad >= p.y) && + (cmin.z - srad <= p.z) && (cmax.z + srad >= p.z) ); + + if (bNodeWithinSRad && RayIntersectsNode(p, n)) { int nc = (int)children.size(); // if this node has children, search them for intersections if (nc) { for (int ic=0; ic& sel) +void FEOctree::FindCandidateSurfaceElements(vec3d p, vec3d n, set& sel, double srad) { - root.FindIntersectedLeaves(p, n, sel); + root.FindIntersectedLeaves(p, n, sel, srad); } diff --git a/FECore/FEOctree.h b/FECore/FEOctree.h index 6ee2cf30f..ba424ac4e 100644 --- a/FECore/FEOctree.h +++ b/FECore/FEOctree.h @@ -48,7 +48,7 @@ class OTnode bool ElementIntersectsNode(const int j); void PrintNodeContent(); bool RayIntersectsNode(vec3d p, vec3d n); - void FindIntersectedLeaves(vec3d p, vec3d n, std::set& sel); + void FindIntersectedLeaves(vec3d p, vec3d n, std::set& sel, double srad); void CountNodes(int& nnode, int& nlevel); public: @@ -76,7 +76,7 @@ class FECORE_API FEOctree void Init(const double stol); //! find all candidate surface elements intersected by ray - void FindCandidateSurfaceElements(vec3d p, vec3d n, std::set& sel); + void FindCandidateSurfaceElements(vec3d p, vec3d n, std::set& sel, double srad); protected: FESurface* m_ps; //!< the surface to search diff --git a/FECore/FEScalarValuator.cpp b/FECore/FEScalarValuator.cpp index c046493c0..050631f95 100644 --- a/FECore/FEScalarValuator.cpp +++ b/FECore/FEScalarValuator.cpp @@ -34,69 +34,25 @@ SOFTWARE.*/ #include "DumpStream.h" #include "log.h" -REGISTER_SUPER_CLASS(FEScalarValuator, FESCALARGENERATOR_ID); - -//============================================================================= -BEGIN_FECORE_CLASS(FEConstValue, FEScalarValuator) - ADD_PARAMETER(m_val, "const"); -END_FECORE_CLASS(); - //============================================================================= - -BEGIN_FECORE_CLASS(FEMathValue, FEScalarValuator) - ADD_PARAMETER(m_expr, "math"); -END_FECORE_CLASS(); - -void FEMathValue::setMathString(const std::string& s) -{ - m_expr = s; -} - -bool FEMathValue::Init() -{ - return create(); -} - -void FEMathValue::Serialize(DumpStream& ar) -{ - FEScalarValuator::Serialize(ar); - if (ar.IsShallow()) return; - if (ar.IsLoading()) create(); -} - -bool FEMathValue::create(FECoreBase* pc) +bool FEMathExpression::Init(const std::string& expr, FECoreBase* pc) { - // see if this is already initialized - if (m_math.Variables() > 0) return true; + Clear(); - m_math.AddVariable("X"); - m_math.AddVariable("Y"); - m_math.AddVariable("Z"); - m_math.AddVariable("t"); - bool b = m_math.Create(m_expr, true); + AddVariable("X"); + AddVariable("Y"); + AddVariable("Z"); + AddVariable("t"); + bool b = Create(expr, true); // lookup all the other variables. - if (m_math.Variables() > 4) + if (Variables() > 4) { - if (pc == nullptr) - { - // try to find the owner of this parameter - // First, we need the model parameter - FEModelParam* param = GetModelParam(); - if (param == nullptr) return false; - - // we'll need the model for this - FEModel* fem = GetFEModel(); - if (fem == nullptr) return false; - - // Now try to find the owner of this parameter - pc = fem->FindParameterOwner(param); - if (pc == nullptr) return false; - } + if (pc == nullptr) return false; - for (int i = 4; i < m_math.Variables(); ++i) + for (int i = 4; i < Variables(); ++i) { - MVariable* vari = m_math.Variable(i); + MVariable* vari = Variable(i); ParamString ps(vari->Name().c_str()); FEParam* p = pc->FindParameter(ps); @@ -112,7 +68,7 @@ bool FEMathValue::create(FECoreBase* pc) else { // see if it's a data map - FEModel* fem = GetFEModel(); + FEModel* fem = pc->GetFEModel(); // see if it's a global variable p = fem->FindParameter(ps); @@ -129,14 +85,14 @@ bool FEMathValue::create(FECoreBase* pc) FEDataMap* map = mesh.FindDataMap(vari->Name()); assert(map); - if (map == nullptr) { + if (map == nullptr) { const char* szvar = vari->Name().c_str(); - feLogError("Don't understand variable name \"%s\" in math expression.", szvar); + feLogErrorEx(fem, "Don't understand variable name \"%s\" in math expression.", szvar); return false; } if (map->DataType() != FEDataType::FE_DOUBLE) { const char* szvar = vari->Name().c_str(); - feLogError("Variable \"%s\" is not a scalar variable.", szvar); + feLogErrorEx(fem, "Variable \"%s\" is not a scalar variable.", szvar); return false; } @@ -153,26 +109,19 @@ bool FEMathValue::create(FECoreBase* pc) return b; } -FEMathValue::~FEMathValue() -{ -} - -FEScalarValuator* FEMathValue::copy() +void FEMathExpression::operator = (const FEMathExpression& me) { - FEMathValue* newExpr = new FEMathValue(GetFEModel()); - newExpr->m_expr = m_expr; - newExpr->m_math = m_math; - newExpr->m_vars = m_vars; - return newExpr; + MSimpleExpression::operator=(me); + m_vars = me.m_vars; } -double FEMathValue::operator()(const FEMaterialPoint& pt) +double FEMathExpression::value(FEModel* fem, const FEMaterialPoint& pt) { std::vector var(4 + m_vars.size()); var[0] = pt.m_r0.x; var[1] = pt.m_r0.y; var[2] = pt.m_r0.z; - var[3] = GetFEModel()->GetTime().currentTime; + var[3] = fem->GetTime().currentTime; if (m_vars.empty() == false) { for (int i = 0; i < (int)m_vars.size(); ++i) @@ -191,11 +140,84 @@ double FEMathValue::operator()(const FEMaterialPoint& pt) else { FEDataMap& map = *mp.map; - var[4+i] = map.value(pt); + var[4 + i] = map.value(pt); } } } - return m_math.value_s(var); + return value_s(var); +} + +//============================================================================= +REGISTER_SUPER_CLASS(FEScalarValuator, FESCALARGENERATOR_ID); + +//============================================================================= +BEGIN_FECORE_CLASS(FEConstValue, FEScalarValuator) + ADD_PARAMETER(m_val, "const"); +END_FECORE_CLASS(); + +//============================================================================= + +BEGIN_FECORE_CLASS(FEMathValue, FEScalarValuator) + ADD_PARAMETER(m_expr, "math"); +END_FECORE_CLASS(); + +void FEMathValue::setMathString(const std::string& s) +{ + m_expr = s; +} + +bool FEMathValue::Init() +{ + return create(); +} + +void FEMathValue::Serialize(DumpStream& ar) +{ + FEScalarValuator::Serialize(ar); + if (ar.IsShallow()) return; + if (ar.IsLoading()) create(); +} + +bool FEMathValue::create(FECoreBase* pc) +{ + // see if this is already initialized + if (m_math.Variables() > 0) return true; + + if (pc == nullptr) + { + // try to find the owner of this parameter + // First, we need the model parameter + FEModelParam* param = GetModelParam(); + if (param == nullptr) return false; + + // we'll need the model for this + FEModel* fem = GetFEModel(); + if (fem == nullptr) return false; + + // Now try to find the owner of this parameter + pc = fem->FindParameterOwner(param); + } + + // initialize the math expression + bool b = m_math.Init(m_expr, pc); + return b; +} + +FEMathValue::~FEMathValue() +{ +} + +FEScalarValuator* FEMathValue::copy() +{ + FEMathValue* newExpr = new FEMathValue(GetFEModel()); + newExpr->m_expr = m_expr; + newExpr->m_math = m_math; + return newExpr; +} + +double FEMathValue::operator()(const FEMaterialPoint& pt) +{ + return m_math.value(GetFEModel(), pt); } //--------------------------------------------------------------------------------------- diff --git a/FECore/FEScalarValuator.h b/FECore/FEScalarValuator.h index 9d57834a4..e1c69038c 100644 --- a/FECore/FEScalarValuator.h +++ b/FECore/FEScalarValuator.h @@ -75,15 +75,29 @@ class FECORE_API FEConstValue : public FEScalarValuator }; //--------------------------------------------------------------------------------------- -class FECORE_API FEMathValue : public FEScalarValuator +class FEMathExpression : public MSimpleExpression { struct MathParam { int type; // 0 = param, 1 = map - FEParam* pp; - FEDataMap* map; + FEParam* pp; + FEDataMap* map; }; +public: + bool Init(const std::string& expr, FECoreBase* pc = nullptr); + + void operator = (const FEMathExpression& me); + + double value(FEModel* fem, const FEMaterialPoint& pt); + +private: + std::vector m_vars; +}; + +//--------------------------------------------------------------------------------------- +class FECORE_API FEMathValue : public FEScalarValuator +{ public: FEMathValue(FEModel* fem) : FEScalarValuator(fem) {} ~FEMathValue(); @@ -101,8 +115,7 @@ class FECORE_API FEMathValue : public FEScalarValuator private: std::string m_expr; - MSimpleExpression m_math; - std::vector m_vars; + FEMathExpression m_math; DECLARE_FECORE_CLASS(); }; diff --git a/FECore/FESurface.cpp b/FECore/FESurface.cpp index 20d2eae26..09098decf 100644 --- a/FECore/FESurface.cpp +++ b/FECore/FESurface.cpp @@ -77,7 +77,6 @@ void FESurface::Create(const FEFacetSet& set) assert(m_surf == &set); FEMesh& m = *GetMesh(); - int NN = m.Nodes(); // count nr of faces int faces = set.Faces(); @@ -118,7 +117,6 @@ void FESurface::CreateMaterialPointData() { FESurfaceElement& el = m_el[i]; int nint = el.GaussPoints(); - int neln = el.Nodes(); el.ClearData(); for (int n = 0; n < nint; ++n) { @@ -342,7 +340,6 @@ FEElement* FESurface::FindElement(FESurfaceElement& el) FEMesh& mesh = *GetMesh(); FENodeElemList& NEL = mesh.NodeElementList(); - vector& sf = el.m_node; int node = el.m_node[0]; int nval = NEL.Valence(node); FEElement** ppe = NEL.ElementList(node); @@ -699,6 +696,33 @@ vec3d FESurface::Position(FESurfaceElement& el, double r, double s) return q; } +//----------------------------------------------------------------------------- +//! This function calculates the position of integration point n + +vec3d FESurface::Position(FESurfaceElement &el, int n) +{ + // get the mesh to which this surface belongs + FEMesh& mesh = *m_pMesh; + + // number of element nodes + int ne = el.Nodes(); + + // get the elements nodal positions + vec3d y[FEElement::MAX_NODES]; + if (!m_bshellb) for (int i = 0; i(*el.GetMaterialPoint(n)); - double* N = el.H(n); double* Gr = el.Gr(n); double* Gs = el.Gs(n); diff --git a/FECore/FESurface.h b/FECore/FESurface.h index 4c6bf94bd..712175674 100644 --- a/FECore/FESurface.h +++ b/FECore/FESurface.h @@ -166,6 +166,9 @@ class FECORE_API FESurface : public FEMeshPartition //! Get the spatial position given natural coordinates vec3d Position(FESurfaceElement& el, double r, double s); + //! Get the spatial position of an integration point + vec3d Position(FESurfaceElement& el, int n); + //! Get the nodal coordinates of an element void NodalCoordinates(FESurfaceElement& el, vec3d* re); diff --git a/FECore/FESurfaceLoad.cpp b/FECore/FESurfaceLoad.cpp index ef7a2f6eb..23992c47d 100644 --- a/FECore/FESurfaceLoad.cpp +++ b/FECore/FESurfaceLoad.cpp @@ -71,5 +71,12 @@ void FESurfaceLoad::Serialize(DumpStream& ar) if (ar.IsShallow()) return; ar & m_dof; - ar & m_psurf; + + // the mesh manages surfaces for surface loads + if (m_psurf && ar.IsLoading()) + { + FEMesh* pm = m_psurf->GetMesh(); + assert(pm->FindSurface(m_psurf->GetName()) == nullptr); + pm->AddSurface(m_psurf); + } } diff --git a/FECore/FESurfaceMap.cpp b/FECore/FESurfaceMap.cpp index 6cd530c23..656ab0a12 100644 --- a/FECore/FESurfaceMap.cpp +++ b/FECore/FESurfaceMap.cpp @@ -169,6 +169,17 @@ void FESurfaceMap::Serialize(DumpStream& ar) if (ar.IsShallow()) return; ar & m_maxFaceNodes; ar & m_name; + if (ar.IsSaving()) + { + FEFacetSet* fs = const_cast(m_surf); + ar << fs; + } + else + { + FEFacetSet* fs = nullptr; + ar >> fs; + m_surf = fs; + } } //----------------------------------------------------------------------------- diff --git a/FECore/fecore_debug.cpp b/FECore/fecore_debug.cpp index e0b42f4fe..377e3769c 100644 --- a/FECore/fecore_debug.cpp +++ b/FECore/fecore_debug.cpp @@ -30,6 +30,8 @@ SOFTWARE.*/ #include "fecore_debug.h" #include #include +#include +#include "version.h" using namespace std; std::list FECoreDebugger::m_var; @@ -113,6 +115,28 @@ void FECoreDebugger::Remove(FECoreDebugger::Variable* pvar) } } +FILE* FECoreDebugger::m_fp = nullptr; + +void FECoreDebugger::Print(const char* szformat, ...) +{ + if (m_fp == nullptr) + { + char fileName[256] = { 0 }; + sprintf(fileName, "febio_%d.%d.%d_debug.log", FE_SDK_MAJOR_VERSION, FE_SDK_SUB_VERSION, FE_SDK_SUBSUB_VERSION); + m_fp = fopen(fileName, "wt"); assert(m_fp); + } + + if (m_fp) + { + va_list args; + va_start(args, szformat); + vfprintf(m_fp, szformat, args); + va_end(args); + + fflush(m_fp); + } +} + template <> void fecore_print_T(matrix* pd) { matrix& m = *pd; diff --git a/FECore/fecore_debug_t.h b/FECore/fecore_debug_t.h index f38ac6556..f2ce3a77d 100644 --- a/FECore/fecore_debug_t.h +++ b/FECore/fecore_debug_t.h @@ -31,6 +31,7 @@ SOFTWARE.*/ #include #include #include +#include #include "matrix.h" #include "mat3d.h" #include "vec3d.h" @@ -40,7 +41,7 @@ SOFTWARE.*/ // Don't use anything in here or include this file directly. // Instead include fecore_debug.h and use the user functions and macros defined there. -template void fecore_print_T(T* pd) { cout << (*pd); } +template void fecore_print_T(T* pd) { std::cout << (*pd); } template <> void fecore_print_T(matrix* pd); template <> void fecore_print_T(mat3d* pd); template <> void fecore_print_T(mat3ds* pd); @@ -60,13 +61,13 @@ class FECoreDebugger { public: Variable(void* pd, const char* sz) { m_pd = pd; m_szname = sz; } - virtual ~Variable(){} + virtual ~Variable() {} virtual void print() = 0; public: - void* m_pd; - const char* m_szname; + void* m_pd; + const char* m_szname; }; template class Variable_T : public Variable @@ -74,10 +75,10 @@ class FECoreDebugger public: Variable_T(void* pd, const char* sz) : Variable(pd, sz) {} - void print() - { - cout << typeid(T).name() << endl; - fecore_print_T((T*)m_pd); + void print() + { + std::cout << typeid(T).name() << std::endl; + fecore_print_T((T*)m_pd); } }; @@ -87,10 +88,13 @@ class FECoreDebugger static void Add(Variable* pvar); static void Remove(Variable* pvar); + static void Print(const char* szformat, ...); + private: - FECoreDebugger(){} + FECoreDebugger() {} static std::list m_var; + static FILE* m_fp; }; class FECoreWatchVariable diff --git a/FECore/version.h b/FECore/version.h index 48d27b01d..2bdcb0546 100644 --- a/FECore/version.h +++ b/FECore/version.h @@ -34,7 +34,7 @@ SOFTWARE.*/ // of FEBio. In that case, plugins need to be recompiled to be usable with the // newer version of FEBio. #define FE_SDK_MAJOR_VERSION 3 -#define FE_SDK_SUB_VERSION 7 +#define FE_SDK_SUB_VERSION 8 #define FE_SDK_SUBSUB_VERSION 0 //----------------------------------------------------------------------------- diff --git a/NumCore/AccelerateSparseSolver.cpp b/NumCore/AccelerateSparseSolver.cpp index abadb0b63..0f66afc53 100644 --- a/NumCore/AccelerateSparseSolver.cpp +++ b/NumCore/AccelerateSparseSolver.cpp @@ -33,7 +33,7 @@ #include "MatrixTools.h" #include -#ifdef __APPLE__ +#ifdef HAS_ACCEL ////////////////////////////////////////////////////////////// // AccelerateSparseSolver diff --git a/NumCore/BoomerAMGSolver.cpp b/NumCore/BoomerAMGSolver.cpp index ea4e497b0..dd853483a 100644 --- a/NumCore/BoomerAMGSolver.cpp +++ b/NumCore/BoomerAMGSolver.cpp @@ -71,6 +71,7 @@ class BoomerAMGSolver::Implementation int m_interpType; int m_PMaxElmts; int m_NumSweeps; + int m_AggInterpType; int m_AggNumLevels; double m_strong_threshold; int m_nodal; @@ -93,6 +94,7 @@ class BoomerAMGSolver::Implementation m_strong_threshold = 0.5; m_PMaxElmts = 4; m_NumSweeps = 1; + m_AggInterpType = 0; m_AggNumLevels = 0; m_nodal = 0; m_jacobi_pc = false; @@ -239,6 +241,7 @@ class BoomerAMGSolver::Implementation HYPRE_BoomerAMGSetNumSweeps(m_solver, m_NumSweeps); /* Sweeps on each level */ HYPRE_BoomerAMGSetMaxLevels(m_solver, m_maxLevels); /* maximum number of levels */ HYPRE_BoomerAMGSetStrongThreshold(m_solver, m_strong_threshold); + HYPRE_BoomerAMGSetAggInterpType(m_solver, m_AggInterpType); HYPRE_BoomerAMGSetAggNumLevels(m_solver, m_AggNumLevels); HYPRE_BoomerAMGSetMaxIter(m_solver, m_maxIter); HYPRE_BoomerAMGSetTol(m_solver, m_tol); /* conv. tolerance */ @@ -324,6 +327,7 @@ BEGIN_FECORE_CLASS(BoomerAMGSolver, LinearSolver) ADD_PARAMETER(imp->m_strong_threshold, "strong_threshold"); ADD_PARAMETER(imp->m_PMaxElmts , "p_max_elmts"); ADD_PARAMETER(imp->m_NumSweeps , "num_sweeps"); + ADD_PARAMETER(imp->m_AggInterpType , "agg_interp_type"); ADD_PARAMETER(imp->m_AggNumLevels , "agg_num_levels"); ADD_PARAMETER(imp->m_nodal , "nodal"); ADD_PARAMETER(imp->m_jacobi_pc , "do_jacobi"); @@ -396,6 +400,11 @@ void BoomerAMGSolver::SetNumSweeps(int nswp) imp->m_NumSweeps = nswp; } +void BoomerAMGSolver::SetAggInterpType(int aggit) +{ + imp->m_AggInterpType = aggit; +} + void BoomerAMGSolver::SetAggNumLevels(int anlv) { imp->m_AggNumLevels = anlv; @@ -515,6 +524,7 @@ void BoomerAMGSolver::SetInterpType(int inptyp) {} void BoomerAMGSolver::SetStrongThreshold(double thresh) {} void BoomerAMGSolver::SetPMaxElmts(int pmax) {} void BoomerAMGSolver::SetNumSweeps(int nswp) {} +void BoomerAMGSolver::SetAggInterpType(int aggit) {} void BoomerAMGSolver::SetAggNumLevels(int anlv) {} SparseMatrix* BoomerAMGSolver::CreateSparseMatrix(Matrix_Type ntype) { return nullptr; } bool BoomerAMGSolver::SetSparseMatrix(SparseMatrix* pA) { return false; } diff --git a/NumCore/BoomerAMGSolver.h b/NumCore/BoomerAMGSolver.h index 950d1be5a..dc1d4ba4b 100644 --- a/NumCore/BoomerAMGSolver.h +++ b/NumCore/BoomerAMGSolver.h @@ -81,6 +81,12 @@ class BoomerAMGSolver : public LinearSolver COARSE_OP_21 = 21, // CGC coarsening COARSE_OP_22 = 22, // CGC-E coarsening }; + + enum AggInterpType { + AGG_IT_0 = 0, + AGG_IT_5 = 5, + AGG_IT_7 = 7, + }; public: BoomerAMGSolver(FEModel* fem); @@ -98,6 +104,7 @@ class BoomerAMGSolver : public LinearSolver void SetStrongThreshold(double thresh); void SetPMaxElmts(int pmax); void SetNumSweeps(int nswp); + void SetAggInterpType(int aggit); void SetAggNumLevels(int anlv); void SetNodal(int nodal); void SetJacobiPC(bool b); diff --git a/NumCore/MKLDSSolver.cpp b/NumCore/MKLDSSolver.cpp new file mode 100644 index 000000000..1f209dc5c --- /dev/null +++ b/NumCore/MKLDSSolver.cpp @@ -0,0 +1,171 @@ +/*This file is part of the FEBio source code and is licensed under the MIT license +listed below. + +See Copyright-FEBio.txt for details. + +Copyright (c) 2021 University of Utah, The Trustees of Columbia University in +the City of New York, and others. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE.*/ +#include "stdafx.h" +#include +#include +#include "MKLDSSolver.h" +#include "CompactUnSymmMatrix.h" +#include "CompactSymmMatrix.h" + +#ifdef PARDISO +#undef PARDISO +#include "mkl_dss.h" +#include "mkl_types.h" +#define PARDISO +#endif + +BEGIN_FECORE_CLASS(MKLDSSolver, LinearSolver) +END_FECORE_CLASS(); + +#ifdef PARDISO +class MKLDSSolver::Imp +{ +public: + CompactMatrix* A = nullptr; + _MKL_DSS_HANDLE_t handle = nullptr; + MKL_INT opt = MKL_DSS_DEFAULTS; + MKL_INT sym = MKL_DSS_SYMMETRIC; +}; + +//----------------------------------------------------------------------------- +MKLDSSolver::MKLDSSolver(FEModel* fem) : LinearSolver(fem), m(new MKLDSSolver::Imp) +{ +} + +//----------------------------------------------------------------------------- +MKLDSSolver::~MKLDSSolver() +{ + Destroy(); + delete m; +} + +//----------------------------------------------------------------------------- +SparseMatrix* MKLDSSolver::CreateSparseMatrix(Matrix_Type ntype) +{ + // allocate the correct matrix format depending on matrix symmetry type + switch (ntype) + { + case REAL_SYMMETRIC : m->A = new CompactSymmMatrix(1); m->sym = MKL_DSS_SYMMETRIC; break; + case REAL_UNSYMMETRIC : m->A = new CRSSparseMatrix(1); m->sym = MKL_DSS_NON_SYMMETRIC; break; + case REAL_SYMM_STRUCTURE: m->A = new CRSSparseMatrix(1); m->sym = MKL_DSS_SYMMETRIC_STRUCTURE; break; + default: + assert(false); + m->A = nullptr; + } + + return m->A; +} + +//----------------------------------------------------------------------------- +bool MKLDSSolver::SetSparseMatrix(SparseMatrix* pA) +{ + return (m->A != nullptr); +} + +//----------------------------------------------------------------------------- +bool MKLDSSolver::PreProcess() +{ + SparseMatrix& A = *m->A; + int nRows = A.Rows(); + int nCols = A.Columns(); + int nnz = A.NonZeroes(); + + // create handle + if (m->handle == nullptr) + { + MKL_INT error = dss_create(m->handle, m->opt); + if (error != MKL_DSS_SUCCESS) return false; + } + + // define the structure + MKL_INT error = dss_define_structure(m->handle, m->sym, A.Pointers(), nRows, nCols, A.Indices(), nnz); + if (error != MKL_DSS_SUCCESS) return false; + + // reorder +// m->opt = MKL_DSS_AUTO_ORDER; +// m->opt = MKL_DSS_METIS_ORDER; + m->opt = MKL_DSS_METIS_OPENMP_ORDER; + error = dss_reorder(m->handle, m->opt, 0); + if (error != MKL_DSS_SUCCESS) return false; + + return LinearSolver::PreProcess(); +} + +//----------------------------------------------------------------------------- +bool MKLDSSolver::Factor() +{ + // make sure we have work to do + if ((m->A == nullptr) || (m->A->Rows() == 0)) return true; + + SparseMatrix& A = *m->A; +// MKL_INT type = MKL_DSS_POSITIVE_DEFINITE; + MKL_INT type = MKL_DSS_INDEFINITE; + MKL_INT error = dss_factor_real(m->handle, type, A.Values()); + if (error != MKL_DSS_SUCCESS) return false; + + return true; +} + +//----------------------------------------------------------------------------- +bool MKLDSSolver::BackSolve(double* x, double* b) +{ + // make sure we have work to do + if ((m->A == nullptr) || (m->A->Rows() == 0)) return true; + + int nRhs = 1; + m->opt = MKL_DSS_DEFAULTS; + MKL_INT error = dss_solve_real(m->handle, m->opt, b, nRhs, x); + if (error != MKL_DSS_SUCCESS) return false; + + // update stats + UpdateStats(1); + + return true; +} + +//----------------------------------------------------------------------------- +void MKLDSSolver::Destroy() +{ + if (m->handle) + { + m->opt = MKL_DSS_DEFAULTS; + MKL_INT error = dss_delete(m->handle, m->opt); + assert(error == MKL_DSS_SUCCESS); + m->handle = nullptr; + } +} +#else +//----------------------------------------------------------------------------- +class MKLDSSolver::Imp {}; +MKLDSSolver::MKLDSSolver(FEModel* fem) : LinearSolver(fem), m(nullptr) {} +MKLDSSolver::~MKLDSSolver() {} +SparseMatrix* MKLDSSolver::CreateSparseMatrix(Matrix_Type ntype) { return nullptr; } +bool MKLDSSolver::SetSparseMatrix(SparseMatrix* pA) { return false; } +bool MKLDSSolver::PreProcess() { return false; } +bool MKLDSSolver::Factor() { return false; } +bool MKLDSSolver::BackSolve(double* x, double* b) { return false; } +void MKLDSSolver::Destroy() {} +#endif diff --git a/NumCore/MKLDSSolver.h b/NumCore/MKLDSSolver.h new file mode 100644 index 000000000..387b2fe70 --- /dev/null +++ b/NumCore/MKLDSSolver.h @@ -0,0 +1,49 @@ +/*This file is part of the FEBio source code and is licensed under the MIT license +listed below. + +See Copyright-FEBio.txt for details. + +Copyright (c) 2021 University of Utah, The Trustees of Columbia University in +the City of New York, and others. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE.*/ +#pragma once +#include + +class MKLDSSolver : public LinearSolver +{ + class Imp; + +public: + MKLDSSolver(FEModel* fem); + ~MKLDSSolver(); + bool PreProcess() override; + bool Factor() override; + bool BackSolve(double* x, double* y) override; + void Destroy() override; + + SparseMatrix* CreateSparseMatrix(Matrix_Type ntype) override; + bool SetSparseMatrix(SparseMatrix* pA) override; + +protected: + Imp* m; + + DECLARE_FECORE_CLASS(); +}; + diff --git a/NumCore/NumCore.cpp b/NumCore/NumCore.cpp index c0feec992..ffba9f58d 100644 --- a/NumCore/NumCore.cpp +++ b/NumCore/NumCore.cpp @@ -47,6 +47,7 @@ SOFTWARE.*/ #include "FEASTEigenSolver.h" #include "TestSolver.h" #include "AccelerateSparseSolver.h" +#include "MKLDSSolver.h" //============================================================================= // Call this to initialize the NumCore module @@ -68,6 +69,7 @@ void NumCore::InitModule() REGISTER_FECORE_CLASS(StrategySolver , "strategy"); REGISTER_FECORE_CLASS(TestSolver , "test"); REGISTER_FECORE_CLASS(AccelerateSparseSolver, "accelerate"); + REGISTER_FECORE_CLASS(MKLDSSolver, "mkl_dss"); // register preconditioners REGISTER_FECORE_CLASS(ILU0_Preconditioner, "ilu0");