223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'typehints_orientation_rotation' of git.damask.mpie.de:damask/DAMASK into typehints_orientation_rotation

This commit is contained in:
Daniel Otto de Mentock 2022-01-28 13:40:47 +01:00
parent 53a0de2271 c3eb5a41b3
commit 83145b815f
40 changed files with 1263 additions and 816 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.github/workflows/Python.yml vendored
View File

@ -9,7 +9,7 @@ jobs:
strategy: strategy:
matrix: matrix:
python-version: ['3.7', '3.8', '3.9'] #, '3.10'] python-version: ['3.8', '3.9'] #, '3.10']
os: [ubuntu-latest, macos-latest] os: [ubuntu-latest, macos-latest]
steps: steps:

22
.gitlab-ci.yml
View File

@ -46,7 +46,7 @@ variables:
# ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1" PETSC_GNU: "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0" PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
PETSC_INTEL: "Libraries/PETSc/3.16.2/Intel-2022.0.1-IntelMPI-2021.5.0" PETSC_INTEL: "Libraries/PETSc/3.16.3/Intel-2022.0.1-IntelMPI-2021.5.0"
# ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MSC: "FEM/MSC/2021.3.1" MSC: "FEM/MSC/2021.3.1"
IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020" IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
@ -79,63 +79,63 @@ mypy:
################################################################################################### ###################################################################################################
test_grid_GNU: grid_GNU:
stage: compile stage: compile
script: script:
- module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU} - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU}
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU
test_mesh_GNU: mesh_GNU:
stage: compile stage: compile
script: script:
- module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU} - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU}
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU
test_grid_GNU-64bit: grid_GNU-64bit:
stage: compile stage: compile
script: script:
- module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit - module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU-64bit - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU-64bit
test_mesh_GNU-64bit: mesh_GNU-64bit:
stage: compile stage: compile
script: script:
- module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit - module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU-64bit - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU-64bit
test_grid_IntelLLVM: grid_IntelLLVM:
stage: compile stage: compile
script: script:
- module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM} - module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM}
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_IntelLLVM - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_IntelLLVM
test_mesh_IntelLLVM: mesh_IntelLLVM:
stage: compile stage: compile
script: script:
- module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM} - module load ${COMPILER_INTELLLVM} ${MPI_INTELLLVM} ${PETSC_INTELLLVM}
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_IntelLLVM - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_IntelLLVM
test_grid_Intel: grid_Intel:
stage: compile stage: compile
script: script:
- module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL} - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_Intel - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_Intel
test_mesh_Intel: mesh_Intel:
stage: compile stage: compile
script: script:
- module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL} - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_Intel - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_Intel
test_Marc: Marc_Intel:
stage: compile stage: compile
script: script:
- module load $IntelMarc $HDF5Marc $MSC - module load $IntelMarc $HDF5Marc $MSC
@ -158,7 +158,7 @@ setup_mesh:
- cmake -DDAMASK_SOLVER=MESH -DCMAKE_INSTALL_PREFIX=${TESTROOT} ${CI_PROJECT_DIR} - cmake -DDAMASK_SOLVER=MESH -DCMAKE_INSTALL_PREFIX=${TESTROOT} ${CI_PROJECT_DIR}
- make -j2 all install - make -j2 all install
compile_Marc: setup_Marc:
stage: compile stage: compile
script: script:
- module load $IntelMarc $HDF5Marc $MSC - module load $IntelMarc $HDF5Marc $MSC

110
COPYING
View File

@ -1,23 +1,21 @@
GNU GENERAL PUBLIC LICENSE GNU AFFERO GENERAL PUBLIC LICENSE
Version 3, 29 June 2007 Version 3, 19 November 2007
Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/> Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
Everyone is permitted to copy and distribute verbatim copies Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed. of this license document, but changing it is not allowed.
Preamble Preamble
The GNU General Public License is a free, copyleft license for The GNU Affero General Public License is a free, copyleft license for
software and other kinds of works. software and other kinds of works, specifically designed to ensure
cooperation with the community in the case of network server software.
The licenses for most software and other practical works are designed The licenses for most software and other practical works are designed
to take away your freedom to share and change the works. By contrast, to take away your freedom to share and change the works. By contrast,
the GNU General Public License is intended to guarantee your freedom to our General Public Licenses are intended to guarantee your freedom to
share and change all versions of a program--to make sure it remains free share and change all versions of a program--to make sure it remains free
software for all its users. We, the Free Software Foundation, use the software for all its users.
GNU General Public License for most of our software; it applies also to
any other work released this way by its authors. You can apply it to
your programs, too.
When we speak of free software, we are referring to freedom, not When we speak of free software, we are referring to freedom, not
price. Our General Public Licenses are designed to make sure that you price. Our General Public Licenses are designed to make sure that you
@ -26,44 +24,34 @@ them if you wish), that you receive source code or can get it if you
want it, that you can change the software or use pieces of it in new want it, that you can change the software or use pieces of it in new
free programs, and that you know you can do these things. free programs, and that you know you can do these things.
To protect your rights, we need to prevent others from denying you Developers that use our General Public Licenses protect your rights
these rights or asking you to surrender the rights. Therefore, you have with two steps: (1) assert copyright on the software, and (2) offer
certain responsibilities if you distribute copies of the software, or if you this License which gives you legal permission to copy, distribute
you modify it: responsibilities to respect the freedom of others. and/or modify the software.
For example, if you distribute copies of such a program, whether A secondary benefit of defending all users' freedom is that
gratis or for a fee, you must pass on to the recipients the same improvements made in alternate versions of the program, if they
freedoms that you received. You must make sure that they, too, receive receive widespread use, become available for other developers to
or can get the source code. And you must show them these terms so they incorporate. Many developers of free software are heartened and
know their rights. encouraged by the resulting cooperation. However, in the case of
software used on network servers, this result may fail to come about.
The GNU General Public License permits making a modified version and
letting the public access it on a server without ever releasing its
source code to the public.
Developers that use the GNU GPL protect your rights with two steps: The GNU Affero General Public License is designed specifically to
(1) assert copyright on the software, and (2) offer you this License ensure that, in such cases, the modified source code becomes available
giving you legal permission to copy, distribute and/or modify it. to the community. It requires the operator of a network server to
provide the source code of the modified version running there to the
users of that server. Therefore, public use of a modified version, on
a publicly accessible server, gives the public access to the source
code of the modified version.
For the developers' and authors' protection, the GPL clearly explains An older license, called the Affero General Public License and
that there is no warranty for this free software. For both users' and published by Affero, was designed to accomplish similar goals. This is
authors' sake, the GPL requires that modified versions be marked as a different license, not a version of the Affero GPL, but Affero has
changed, so that their problems will not be attributed erroneously to released a new version of the Affero GPL which permits relicensing under
authors of previous versions. this license.
Some devices are designed to deny users access to install or run
modified versions of the software inside them, although the manufacturer
can do so. This is fundamentally incompatible with the aim of
protecting users' freedom to change the software. The systematic
pattern of such abuse occurs in the area of products for individuals to
use, which is precisely where it is most unacceptable. Therefore, we
have designed this version of the GPL to prohibit the practice for those
products. If such problems arise substantially in other domains, we
stand ready to extend this provision to those domains in future versions
of the GPL, as needed to protect the freedom of users.
Finally, every program is threatened constantly by software patents.
States should not allow patents to restrict development and use of
software on general-purpose computers, but in those that do, we wish to
avoid the special danger that patents applied to a free program could
make it effectively proprietary. To prevent this, the GPL assures that
patents cannot be used to render the program non-free.
The precise terms and conditions for copying, distribution and The precise terms and conditions for copying, distribution and
modification follow. modification follow.
@ -72,7 +60,7 @@ modification follow.
0. Definitions. 0. Definitions.
"This License" refers to version 3 of the GNU General Public License. "This License" refers to version 3 of the GNU Affero General Public License.
"Copyright" also means copyright-like laws that apply to other kinds of "Copyright" also means copyright-like laws that apply to other kinds of
works, such as semiconductor masks. works, such as semiconductor masks.
@ -549,35 +537,45 @@ to collect a royalty for further conveying from those to whom you convey
the Program, the only way you could satisfy both those terms and this the Program, the only way you could satisfy both those terms and this
License would be to refrain entirely from conveying the Program. License would be to refrain entirely from conveying the Program.
13. Use with the GNU Affero General Public License. 13. Remote Network Interaction; Use with the GNU General Public License.
Notwithstanding any other provision of this License, if you modify the
Program, your modified version must prominently offer all users
interacting with it remotely through a computer network (if your version
supports such interaction) an opportunity to receive the Corresponding
Source of your version by providing access to the Corresponding Source
from a network server at no charge, through some standard or customary
means of facilitating copying of software. This Corresponding Source
shall include the Corresponding Source for any work covered by version 3
of the GNU General Public License that is incorporated pursuant to the
following paragraph.
Notwithstanding any other provision of this License, you have Notwithstanding any other provision of this License, you have
permission to link or combine any covered work with a work licensed permission to link or combine any covered work with a work licensed
under version 3 of the GNU Affero General Public License into a single under version 3 of the GNU General Public License into a single
combined work, and to convey the resulting work. The terms of this combined work, and to convey the resulting work. The terms of this
License will continue to apply to the part which is the covered work, License will continue to apply to the part which is the covered work,
but the special requirements of the GNU Affero General Public License, but the work with which it is combined will remain governed by version
section 13, concerning interaction through a network will apply to the 3 of the GNU General Public License.
combination as such.
14. Revised Versions of this License. 14. Revised Versions of this License.
The Free Software Foundation may publish revised and/or new versions of The Free Software Foundation may publish revised and/or new versions of
the GNU General Public License from time to time. Such new versions will the GNU Affero General Public License from time to time. Such new versions
be similar in spirit to the present version, but may differ in detail to will be similar in spirit to the present version, but may differ in detail to
address new problems or concerns. address new problems or concerns.
Each version is given a distinguishing version number. If the Each version is given a distinguishing version number. If the
Program specifies that a certain numbered version of the GNU General Program specifies that a certain numbered version of the GNU Affero General
Public License "or any later version" applies to it, you have the Public License "or any later version" applies to it, you have the
option of following the terms and conditions either of that numbered option of following the terms and conditions either of that numbered
version or of any later version published by the Free Software version or of any later version published by the Free Software
Foundation. If the Program does not specify a version number of the Foundation. If the Program does not specify a version number of the
GNU General Public License, you may choose any version ever published GNU Affero General Public License, you may choose any version ever published
by the Free Software Foundation. by the Free Software Foundation.
If the Program specifies that a proxy can decide which future If the Program specifies that a proxy can decide which future
versions of the GNU General Public License can be used, that proxy's versions of the GNU Affero General Public License can be used, that proxy's
public statement of acceptance of a version permanently authorizes you public statement of acceptance of a version permanently authorizes you
to choose that version for the Program. to choose that version for the Program.

10
LICENSE
View File

@ -1,14 +1,14 @@
Copyright 2011-2021 Max-Planck-Institut für Eisenforschung GmbH Copyright 2011-2022 Max-Planck-Institut für Eisenforschung GmbH
DAMASK is free software: you can redistribute it and/or modify DAMASK is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU Affero General Public License as published by
the Free Software Foundation, either version 3 of the License, or the Free Software Foundation, either version 3 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU Affero General Public License for more details.
You should have received a copy of the GNU General Public License You should have received a copy of the GNU Affero General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>. along with this program. If not, see <https://www.gnu.org/licenses/>.

2
PRIVATE

@ -1 +1 @@
Subproject commit b898a8b5552bd9d1c555edc3d8134564dd32fe53 Subproject commit 5598ec4892b0921fccf63e8570f9fcd8e0365716

12
examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml
View File

@ -6,7 +6,7 @@ references:
output: [xi_sl, xi_tw] output: [xi_sl, xi_tw]
N_sl: [3, 3, 6, 0, 6] # basal, prism, -, 1. pyr<a>, -, 2. pyr<c+a> N_sl: [3, 3, 6, 0, 6] # basal, prism, 1. pyr<a>, -, 2. pyr<c+a>
N_tw: [6, 0, 6] # tension, -, compression N_tw: [6, 0, 6] # tension, -, compression
xi_0_sl: [10.e+6, 55.e+6, 60.e+6, 0., 60.e+6] xi_0_sl: [10.e+6, 55.e+6, 60.e+6, 0., 60.e+6]
@ -23,9 +23,13 @@ f_sat_sl-tw: 10.0
h_0_sl-sl: 500.0e+6 h_0_sl-sl: 500.0e+6
h_0_tw-tw: 50.0e+6 h_0_tw-tw: 50.0e+6
h_0_tw-sl: 150.0e+6 h_0_tw-sl: 150.0e+6
h_sl-sl: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, +1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0,
1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, +1.0, 1.0] # unused entries are indicated by -1.0 -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0,
-1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # unused entries are indicated by -1.0
h_tw-tw: [+1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, h_tw-tw: [+1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0,
-1.0, 1.0] # unused entries are indicated by -1.0 -1.0, 1.0] # unused entries are indicated by -1.0
h_tw-sl: [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, h_tw-sl: [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0,

9
examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml
View File

@ -8,7 +8,7 @@ references:
https://doi.org/10.1016/j.actamat.2017.05.015 https://doi.org/10.1016/j.actamat.2017.05.015
output: [gamma_sl] output: [gamma_sl]
N_sl: [3, 3, 0, 12] # basal, 1. prism, -, 1. pyr<c+a> N_sl: [3, 3, 0, 12] # basal, prism, -, 1. pyr<c+a>
n_sl: 20 n_sl: 20
a_sl: 2.0 a_sl: 2.0
dot_gamma_0_sl: 0.001 dot_gamma_0_sl: 0.001
@ -20,5 +20,8 @@ xi_inf_sl: [568.e+6, 150.e+7, 0.0, 3420.e+6]
# L. Wang et al. : # L. Wang et al. :
# xi_0_sl: [127.e+6, 96.e+6, 0.0, 240.e+6] # xi_0_sl: [127.e+6, 96.e+6, 0.0, 240.e+6]
h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, h_sl-sl: [+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0,
-1.0, -1.0, +1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0] # unused entries are indicated by -1.0 -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # unused entries are indicated by -1.0

358
img/COPYING
View File

@ -1,94 +1,312 @@
Creative Commons Attribution-NoDerivatives 4.0 International Public License Creative Commons Attribution-NoDerivatives 4.0 International Public
License
By exercising the Licensed Rights (defined below), You accept and agree to be bound by the terms and conditions of this Creative Commons Attribution-NoDerivatives 4.0 International Public License ("Public License"). To the extent this Public License may be interpreted as a contract, You are granted the Licensed Rights in consideration of Your acceptance of these terms and conditions, and the Licensor grants You such rights in consideration of benefits the Licensor receives from making the Licensed Material available under these terms and conditions. By exercising the Licensed Rights (defined below), You accept and agree
to be bound by the terms and conditions of this Creative Commons
Attribution-NoDerivatives 4.0 International Public License ("Public
License"). To the extent this Public License may be interpreted as a
contract, You are granted the Licensed Rights in consideration of Your
acceptance of these terms and conditions, and the Licensor grants You
such rights in consideration of benefits the Licensor receives from
making the Licensed Material available under these terms and
conditions.
Section 1 – Definitions.
Adapted Material means material subject to Copyright and Similar Rights that is derived from or based upon the Licensed Material and in which the Licensed Material is translated, altered, arranged, transformed, or otherwise modified in a manner requiring permission under the Copyright and Similar Rights held by the Licensor. For purposes of this Public License, where the Licensed Material is a musical work, performance, or sound recording, Adapted Material is always produced where the Licensed Material is synched in timed relation with a moving image. Section 1 -- Definitions.
Copyright and Similar Rights means copyright and/or similar rights closely related to copyright including, without limitation, performance, broadcast, sound recording, and Sui Generis Database Rights, without regard to how the rights are labeled or categorized. For purposes of this Public License, the rights specified in Section 2(b)(1)-(2) are not Copyright and Similar Rights.
Effective Technological Measures means those measures that, in the absence of proper authority, may not be circumvented under laws fulfilling obligations under Article 11 of the WIPO Copyright Treaty adopted on December 20, 1996, and/or similar international agreements.
Exceptions and Limitations means fair use, fair dealing, and/or any other exception or limitation to Copyright and Similar Rights that applies to Your use of the Licensed Material.
Licensed Material means the artistic or literary work, database, or other material to which the Licensor applied this Public License.
Licensed Rights means the rights granted to You subject to the terms and conditions of this Public License, which are limited to all Copyright and Similar Rights that apply to Your use of the Licensed Material and that the Licensor has authority to license.
Licensor means the individual(s) or entity(ies) granting rights under this Public License.
Share means to provide material to the public by any means or process that requires permission under the Licensed Rights, such as reproduction, public display, public performance, distribution, dissemination, communication, or importation, and to make material available to the public including in ways that members of the public may access the material from a place and at a time individually chosen by them.
Sui Generis Database Rights means rights other than copyright resulting from Directive 96/9/EC of the European Parliament and of the Council of 11 March 1996 on the legal protection of databases, as amended and/or succeeded, as well as other essentially equivalent rights anywhere in the world.
You means the individual or entity exercising the Licensed Rights under this Public License. Your has a corresponding meaning.
Section 2 – Scope. a. Adapted Material means material subject to Copyright and Similar
Rights that is derived from or based upon the Licensed Material
and in which the Licensed Material is translated, altered,
arranged, transformed, or otherwise modified in a manner requiring
permission under the Copyright and Similar Rights held by the
Licensor. For purposes of this Public License, where the Licensed
Material is a musical work, performance, or sound recording,
Adapted Material is always produced where the Licensed Material is
synched in timed relation with a moving image.
License grant. b. Copyright and Similar Rights means copyright and/or similar rights
Subject to the terms and conditions of this Public License, the Licensor hereby grants You a worldwide, royalty-free, non-sublicensable, non-exclusive, irrevocable license to exercise the Licensed Rights in the Licensed Material to: closely related to copyright including, without limitation,
reproduce and Share the Licensed Material, in whole or in part; and performance, broadcast, sound recording, and Sui Generis Database
produce and reproduce, but not Share, Adapted Material. Rights, without regard to how the rights are labeled or
Exceptions and Limitations. For the avoidance of doubt, where Exceptions and Limitations apply to Your use, this Public License does not apply, and You do not need to comply with its terms and conditions. categorized. For purposes of this Public License, the rights
Term. The term of this Public License is specified in Section 6(a). specified in Section 2(b)(1)-(2) are not Copyright and Similar
Media and formats; technical modifications allowed. The Licensor authorizes You to exercise the Licensed Rights in all media and formats whether now known or hereafter created, and to make technical modifications necessary to do so. The Licensor waives and/or agrees not to assert any right or authority to forbid You from making technical modifications necessary to exercise the Licensed Rights, including technical modifications necessary to circumvent Effective Technological Measures. For purposes of this Public License, simply making modifications authorized by this Section 2(a)(4) never produces Adapted Material. Rights.
Downstream recipients.
Offer from the Licensor – Licensed Material. Every recipient of the Licensed Material automatically receives an offer from the Licensor to exercise the Licensed Rights under the terms and conditions of this Public License.
No downstream restrictions. You may not offer or impose any additional or different terms or conditions on, or apply any Effective Technological Measures to, the Licensed Material if doing so restricts exercise of the Licensed Rights by any recipient of the Licensed Material.
No endorsement. Nothing in this Public License constitutes or may be construed as permission to assert or imply that You are, or that Your use of the Licensed Material is, connected with, or sponsored, endorsed, or granted official status by, the Licensor or others designated to receive attribution as provided in Section 3(a)(1)(A)(i).
Other rights. c. Effective Technological Measures means those measures that, in the
Moral rights, such as the right of integrity, are not licensed under this Public License, nor are publicity, privacy, and/or other similar personality rights; however, to the extent possible, the Licensor waives and/or agrees not to assert any such rights held by the Licensor to the limited extent necessary to allow You to exercise the Licensed Rights, but not otherwise. absence of proper authority, may not be circumvented under laws
Patent and trademark rights are not licensed under this Public License. fulfilling obligations under Article 11 of the WIPO Copyright
To the extent possible, the Licensor waives any right to collect royalties from You for the exercise of the Licensed Rights, whether directly or through a collecting society under any voluntary or waivable statutory or compulsory licensing scheme. In all other cases the Licensor expressly reserves any right to collect such royalties. Treaty adopted on December 20, 1996, and/or similar international
agreements.
Section 3 – License Conditions. d. Exceptions and Limitations means fair use, fair dealing, and/or
any other exception or limitation to Copyright and Similar Rights
that applies to Your use of the Licensed Material.
Your exercise of the Licensed Rights is expressly made subject to the following conditions. e. Licensed Material means the artistic or literary work, database,
or other material to which the Licensor applied this Public
License.
Attribution. f. Licensed Rights means the rights granted to You subject to the
terms and conditions of this Public License, which are limited to
all Copyright and Similar Rights that apply to Your use of the
Licensed Material and that the Licensor has authority to license.
If You Share the Licensed Material, You must: g. Licensor means the individual(s) or entity(ies) granting rights
retain the following if it is supplied by the Licensor with the Licensed Material: under this Public License.
identification of the creator(s) of the Licensed Material and any others designated to receive attribution, in any reasonable manner requested by the Licensor (including by pseudonym if designated);
a copyright notice;
a notice that refers to this Public License;
a notice that refers to the disclaimer of warranties;
a URI or hyperlink to the Licensed Material to the extent reasonably practicable;
indicate if You modified the Licensed Material and retain an indication of any previous modifications; and
indicate the Licensed Material is licensed under this Public License, and include the text of, or the URI or hyperlink to, this Public License.
For the avoidance of doubt, You do not have permission under this Public License to Share Adapted Material.
You may satisfy the conditions in Section 3(a)(1) in any reasonable manner based on the medium, means, and context in which You Share the Licensed Material. For example, it may be reasonable to satisfy the conditions by providing a URI or hyperlink to a resource that includes the required information.
If requested by the Licensor, You must remove any of the information required by Section 3(a)(1)(A) to the extent reasonably practicable.
Section 4 – Sui Generis Database Rights. h. Share means to provide material to the public by any means or
process that requires permission under the Licensed Rights, such
as reproduction, public display, public performance, distribution,
dissemination, communication, or importation, and to make material
available to the public including in ways that members of the
public may access the material from a place and at a time
individually chosen by them.
Where the Licensed Rights include Sui Generis Database Rights that apply to Your use of the Licensed Material: i. Sui Generis Database Rights means rights other than copyright
resulting from Directive 96/9/EC of the European Parliament and of
the Council of 11 March 1996 on the legal protection of databases,
as amended and/or succeeded, as well as other essentially
equivalent rights anywhere in the world.
for the avoidance of doubt, Section 2(a)(1) grants You the right to extract, reuse, reproduce, and Share all or a substantial portion of the contents of the database, provided You do not Share Adapted Material; j. You means the individual or entity exercising the Licensed Rights
if You include all or a substantial portion of the database contents in a database in which You have Sui Generis Database Rights, then the database in which You have Sui Generis Database Rights (but not its individual contents) is Adapted Material; and under this Public License. Your has a corresponding meaning.
You must comply with the conditions in Section 3(a) if You Share all or a substantial portion of the contents of the database.
For the avoidance of doubt, this Section 4 supplements and does not replace Your obligations under this Public License where the Licensed Rights include other Copyright and Similar Rights.
Section 5 – Disclaimer of Warranties and Limitation of Liability. Section 2 -- Scope.
Unless otherwise separately undertaken by the Licensor, to the extent possible, the Licensor offers the Licensed Material as-is and as-available, and makes no representations or warranties of any kind concerning the Licensed Material, whether express, implied, statutory, or other. This includes, without limitation, warranties of title, merchantability, fitness for a particular purpose, non-infringement, absence of latent or other defects, accuracy, or the presence or absence of errors, whether or not known or discoverable. Where disclaimers of warranties are not allowed in full or in part, this disclaimer may not apply to You. a. License grant.
To the extent possible, in no event will the Licensor be liable to You on any legal theory (including, without limitation, negligence) or otherwise for any direct, special, indirect, incidental, consequential, punitive, exemplary, or other losses, costs, expenses, or damages arising out of this Public License or use of the Licensed Material, even if the Licensor has been advised of the possibility of such losses, costs, expenses, or damages. Where a limitation of liability is not allowed in full or in part, this limitation may not apply to You.
The disclaimer of warranties and limitation of liability provided above shall be interpreted in a manner that, to the extent possible, most closely approximates an absolute disclaimer and waiver of all liability. 1. Subject to the terms and conditions of this Public License,
the Licensor hereby grants You a worldwide, royalty-free,
non-sublicensable, non-exclusive, irrevocable license to
exercise the Licensed Rights in the Licensed Material to:
Section 6 – Term and Termination. a. reproduce and Share the Licensed Material, in whole or
in part; and
This Public License applies for the term of the Copyright and Similar Rights licensed here. However, if You fail to comply with this Public License, then Your rights under this Public License terminate automatically. b. produce and reproduce, but not Share, Adapted Material.
Where Your right to use the Licensed Material has terminated under Section 6(a), it reinstates: 2. Exceptions and Limitations. For the avoidance of doubt, where
automatically as of the date the violation is cured, provided it is cured within 30 days of Your discovery of the violation; or Exceptions and Limitations apply to Your use, this Public
upon express reinstatement by the Licensor. License does not apply, and You do not need to comply with
For the avoidance of doubt, this Section 6(b) does not affect any right the Licensor may have to seek remedies for Your violations of this Public License. its terms and conditions.
For the avoidance of doubt, the Licensor may also offer the Licensed Material under separate terms or conditions or stop distributing the Licensed Material at any time; however, doing so will not terminate this Public License.
Sections 1, 5, 6, 7, and 8 survive termination of this Public License.
Section 7 – Other Terms and Conditions. 3. Term. The term of this Public License is specified in Section
6(a).
The Licensor shall not be bound by any additional or different terms or conditions communicated by You unless expressly agreed. 4. Media and formats; technical modifications allowed. The
Any arrangements, understandings, or agreements regarding the Licensed Material not stated herein are separate from and independent of the terms and conditions of this Public License. Licensor authorizes You to exercise the Licensed Rights in
all media and formats whether now known or hereafter created,
and to make technical modifications necessary to do so. The
Licensor waives and/or agrees not to assert any right or
authority to forbid You from making technical modifications
necessary to exercise the Licensed Rights, including
technical modifications necessary to circumvent Effective
Technological Measures. For purposes of this Public License,
simply making modifications authorized by this Section 2(a)
(4) never produces Adapted Material.
Section 8 – Interpretation. 5. Downstream recipients.
For the avoidance of doubt, this Public License does not, and shall not be interpreted to, reduce, limit, restrict, or impose conditions on any use of the Licensed Material that could lawfully be made without permission under this Public License. a. Offer from the Licensor -- Licensed Material. Every
To the extent possible, if any provision of this Public License is deemed unenforceable, it shall be automatically reformed to the minimum extent necessary to make it enforceable. If the provision cannot be reformed, it shall be severed from this Public License without affecting the enforceability of the remaining terms and conditions. recipient of the Licensed Material automatically
No term or condition of this Public License will be waived and no failure to comply consented to unless expressly agreed to by the Licensor. receives an offer from the Licensor to exercise the
Nothing in this Public License constitutes or may be interpreted as a limitation upon, or waiver of, any privileges and immunities that apply to the Licensor or You, including from the legal processes of any jurisdiction or authority. Licensed Rights under the terms and conditions of this
Public License.
b. No downstream restrictions. You may not offer or impose
any additional or different terms or conditions on, or
apply any Effective Technological Measures to, the
Licensed Material if doing so restricts exercise of the
Licensed Rights by any recipient of the Licensed
Material.
6. No endorsement. Nothing in this Public License constitutes or
may be construed as permission to assert or imply that You
are, or that Your use of the Licensed Material is, connected
with, or sponsored, endorsed, or granted official status by,
the Licensor or others designated to receive attribution as
provided in Section 3(a)(1)(A)(i).
b. Other rights.
1. Moral rights, such as the right of integrity, are not
licensed under this Public License, nor are publicity,
privacy, and/or other similar personality rights; however, to
the extent possible, the Licensor waives and/or agrees not to
assert any such rights held by the Licensor to the limited
extent necessary to allow You to exercise the Licensed
Rights, but not otherwise.
2. Patent and trademark rights are not licensed under this
Public License.
3. To the extent possible, the Licensor waives any right to
collect royalties from You for the exercise of the Licensed
Rights, whether directly or through a collecting society
under any voluntary or waivable statutory or compulsory
licensing scheme. In all other cases the Licensor expressly
reserves any right to collect such royalties.
Section 3 -- License Conditions.
Your exercise of the Licensed Rights is expressly made subject to the
following conditions.
a. Attribution.
1. If You Share the Licensed Material, You must:
a. retain the following if it is supplied by the Licensor
with the Licensed Material:
i. identification of the creator(s) of the Licensed
Material and any others designated to receive
attribution, in any reasonable manner requested by
the Licensor (including by pseudonym if
designated);
ii. a copyright notice;
iii. a notice that refers to this Public License;
iv. a notice that refers to the disclaimer of
warranties;
v. a URI or hyperlink to the Licensed Material to the
extent reasonably practicable;
b. indicate if You modified the Licensed Material and
retain an indication of any previous modifications; and
c. indicate the Licensed Material is licensed under this
Public License, and include the text of, or the URI or
hyperlink to, this Public License.
For the avoidance of doubt, You do not have permission under
this Public License to Share Adapted Material.
2. You may satisfy the conditions in Section 3(a)(1) in any
reasonable manner based on the medium, means, and context in
which You Share the Licensed Material. For example, it may be
reasonable to satisfy the conditions by providing a URI or
hyperlink to a resource that includes the required
information.
3. If requested by the Licensor, You must remove any of the
information required by Section 3(a)(1)(A) to the extent
reasonably practicable.
Section 4 -- Sui Generis Database Rights.
Where the Licensed Rights include Sui Generis Database Rights that
apply to Your use of the Licensed Material:
a. for the avoidance of doubt, Section 2(a)(1) grants You the right
to extract, reuse, reproduce, and Share all or a substantial
portion of the contents of the database, provided You do not Share
Adapted Material;
b. if You include all or a substantial portion of the database
contents in a database in which You have Sui Generis Database
Rights, then the database in which You have Sui Generis Database
Rights (but not its individual contents) is Adapted Material; and
c. You must comply with the conditions in Section 3(a) if You Share
all or a substantial portion of the contents of the database.
For the avoidance of doubt, this Section 4 supplements and does not
replace Your obligations under this Public License where the Licensed
Rights include other Copyright and Similar Rights.
Section 5 -- Disclaimer of Warranties and Limitation of Liability.
a. UNLESS OTHERWISE SEPARATELY UNDERTAKEN BY THE LICENSOR, TO THE
EXTENT POSSIBLE, THE LICENSOR OFFERS THE LICENSED MATERIAL AS-IS
AND AS-AVAILABLE, AND MAKES NO REPRESENTATIONS OR WARRANTIES OF
ANY KIND CONCERNING THE LICENSED MATERIAL, WHETHER EXPRESS,
IMPLIED, STATUTORY, OR OTHER. THIS INCLUDES, WITHOUT LIMITATION,
WARRANTIES OF TITLE, MERCHANTABILITY, FITNESS FOR A PARTICULAR
PURPOSE, NON-INFRINGEMENT, ABSENCE OF LATENT OR OTHER DEFECTS,
ACCURACY, OR THE PRESENCE OR ABSENCE OF ERRORS, WHETHER OR NOT
KNOWN OR DISCOVERABLE. WHERE DISCLAIMERS OF WARRANTIES ARE NOT
ALLOWED IN FULL OR IN PART, THIS DISCLAIMER MAY NOT APPLY TO YOU.
b. TO THE EXTENT POSSIBLE, IN NO EVENT WILL THE LICENSOR BE LIABLE
TO YOU ON ANY LEGAL THEORY (INCLUDING, WITHOUT LIMITATION,
NEGLIGENCE) OR OTHERWISE FOR ANY DIRECT, SPECIAL, INDIRECT,
INCIDENTAL, CONSEQUENTIAL, PUNITIVE, EXEMPLARY, OR OTHER LOSSES,
COSTS, EXPENSES, OR DAMAGES ARISING OUT OF THIS PUBLIC LICENSE OR
USE OF THE LICENSED MATERIAL, EVEN IF THE LICENSOR HAS BEEN
ADVISED OF THE POSSIBILITY OF SUCH LOSSES, COSTS, EXPENSES, OR
DAMAGES. WHERE A LIMITATION OF LIABILITY IS NOT ALLOWED IN FULL OR
IN PART, THIS LIMITATION MAY NOT APPLY TO YOU.
c. The disclaimer of warranties and limitation of liability provided
above shall be interpreted in a manner that, to the extent
possible, most closely approximates an absolute disclaimer and
waiver of all liability.
Section 6 -- Term and Termination.
a. This Public License applies for the term of the Copyright and
Similar Rights licensed here. However, if You fail to comply with
this Public License, then Your rights under this Public License
terminate automatically.
b. Where Your right to use the Licensed Material has terminated under
Section 6(a), it reinstates:
1. automatically as of the date the violation is cured, provided
it is cured within 30 days of Your discovery of the
violation; or
2. upon express reinstatement by the Licensor.
For the avoidance of doubt, this Section 6(b) does not affect any
right the Licensor may have to seek remedies for Your violations
of this Public License.
c. For the avoidance of doubt, the Licensor may also offer the
Licensed Material under separate terms or conditions or stop
distributing the Licensed Material at any time; however, doing so
will not terminate this Public License.
d. Sections 1, 5, 6, 7, and 8 survive termination of this Public
License.
Section 7 -- Other Terms and Conditions.
a. The Licensor shall not be bound by any additional or different
terms or conditions communicated by You unless expressly agreed.
b. Any arrangements, understandings, or agreements regarding the
Licensed Material not stated herein are separate from and
independent of the terms and conditions of this Public License.
Section 8 -- Interpretation.
a. For the avoidance of doubt, this Public License does not, and
shall not be interpreted to, reduce, limit, restrict, or impose
conditions on any use of the Licensed Material that could lawfully
be made without permission under this Public License.
b. To the extent possible, if any provision of this Public License is
deemed unenforceable, it shall be automatically reformed to the
minimum extent necessary to make it enforceable. If the provision
cannot be reformed, it shall be severed from this Public License
without affecting the enforceability of the remaining terms and
conditions.
c. No term or condition of this Public License will be waived and no
failure to comply consented to unless expressly agreed to by the
Licensor.
d. Nothing in this Public License constitutes or may be interpreted
as a limitation upon, or waiver of, any privileges and immunities
that apply to the Licensor or You, including from the legal
processes of any jurisdiction or authority.

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha5-389-ga000e477c v3.0.0-alpha5-518-g4fa97b9a3

2
python/damask/__init__.py
View File

@ -21,8 +21,8 @@ from ._rotation import Rotation # noqa
from ._crystal import Crystal # noqa from ._crystal import Crystal # noqa
from ._orientation import Orientation # noqa from ._orientation import Orientation # noqa
from ._table import Table # noqa from ._table import Table # noqa
from ._vtk import VTK # noqa
from ._colormap import Colormap # noqa from ._colormap import Colormap # noqa
from ._vtk import VTK # noqa
from ._config import Config # noqa from ._config import Config # noqa
from ._configmaterial import ConfigMaterial # noqa from ._configmaterial import ConfigMaterial # noqa
from ._grid import Grid # noqa from ._grid import Grid # noqa

10
python/damask/_colormap.py
View File

@ -63,6 +63,14 @@ class Colormap(mpl.colors.ListedColormap):
"""Concatenate (in-place).""" """Concatenate (in-place)."""
return self.__add__(other) return self.__add__(other)
def __mul__(self, factor: int) -> 'Colormap':
"""Repeat."""
return Colormap(np.vstack([self.colors]*factor),f'{self.name}*{factor}')
def __imul__(self, factor: int) -> 'Colormap':
"""Repeat (in-place)."""
return self.__mul__(factor)
def __invert__(self) -> 'Colormap': def __invert__(self) -> 'Colormap':
"""Reverse.""" """Reverse."""
return self.reversed() return self.reversed()
@ -250,7 +258,7 @@ class Colormap(mpl.colors.ListedColormap):
np.logical_or (np.isnan(field), field == gap)) # mask NaN (and gap if present) np.logical_or (np.isnan(field), field == gap)) # mask NaN (and gap if present)
l,r = (field[mask].min(),field[mask].max()) if bounds is None else \ l,r = (field[mask].min(),field[mask].max()) if bounds is None else \
np.array(bounds,float)[:2] (bounds[0],bounds[1])
delta,avg = r-l,0.5*abs(r+l) delta,avg = r-l,0.5*abs(r+l)

71
python/damask/_grid.py
View File

@ -4,7 +4,7 @@ import warnings
import multiprocessing as mp import multiprocessing as mp
from functools import partial from functools import partial
import typing import typing
from typing import Union, Optional, TextIO, List, Sequence from typing import Union, Optional, TextIO, List, Sequence, Literal
from pathlib import Path from pathlib import Path
import numpy as np import numpy as np
@ -70,7 +70,7 @@ class Grid:
]) ])
def __copy__(self) -> "Grid": def __copy__(self) -> 'Grid':
"""Create deep copy.""" """Create deep copy."""
return copy.deepcopy(self) return copy.deepcopy(self)
@ -161,7 +161,7 @@ class Grid:
@staticmethod @staticmethod
def load(fname: Union[str, Path]) -> "Grid": def load(fname: Union[str, Path]) -> 'Grid':
""" """
Load from VTK image data file. Load from VTK image data file.
@ -190,7 +190,7 @@ class Grid:
@typing. no_type_check @typing. no_type_check
@staticmethod @staticmethod
def load_ASCII(fname)-> "Grid": def load_ASCII(fname)-> 'Grid':
""" """
Load from geom file. Load from geom file.
@ -264,7 +264,7 @@ class Grid:
@staticmethod @staticmethod
def load_Neper(fname: Union[str, Path]) -> "Grid": def load_Neper(fname: Union[str, Path]) -> 'Grid':
""" """
Load from Neper VTK file. Load from Neper VTK file.
@ -279,7 +279,7 @@ class Grid:
Grid-based geometry from file. Grid-based geometry from file.
""" """
v = VTK.load(fname,'vtkImageData') v = VTK.load(fname,'ImageData')
cells = np.array(v.vtk_data.GetDimensions())-1 cells = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
@ -292,7 +292,7 @@ class Grid:
def load_DREAM3D(fname: Union[str, Path], def load_DREAM3D(fname: Union[str, Path],
feature_IDs: str = None, cell_data: str = None, feature_IDs: str = None, cell_data: str = None,
phases: str = 'Phases', Euler_angles: str = 'EulerAngles', phases: str = 'Phases', Euler_angles: str = 'EulerAngles',
base_group: str = None) -> "Grid": base_group: str = None) -> 'Grid':
""" """
Load DREAM.3D (HDF5) file. Load DREAM.3D (HDF5) file.
@ -354,7 +354,7 @@ class Grid:
@staticmethod @staticmethod
def from_table(table: Table, def from_table(table: Table,
coordinates: str, coordinates: str,
labels: Union[str, Sequence[str]]) -> "Grid": labels: Union[str, Sequence[str]]) -> 'Grid':
""" """
Create grid from ASCII table. Create grid from ASCII table.
@ -422,13 +422,14 @@ class Grid:
Grid-based geometry from tessellation. Grid-based geometry from tessellation.
""" """
weights_p: FloatSequence
if periodic: if periodic:
weights_p = np.tile(weights,27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3) weights_p = np.tile(weights,27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.]))) seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]]))) seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
else: else:
weights_p = weights weights_p = np.array(weights,float)
seeds_p = seeds seeds_p = seeds
coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3) coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3)
@ -452,7 +453,7 @@ class Grid:
size: FloatSequence, size: FloatSequence,
seeds: np.ndarray, seeds: np.ndarray,
material: IntSequence = None, material: IntSequence = None,
periodic: bool = True) -> "Grid": periodic: bool = True) -> 'Grid':
""" """
Create grid from Voronoi tessellation. Create grid from Voronoi tessellation.
@ -538,7 +539,7 @@ class Grid:
surface: str, surface: str,
threshold: float = 0.0, threshold: float = 0.0,
periods: int = 1, periods: int = 1,
materials: IntSequence = (0,1)) -> "Grid": materials: IntSequence = (0,1)) -> 'Grid':
""" """
Create grid from definition of triply periodic minimal surface. Create grid from definition of triply periodic minimal surface.
@ -674,7 +675,7 @@ class Grid:
def show(self) -> None: def show(self) -> None:
"""Show on screen.""" """Show on screen."""
VTK.from_rectilinear_grid(self.cells,self.size,self.origin).show() VTK.from_image_data(self.cells,self.size,self.origin).show()
def add_primitive(self, def add_primitive(self,
@ -684,7 +685,7 @@ class Grid:
fill: int = None, fill: int = None,
R: Rotation = Rotation(), R: Rotation = Rotation(),
inverse: bool = False, inverse: bool = False,
periodic: bool = True) -> "Grid": periodic: bool = True) -> 'Grid':
""" """
Insert a primitive geometric object at a given position. Insert a primitive geometric object at a given position.
@ -769,7 +770,7 @@ class Grid:
) )
def mirror(self, directions: Sequence[str], reflect: bool = False) -> "Grid": def mirror(self, directions: Sequence[str], reflect: bool = False) -> 'Grid':
""" """
Mirror grid along given directions. Mirror grid along given directions.
@ -821,7 +822,7 @@ class Grid:
) )
def flip(self, directions: Sequence[str]) -> "Grid": def flip(self, directions: Union[Literal['x', 'y', 'z'], Sequence[Literal['x', 'y', 'z']]]) -> 'Grid':
""" """
Flip grid along given directions. Flip grid along given directions.
@ -841,7 +842,8 @@ class Grid:
if not set(directions).issubset(valid): if not set(directions).issubset(valid):
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified') raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
mat = np.flip(self.material, (valid.index(d) for d in directions if d in valid))
mat = np.flip(self.material, [valid.index(d) for d in directions if d in valid])
return Grid(material = mat, return Grid(material = mat,
size = self.size, size = self.size,
@ -850,7 +852,7 @@ class Grid:
) )
def scale(self, cells: IntSequence, periodic: bool = True) -> "Grid": def scale(self, cells: IntSequence, periodic: bool = True) -> 'Grid':
""" """
Scale grid to new cells. Scale grid to new cells.
@ -897,7 +899,7 @@ class Grid:
def clean(self, def clean(self,
stencil: int = 3, stencil: int = 3,
selection: IntSequence = None, selection: IntSequence = None,
periodic: bool = True) -> "Grid": periodic: bool = True) -> 'Grid':
""" """
Smooth grid by selecting most frequent material index within given stencil at each location. Smooth grid by selecting most frequent material index within given stencil at each location.
@ -937,7 +939,7 @@ class Grid:
) )
def renumber(self) -> "Grid": def renumber(self) -> 'Grid':
""" """
Renumber sorted material indices as 0,...,N-1. Renumber sorted material indices as 0,...,N-1.
@ -956,7 +958,7 @@ class Grid:
) )
def rotate(self, R: Rotation, fill: int = None) -> "Grid": def rotate(self, R: Rotation, fill: int = None) -> 'Grid':
""" """
Rotate grid (pad if required). Rotate grid (pad if required).
@ -973,14 +975,13 @@ class Grid:
Updated grid-based geometry. Updated grid-based geometry.
""" """
if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if isinstance(fill,float) or self.material.dtype in np.sctypes['float'] else int
material = self.material material = self.material
# These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'') # These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
# see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf # see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
for angle,axes in zip(R.as_Euler_angles(degrees=True)[::-1], [(0,1),(1,2),(0,1)]): for angle,axes in zip(R.as_Euler_angles(degrees=True)[::-1], [(0,1),(1,2),(0,1)]):
material_temp = ndimage.rotate(material,angle,axes,order=0,prefilter=False,output=dtype,cval=fill) material_temp = ndimage.rotate(material,angle,axes,order=0,prefilter=False,
output=self.material.dtype,
cval=np.nanmax(self.material) + 1 if fill is None else fill)
# avoid scipy interpolation errors for rotations close to multiples of 90° # avoid scipy interpolation errors for rotations close to multiples of 90°
material = material_temp if np.prod(material_temp.shape) != np.prod(material.shape) else \ material = material_temp if np.prod(material_temp.shape) != np.prod(material.shape) else \
np.rot90(material,k=np.rint(angle/90.).astype(int),axes=axes) np.rot90(material,k=np.rint(angle/90.).astype(int),axes=axes)
@ -997,7 +998,7 @@ class Grid:
def canvas(self, def canvas(self,
cells: IntSequence = None, cells: IntSequence = None,
offset: IntSequence = None, offset: IntSequence = None,
fill: int = None) -> "Grid": fill: int = None) -> 'Grid':
""" """
Crop or enlarge/pad grid. Crop or enlarge/pad grid.
@ -1031,10 +1032,8 @@ class Grid:
""" """
offset_ = np.array(offset,int) if offset is not None else np.zeros(3,int) offset_ = np.array(offset,int) if offset is not None else np.zeros(3,int)
cells_ = np.array(cells,int) if cells is not None else self.cells cells_ = np.array(cells,int) if cells is not None else self.cells
if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if int(fill) != fill or self.material.dtype in np.sctypes['float'] else int
canvas = np.full(cells_,fill,dtype) canvas = np.full(cells_,np.nanmax(self.material) + 1 if fill is None else fill,self.material.dtype)
LL = np.clip( offset_, 0,np.minimum(self.cells, cells_+offset_)) LL = np.clip( offset_, 0,np.minimum(self.cells, cells_+offset_))
UR = np.clip( offset_+cells_, 0,np.minimum(self.cells, cells_+offset_)) UR = np.clip( offset_+cells_, 0,np.minimum(self.cells, cells_+offset_))
@ -1050,7 +1049,7 @@ class Grid:
) )
def substitute(self, from_material: IntSequence, to_material: IntSequence) -> "Grid": def substitute(self, from_material: IntSequence, to_material: IntSequence) -> 'Grid':
""" """
Substitute material indices. Substitute material indices.
@ -1067,20 +1066,18 @@ class Grid:
Updated grid-based geometry. Updated grid-based geometry.
""" """
def mp(entry, mapper): material = self.material.copy()
return mapper[entry] if entry in mapper else entry for f,t in zip(from_material,to_material): # ToDo Python 3.10 has strict mode for zip
material[self.material==f] = t
mp = np.vectorize(mp) return Grid(material = material,
mapper = dict(zip(from_material,to_material))
return Grid(material = mp(self.material,mapper).reshape(self.cells),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Grid','substitute')], comments = self.comments+[util.execution_stamp('Grid','substitute')],
) )
def sort(self) -> "Grid": def sort(self) -> 'Grid':
""" """
Sort material indices such that min(material) is located at (0,0,0). Sort material indices such that min(material) is located at (0,0,0).
@ -1106,7 +1103,7 @@ class Grid:
vicinity: int = 1, vicinity: int = 1,
offset: int = None, offset: int = None,
trigger: IntSequence = [], trigger: IntSequence = [],
periodic: bool = True) -> "Grid": periodic: bool = True) -> 'Grid':
""" """
Offset material index of points in the vicinity of xxx. Offset material index of points in the vicinity of xxx.

2
python/damask/_result.py
View File

@ -981,7 +981,7 @@ class Result:
t = 'tensor' t = 'tensor'
if o is None: o = 'fro' if o is None: o = 'fro'
else: else:
raise ValueError(f'invalid norm order {ord}') raise ValueError(f'invalid shape of {x["label"]}')
return { return {
'data': np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True), 'data': np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True),

85
python/damask/_table.py
View File

@ -1,15 +1,18 @@
import re import re
import copy import copy
from pathlib import Path
from typing import Union, Optional, Tuple, List
import pandas as pd import pandas as pd
import numpy as np import numpy as np
from ._typehints import FileHandle
from . import util from . import util
class Table: class Table:
"""Manipulate multi-dimensional spreadsheet-like data.""" """Manipulate multi-dimensional spreadsheet-like data."""
def __init__(self,data,shapes,comments=None): def __init__(self, data: np.ndarray, shapes: dict, comments: Optional[Union[str, list]] = None):
""" """
New spreadsheet. New spreadsheet.
@ -30,7 +33,7 @@ class Table:
self._relabel('uniform') self._relabel('uniform')
def __repr__(self): def __repr__(self) -> str:
"""Brief overview.""" """Brief overview."""
self._relabel('shapes') self._relabel('shapes')
data_repr = self.data.__repr__() data_repr = self.data.__repr__()
@ -38,7 +41,7 @@ class Table:
return '\n'.join(['# '+c for c in self.comments])+'\n'+data_repr return '\n'.join(['# '+c for c in self.comments])+'\n'+data_repr
def __getitem__(self,item): def __getitem__(self, item: Union[slice, Tuple[slice, ...]]) -> 'Table':
""" """
Slice the Table according to item. Slice the Table according to item.
@ -85,19 +88,19 @@ class Table:
comments=self.comments) comments=self.comments)
def __len__(self): def __len__(self) -> int:
"""Number of rows.""" """Number of rows."""
return len(self.data) return len(self.data)
def __copy__(self): def __copy__(self) -> 'Table':
"""Create deep copy.""" """Create deep copy."""
return copy.deepcopy(self) return copy.deepcopy(self)
copy = __copy__ copy = __copy__
def _label(self,what,how): def _label(self, what: Union[str, List[str]], how: str) -> List[str]:
""" """
Expand labels according to data shape. Expand labels according to data shape.
@ -105,7 +108,7 @@ class Table:
---------- ----------
what : str or list what : str or list
Labels to expand. Labels to expand.
how : str how : {'uniform, 'shapes', 'linear'}
Mode of labeling. Mode of labeling.
'uniform' ==> v v v 'uniform' ==> v v v
'shapes' ==> 3:v v v 'shapes' ==> 3:v v v
@ -128,30 +131,34 @@ class Table:
return labels return labels
def _relabel(self,how): def _relabel(self, how: str):
""" """
Modify labeling of data in-place. Modify labeling of data in-place.
Parameters Parameters
---------- ----------
how : str how : {'uniform, 'shapes', 'linear'}
Mode of labeling. Mode of labeling.
'uniform' ==> v v v 'uniform' ==> v v v
'shapes' ==> 3:v v v 'shapes' ==> 3:v v v
'linear' ==> 1_v 2_v 3_v 'linear' ==> 1_v 2_v 3_v
""" """
self.data.columns = self._label(self.shapes,how) self.data.columns = self._label(self.shapes,how) #type: ignore
def _add_comment(self,label,shape,info): def _add_comment(self, label: str, shape: Tuple[int, ...], info: Optional[str]):
if info is not None: if info is not None:
specific = f'{label}{" "+str(shape) if np.prod(shape,dtype=int) > 1 else ""}: {info}' specific = f'{label}{" "+str(shape) if np.prod(shape,dtype=int) > 1 else ""}: {info}'
general = util.execution_stamp('Table') general = util.execution_stamp('Table')
self.comments.append(f'{specific} / {general}') self.comments.append(f'{specific} / {general}')
def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): def isclose(self,
other: 'Table',
rtol: float = 1e-5,
atol: float = 1e-8,
equal_nan: bool = True) -> np.ndarray:
""" """
Report where values are approximately equal to corresponding ones of other Table. Report where values are approximately equal to corresponding ones of other Table.
@ -179,7 +186,11 @@ class Table:
equal_nan=equal_nan) equal_nan=equal_nan)
def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): def allclose(self,
other: 'Table',
rtol: float = 1e-5,
atol: float = 1e-8,
equal_nan: bool = True) -> bool:
""" """
Test whether all values are approximately equal to corresponding ones of other Table. Test whether all values are approximately equal to corresponding ones of other Table.
@ -208,7 +219,7 @@ class Table:
@staticmethod @staticmethod
def load(fname): def load(fname: FileHandle) -> 'Table':
""" """
Load from ASCII table file. Load from ASCII table file.
@ -229,10 +240,7 @@ class Table:
Table data from file. Table data from file.
""" """
try: f = open(fname) if isinstance(fname, (str, Path)) else fname
f = open(fname)
except TypeError:
f = fname
f.seek(0) f.seek(0)
comments = [] comments = []
@ -261,7 +269,7 @@ class Table:
@staticmethod @staticmethod
def load_ang(fname): def load_ang(fname: FileHandle) -> 'Table':
""" """
Load from ang file. Load from ang file.
@ -286,10 +294,7 @@ class Table:
Table data from file. Table data from file.
""" """
try: f = open(fname) if isinstance(fname, (str, Path)) else fname
f = open(fname)
except TypeError:
f = fname
f.seek(0) f.seek(0)
content = f.readlines() content = f.readlines()
@ -312,11 +317,11 @@ class Table:
@property @property
def labels(self): def labels(self) -> List[Tuple[int, ...]]:
return list(self.shapes) return list(self.shapes)
def get(self,label): def get(self, label: str) -> np.ndarray:
""" """
Get column data. Get column data.
@ -336,7 +341,7 @@ class Table:
return data.astype(type(data.flatten()[0])) return data.astype(type(data.flatten()[0]))
def set(self,label,data,info=None): def set(self, label: str, data: np.ndarray, info: str = None) -> 'Table':
""" """
Set column data. Set column data.
@ -356,7 +361,7 @@ class Table:
""" """
dup = self.copy() dup = self.copy()
dup._add_comment(label,data.shape[1:],info) dup._add_comment(label, data.shape[1:], info)
m = re.match(r'(.*)\[((\d+,)*(\d+))\]',label) m = re.match(r'(.*)\[((\d+,)*(\d+))\]',label)
if m: if m:
key = m.group(1) key = m.group(1)
@ -369,7 +374,7 @@ class Table:
return dup return dup
def add(self,label,data,info=None): def add(self, label: str, data: np.ndarray, info: str = None) -> 'Table':
""" """
Add column data. Add column data.
@ -401,7 +406,7 @@ class Table:
return dup return dup
def delete(self,label): def delete(self, label: str) -> 'Table':
""" """
Delete column data. Delete column data.
@ -422,7 +427,7 @@ class Table:
return dup return dup
def rename(self,old,new,info=None): def rename(self, old: Union[str, List[str]], new: Union[str, List[str]], info: str = None) -> 'Table':
""" """
Rename column data. Rename column data.
@ -448,7 +453,7 @@ class Table:
return dup return dup
def sort_by(self,labels,ascending=True): def sort_by(self, labels: Union[str, List[str]], ascending: Union[bool, List[bool]] = True) -> 'Table':
""" """
Sort table by values of given labels. Sort table by values of given labels.
@ -481,7 +486,7 @@ class Table:
return dup return dup
def append(self,other): def append(self, other: 'Table') -> 'Table':
""" """
Append other table vertically (similar to numpy.vstack). Append other table vertically (similar to numpy.vstack).
@ -506,7 +511,7 @@ class Table:
return dup return dup
def join(self,other): def join(self, other: 'Table') -> 'Table':
""" """
Append other table horizontally (similar to numpy.hstack). Append other table horizontally (similar to numpy.hstack).
@ -533,7 +538,7 @@ class Table:
return dup return dup
def save(self,fname): def save(self, fname: FileHandle):
""" """
Save as plain text file. Save as plain text file.
@ -543,9 +548,8 @@ class Table:
Filename or file for writing. Filename or file for writing.
""" """
seen = set()
labels = [] labels = []
for l in [x for x in self.data.columns if not (x in seen or seen.add(x))]: for l in list(dict.fromkeys(self.data.columns)):
if self.shapes[l] == (1,): if self.shapes[l] == (1,):
labels.append(f'{l}') labels.append(f'{l}')
elif len(self.shapes[l]) == 1: elif len(self.shapes[l]) == 1:
@ -555,10 +559,7 @@ class Table:
labels += [f'{util.srepr(self.shapes[l],"x")}:{i+1}_{l}' \ labels += [f'{util.srepr(self.shapes[l],"x")}:{i+1}_{l}' \
for i in range(np.prod(self.shapes[l]))] for i in range(np.prod(self.shapes[l]))]
try: f = open(fname,'w',newline='\n') if isinstance(fname, (str, Path)) else fname
fhandle = open(fname,'w',newline='\n')
except TypeError:
fhandle = fname
fhandle.write('\n'.join([f'# {c}' for c in self.comments] + [' '.join(labels)])+'\n') f.write('\n'.join([f'# {c}' for c in self.comments] + [' '.join(labels)])+'\n')
self.data.to_csv(fhandle,sep=' ',na_rep='nan',index=False,header=False) self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False)

52
python/damask/_vtk.py
View File

@ -2,6 +2,7 @@ import os
import warnings import warnings
import multiprocessing as mp import multiprocessing as mp
from pathlib import Path from pathlib import Path
from typing import Union, Literal, List
import numpy as np import numpy as np
import vtk import vtk
@ -9,6 +10,7 @@ from vtk.util.numpy_support import numpy_to_vtk as np_to_vtk
from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray
from vtk.util.numpy_support import vtk_to_numpy as vtk_to_np from vtk.util.numpy_support import vtk_to_numpy as vtk_to_np
from ._typehints import FloatSequence, IntSequence
from . import util from . import util
from . import Table from . import Table
@ -20,7 +22,7 @@ class VTK:
High-level interface to VTK. High-level interface to VTK.
""" """
def __init__(self,vtk_data): def __init__(self, vtk_data: vtk.vtkDataSet):
""" """
New spatial visualization. New spatial visualization.
@ -36,7 +38,7 @@ class VTK:
@staticmethod @staticmethod
def from_image_data(cells,size,origin=np.zeros(3)): def from_image_data(cells: IntSequence, size: FloatSequence, origin: FloatSequence = np.zeros(3)) -> 'VTK':
""" """
Create VTK of type vtk.vtkImageData. Create VTK of type vtk.vtkImageData.
@ -60,13 +62,13 @@ class VTK:
vtk_data = vtk.vtkImageData() vtk_data = vtk.vtkImageData()
vtk_data.SetDimensions(*(np.array(cells)+1)) vtk_data.SetDimensions(*(np.array(cells)+1))
vtk_data.SetOrigin(*(np.array(origin))) vtk_data.SetOrigin(*(np.array(origin)))
vtk_data.SetSpacing(*(size/cells)) vtk_data.SetSpacing(*(np.array(size)/np.array(cells)))
return VTK(vtk_data) return VTK(vtk_data)
@staticmethod @staticmethod
def from_rectilinear_grid(grid,size,origin=np.zeros(3)): def from_rectilinear_grid(grid: np.ndarray, size: FloatSequence, origin: FloatSequence = np.zeros(3)) -> 'VTK':
""" """
Create VTK of type vtk.vtkRectilinearGrid. Create VTK of type vtk.vtkRectilinearGrid.
@ -98,7 +100,7 @@ class VTK:
@staticmethod @staticmethod
def from_unstructured_grid(nodes,connectivity,cell_type): def from_unstructured_grid(nodes: np.ndarray, connectivity: np.ndarray, cell_type: str) -> 'VTK':
""" """
Create VTK of type vtk.vtkUnstructuredGrid. Create VTK of type vtk.vtkUnstructuredGrid.
@ -138,7 +140,7 @@ class VTK:
@staticmethod @staticmethod
def from_poly_data(points): def from_poly_data(points: np.ndarray) -> 'VTK':
""" """
Create VTK of type vtk.polyData. Create VTK of type vtk.polyData.
@ -172,15 +174,17 @@ class VTK:
@staticmethod @staticmethod
def load(fname,dataset_type=None): def load(fname: Union[str, Path],
dataset_type: Literal['ImageData', 'UnstructuredGrid', 'PolyData'] = None) -> 'VTK':
""" """
Load from VTK file. Load from VTK file.
Parameters Parameters
---------- ----------
fname : str or pathlib.Path fname : str or pathlib.Path
Filename for reading. Valid extensions are .vti, .vtr, .vtu, .vtp, and .vtk. Filename for reading.
dataset_type : {'vtkImageData', ''vtkRectilinearGrid', 'vtkUnstructuredGrid', 'vtkPolyData'}, optional Valid extensions are .vti, .vtr, .vtu, .vtp, and .vtk.
dataset_type : {'ImageData', 'UnstructuredGrid', 'PolyData'}, optional
Name of the vtk.vtkDataSet subclass when opening a .vtk file. Name of the vtk.vtkDataSet subclass when opening a .vtk file.
Returns Returns
@ -234,7 +238,7 @@ class VTK:
def _write(writer): def _write(writer):
"""Wrapper for parallel writing.""" """Wrapper for parallel writing."""
writer.Write() writer.Write()
def save(self,fname,parallel=True,compress=True): def save(self, fname: Union[str, Path], parallel: bool = True, compress: bool = True):
""" """
Save as VTK file. Save as VTK file.
@ -280,7 +284,7 @@ class VTK:
# Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data # Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data
# Needs support for damask.Table # Needs support for damask.Table
def add(self,data,label=None): def add(self, data: Union[np.ndarray, np.ma.MaskedArray], label: str = None):
""" """
Add data to either cells or points. Add data to either cells or points.
@ -327,7 +331,7 @@ class VTK:
raise TypeError raise TypeError
def get(self,label): def get(self, label: str) -> np.ndarray:
""" """
Get either cell or point data. Get either cell or point data.
@ -369,7 +373,7 @@ class VTK:
raise ValueError(f'Array "{label}" not found.') raise ValueError(f'Array "{label}" not found.')
def get_comments(self): def get_comments(self) -> List[str]:
"""Return the comments.""" """Return the comments."""
fielddata = self.vtk_data.GetFieldData() fielddata = self.vtk_data.GetFieldData()
for a in range(fielddata.GetNumberOfArrays()): for a in range(fielddata.GetNumberOfArrays()):
@ -379,7 +383,7 @@ class VTK:
return [] return []
def set_comments(self,comments): def set_comments(self, comments: Union[str, List[str]]):
""" """
Set comments. Set comments.
@ -396,7 +400,7 @@ class VTK:
self.vtk_data.GetFieldData().AddArray(s) self.vtk_data.GetFieldData().AddArray(s)
def add_comments(self,comments): def add_comments(self, comments: Union[str, List[str]]):
""" """
Add comments. Add comments.
@ -409,7 +413,7 @@ class VTK:
self.set_comments(self.get_comments() + ([comments] if isinstance(comments,str) else comments)) self.set_comments(self.get_comments() + ([comments] if isinstance(comments,str) else comments))
def __repr__(self): def __repr__(self) -> str:
"""ASCII representation of the VTK data.""" """ASCII representation of the VTK data."""
writer = vtk.vtkDataSetWriter() writer = vtk.vtkDataSetWriter()
writer.SetHeader(f'# {util.execution_stamp("VTK")}') writer.SetHeader(f'# {util.execution_stamp("VTK")}')
@ -419,7 +423,7 @@ class VTK:
return writer.GetOutputString() return writer.GetOutputString()
def show(self) -> None: def show(self):
""" """
Render. Render.
@ -442,22 +446,24 @@ class VTK:
mapper = vtk.vtkDataSetMapper() mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(self.vtk_data) mapper.SetInputData(self.vtk_data)
actor = vtk.vtkActor() actor = vtk.vtkActor()
actor.SetMapper(mapper) actor.SetMapper(mapper)
actor.GetProperty().SetColor(230/255,150/255,68/255)
ren = vtk.vtkRenderer() ren = vtk.vtkRenderer()
ren.AddActor(actor)
ren.SetBackground(67/255,128/255,208/255)
window = vtk.vtkRenderWindow() window = vtk.vtkRenderWindow()
window.AddRenderer(ren) window.AddRenderer(ren)
ren.AddActor(actor)
ren.SetBackground(0.2,0.2,0.2)
window.SetSize(width,height) window.SetSize(width,height)
window.SetWindowName(util.execution_stamp('VTK','show'))
iren = vtk.vtkRenderWindowInteractor() iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(window) iren.SetRenderWindow(window)
if os.name == 'posix' and 'DISPLAY' not in os.environ:
iren.Initialize() print('Found no rendering device')
else:
window.Render() window.Render()
iren.Start() iren.Start()

2
python/damask/seeds.py
View File

@ -79,7 +79,7 @@ def from_Poisson_disc(size: _FloatSequence, N_seeds: int, N_candidates: int, dis
s = 1 s = 1
i = 0 i = 0
progress = _util._ProgressBar(N_seeds+1,'',50) progress = _util.ProgressBar(N_seeds+1,'',50)
while s < N_seeds: while s < N_seeds:
i += 1 i += 1
candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3)) candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3))

167
python/damask/util.py
View File

@ -7,12 +7,16 @@ import subprocess
import shlex import shlex
import re import re
import fractions import fractions
import collections.abc as abc
from functools import reduce from functools import reduce
from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal, Optional
from pathlib import Path
import numpy as np import numpy as np
import h5py import h5py
from . import version from . import version
from ._typehints import FloatSequence
# limit visibility # limit visibility
__all__=[ __all__=[
@ -50,16 +54,16 @@ _colors = {
#################################################################################################### ####################################################################################################
# Functions # Functions
#################################################################################################### ####################################################################################################
def srepr(arg,glue = '\n'): def srepr(msg, glue: str = '\n') -> str:
r""" r"""
Join items with glue string. Join items with glue string.
Parameters Parameters
---------- ----------
arg : iterable msg : object with __repr__ or sequence of objects with __repr__
Items to join. Items to join.
glue : str, optional glue : str, optional
Glue used for joining operation. Defaults to \n. Glue used for joining operation. Defaults to '\n'.
Returns Returns
------- -------
@ -67,21 +71,21 @@ def srepr(arg,glue = '\n'):
String representation of the joined items. String representation of the joined items.
""" """
if (not hasattr(arg, 'strip') and if (not hasattr(msg, 'strip') and
(hasattr(arg, '__getitem__') or (hasattr(msg, '__getitem__') or
hasattr(arg, '__iter__'))): hasattr(msg, '__iter__'))):
return glue.join(str(x) for x in arg) return glue.join(str(x) for x in msg)
else: else:
return arg if isinstance(arg,str) else repr(arg) return msg if isinstance(msg,str) else repr(msg)
def emph(what): def emph(msg) -> str:
""" """
Format with emphasis. Format with emphasis.
Parameters Parameters
---------- ----------
what : object with __repr__ or iterable of objects with __repr__. msg : object with __repr__ or sequence of objects with __repr__
Message to format. Message to format.
Returns Returns
@ -90,15 +94,15 @@ def emph(what):
Formatted string representation of the joined items. Formatted string representation of the joined items.
""" """
return _colors['bold']+srepr(what)+_colors['end_color'] return _colors['bold']+srepr(msg)+_colors['end_color']
def deemph(what): def deemph(msg) -> str:
""" """
Format with deemphasis. Format with deemphasis.
Parameters Parameters
---------- ----------
what : object with __repr__ or iterable of objects with __repr__. msg : object with __repr__ or sequence of objects with __repr__
Message to format. Message to format.
Returns Returns
@ -107,15 +111,15 @@ def deemph(what):
Formatted string representation of the joined items. Formatted string representation of the joined items.
""" """
return _colors['dim']+srepr(what)+_colors['end_color'] return _colors['dim']+srepr(msg)+_colors['end_color']
def warn(what): def warn(msg) -> str:
""" """
Format for warning. Format for warning.
Parameters Parameters
---------- ----------
what : object with __repr__ or iterable of objects with __repr__. msg : object with __repr__ or sequence of objects with __repr__
Message to format. Message to format.
Returns Returns
@ -124,15 +128,15 @@ def warn(what):
Formatted string representation of the joined items. Formatted string representation of the joined items.
""" """
return _colors['warning']+emph(what)+_colors['end_color'] return _colors['warning']+emph(msg)+_colors['end_color']
def strikeout(what): def strikeout(msg) -> str:
""" """
Format as strikeout. Format as strikeout.
Parameters Parameters
---------- ----------
what : object with __repr__ or iterable of objects with __repr__. msg : object with __repr__ or iterable of objects with __repr__
Message to format. Message to format.
Returns Returns
@ -141,10 +145,10 @@ def strikeout(what):
Formatted string representation of the joined items. Formatted string representation of the joined items.
""" """
return _colors['crossout']+srepr(what)+_colors['end_color'] return _colors['crossout']+srepr(msg)+_colors['end_color']
def run(cmd,wd='./',env=None,timeout=None): def run(cmd: str, wd: str = './', env: Dict[str, str] = None, timeout: int = None) -> Tuple[str, str]:
""" """
Run a command. Run a command.
@ -153,7 +157,7 @@ def run(cmd,wd='./',env=None,timeout=None):
cmd : str cmd : str
Command to be executed. Command to be executed.
wd : str, optional wd : str, optional
Working directory of process. Defaults to ./ . Working directory of process. Defaults to './'.
env : dict, optional env : dict, optional
Environment for execution. Environment for execution.
timeout : integer, optional timeout : integer, optional
@ -185,7 +189,7 @@ def run(cmd,wd='./',env=None,timeout=None):
execute = run execute = run
def natural_sort(key): def natural_sort(key: str) -> List[Union[int, str]]:
""" """
Natural sort. Natural sort.
@ -200,7 +204,10 @@ def natural_sort(key):
return [ convert(c) for c in re.split('([0-9]+)', key) ] return [ convert(c) for c in re.split('([0-9]+)', key) ]
def show_progress(iterable,N_iter=None,prefix='',bar_length=50): def show_progress(iterable: Iterable,
N_iter: int = None,
prefix: str = '',
bar_length: int = 50) -> Any:
""" """
Decorate a loop with a progress bar. Decorate a loop with a progress bar.
@ -208,39 +215,49 @@ def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
Parameters Parameters
---------- ----------
iterable : iterable or function with yield statement iterable : iterable
Iterable (or function with yield statement) to be decorated. Iterable to be decorated.
N_iter : int, optional N_iter : int, optional
Total number of iterations. Required unless obtainable as len(iterable). Total number of iterations. Required if iterable is not a sequence.
prefix : str, optional prefix : str, optional
Prefix string. Prefix string.
bar_length : int, optional bar_length : int, optional
Length of progress bar in characters. Defaults to 50. Length of progress bar in characters. Defaults to 50.
""" """
if N_iter in [0,1] or (hasattr(iterable,'__len__') and len(iterable) <= 1): if isinstance(iterable,abc.Sequence):
if N_iter is None:
N = len(iterable)
else:
raise ValueError('N_iter given for sequence')
else:
if N_iter is None:
raise ValueError('N_iter not given')
else:
N = N_iter
if N <= 1:
for item in iterable: for item in iterable:
yield item yield item
else: else:
status = _ProgressBar(N_iter if N_iter is not None else len(iterable),prefix,bar_length) status = ProgressBar(N,prefix,bar_length)
for i,item in enumerate(iterable): for i,item in enumerate(iterable):
yield item yield item
status.update(i) status.update(i)
def scale_to_coprime(v): def scale_to_coprime(v: FloatSequence) -> np.ndarray:
""" """
Scale vector to co-prime (relatively prime) integers. Scale vector to co-prime (relatively prime) integers.
Parameters Parameters
---------- ----------
v : numpy.ndarray of shape (:) v : sequence of float, len (:)
Vector to scale. Vector to scale.
Returns Returns
------- -------
m : numpy.ndarray of shape (:) m : numpy.ndarray, shape (:)
Vector scaled to co-prime numbers. Vector scaled to co-prime numbers.
""" """
@ -257,17 +274,21 @@ def scale_to_coprime(v):
except AttributeError: except AttributeError:
return a * b // np.gcd(a, b) return a * b // np.gcd(a, b)
m = (np.array(v) * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v)) ** 0.5).astype(int) v_ = np.array(v)
m = (v_ * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v_))**0.5).astype(int)
m = m//reduce(np.gcd,m) m = m//reduce(np.gcd,m)
with np.errstate(invalid='ignore'): with np.errstate(invalid='ignore'):
if not np.allclose(np.ma.masked_invalid(v/m),v[np.argmax(abs(v))]/m[np.argmax(abs(v))]): if not np.allclose(np.ma.masked_invalid(v_/m),v_[np.argmax(abs(v_))]/m[np.argmax(abs(v_))]):
raise ValueError(f'Invalid result {m} for input {v}. Insufficient precision?') raise ValueError(f'Invalid result {m} for input {v_}. Insufficient precision?')
return m return m
def project_equal_angle(vector,direction='z',normalize=True,keepdims=False): def project_equal_angle(vector: np.ndarray,
direction: Literal['x', 'y', 'z'] = 'z',
normalize: bool = True,
keepdims: bool = False) -> np.ndarray:
""" """
Apply equal-angle projection to vector. Apply equal-angle projection to vector.
@ -275,9 +296,8 @@ def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
---------- ----------
vector : numpy.ndarray, shape (...,3) vector : numpy.ndarray, shape (...,3)
Vector coordinates to be projected. Vector coordinates to be projected.
direction : str direction : {'x', 'y', 'z'}
Projection direction 'x', 'y', or 'z'. Projection direction. Defaults to 'z'.
Defaults to 'z'.
normalize : bool normalize : bool
Ensure unit length of input vector. Defaults to True. Ensure unit length of input vector. Defaults to True.
keepdims : bool keepdims : bool
@ -309,7 +329,10 @@ def project_equal_angle(vector,direction='z',normalize=True,keepdims=False):
return np.roll(np.block([v[...,:2]/(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]), return np.roll(np.block([v[...,:2]/(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
-shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2] -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
def project_equal_area(vector,direction='z',normalize=True,keepdims=False): def project_equal_area(vector: np.ndarray,
direction: Literal['x', 'y', 'z'] = 'z',
normalize: bool = True,
keepdims: bool = False) -> np.ndarray:
""" """
Apply equal-area projection to vector. Apply equal-area projection to vector.
@ -317,9 +340,8 @@ def project_equal_area(vector,direction='z',normalize=True,keepdims=False):
---------- ----------
vector : numpy.ndarray, shape (...,3) vector : numpy.ndarray, shape (...,3)
Vector coordinates to be projected. Vector coordinates to be projected.
direction : str direction : {'x', 'y', 'z'}
Projection direction 'x', 'y', or 'z'. Projection direction. Defaults to 'z'.
Defaults to 'z'.
normalize : bool normalize : bool
Ensure unit length of input vector. Defaults to True. Ensure unit length of input vector. Defaults to True.
keepdims : bool keepdims : bool
@ -351,15 +373,14 @@ def project_equal_area(vector,direction='z',normalize=True,keepdims=False):
return np.roll(np.block([v[...,:2]/np.sqrt(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]), return np.roll(np.block([v[...,:2]/np.sqrt(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
-shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2] -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
def execution_stamp(class_name: str, function_name: str = None) -> str:
def execution_stamp(class_name,function_name=None):
"""Timestamp the execution of a (function within a) class.""" """Timestamp the execution of a (function within a) class."""
now = datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z') now = datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z')
_function_name = '' if function_name is None else f'.{function_name}' _function_name = '' if function_name is None else f'.{function_name}'
return f'damask.{class_name}{_function_name} v{version} ({now})' return f'damask.{class_name}{_function_name} v{version} ({now})'
def hybrid_IA(dist,N,rng_seed=None): def hybrid_IA(dist: np.ndarray, N: int, rng_seed = None) -> np.ndarray:
""" """
Hybrid integer approximation. Hybrid integer approximation.
@ -387,7 +408,10 @@ def hybrid_IA(dist,N,rng_seed=None):
return np.repeat(np.arange(len(dist)),repeats)[np.random.default_rng(rng_seed).permutation(N_inv_samples)[:N]] return np.repeat(np.arange(len(dist)),repeats)[np.random.default_rng(rng_seed).permutation(N_inv_samples)[:N]]
def shapeshifter(fro,to,mode='left',keep_ones=False): def shapeshifter(fro: Tuple[int, ...],
to: Tuple[int, ...],
mode: Literal['left','right'] = 'left',
keep_ones: bool = False) -> Tuple[Optional[int], ...]:
""" """
Return dimensions that reshape 'fro' to become broadcastable to 'to'. Return dimensions that reshape 'fro' to become broadcastable to 'to'.
@ -398,9 +422,9 @@ def shapeshifter(fro,to,mode='left',keep_ones=False):
to : tuple to : tuple
Target shape of array after broadcasting. Target shape of array after broadcasting.
len(to) cannot be less than len(fro). len(to) cannot be less than len(fro).
mode : str, optional mode : {'left', 'right'}, optional
Indicates whether new axes are preferably added to Indicates whether new axes are preferably added to
either 'left' or 'right' of the original shape. either left or right of the original shape.
Defaults to 'left'. Defaults to 'left'.
keep_ones : bool, optional keep_ones : bool, optional
Treat '1' in fro as literal value instead of dimensional placeholder. Treat '1' in fro as literal value instead of dimensional placeholder.
@ -434,21 +458,22 @@ def shapeshifter(fro,to,mode='left',keep_ones=False):
fro = (1,) if not len(fro) else fro fro = (1,) if not len(fro) else fro
to = (1,) if not len(to) else to to = (1,) if not len(to) else to
try: try:
grp = re.match(beg[mode] match = re.match(beg[mode]
+f',{sep[mode]}'.join(map(lambda x: f'{x}' +f',{sep[mode]}'.join(map(lambda x: f'{x}'
if x>1 or (keep_ones and len(fro)>1) else if x>1 or (keep_ones and len(fro)>1) else
'\\d+',fro)) '\\d+',fro))
+f',{end[mode]}', +f',{end[mode]}',','.join(map(str,to))+',')
','.join(map(str,to))+',').groups() assert match
except AttributeError: grp = match.groups()
except AssertionError:
raise ValueError(f'Shapes can not be shifted {fro} --> {to}') raise ValueError(f'Shapes can not be shifted {fro} --> {to}')
fill = () fill: Tuple[Optional[int], ...] = ()
for g,d in zip(grp,fro+(None,)): for g,d in zip(grp,fro+(None,)):
fill += (1,)*g.count(',')+(d,) fill += (1,)*g.count(',')+(d,)
return fill[:-1] return fill[:-1]
def shapeblender(a,b): def shapeblender(a: Tuple[int, ...], b: Tuple[int, ...]) -> Tuple[int, ...]:
""" """
Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'. Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
@ -476,7 +501,7 @@ def shapeblender(a,b):
return a + b[i:] return a + b[i:]
def extend_docstring(extra_docstring): def extend_docstring(extra_docstring: str) -> Callable:
""" """
Decorator: Append to function's docstring. Decorator: Append to function's docstring.
@ -492,7 +517,7 @@ def extend_docstring(extra_docstring):
return _decorator return _decorator
def extended_docstring(f,extra_docstring): def extended_docstring(f: Callable, extra_docstring: str) -> Callable:
""" """
Decorator: Combine another function's docstring with a given docstring. Decorator: Combine another function's docstring with a given docstring.
@ -510,7 +535,7 @@ def extended_docstring(f,extra_docstring):
return _decorator return _decorator
def DREAM3D_base_group(fname): def DREAM3D_base_group(fname: Union[str, Path]) -> str:
""" """
Determine the base group of a DREAM.3D file. Determine the base group of a DREAM.3D file.
@ -536,7 +561,7 @@ def DREAM3D_base_group(fname):
return base_group return base_group
def DREAM3D_cell_data_group(fname): def DREAM3D_cell_data_group(fname: Union[str, Path]) -> str:
""" """
Determine the cell data group of a DREAM.3D file. Determine the cell data group of a DREAM.3D file.
@ -568,18 +593,18 @@ def DREAM3D_cell_data_group(fname):
return cell_data_group return cell_data_group
def Bravais_to_Miller(*,uvtw=None,hkil=None): def Bravais_to_Miller(*, uvtw: np.ndarray = None, hkil: np.ndarray = None) -> np.ndarray:
""" """
Transform 4 Miller–Bravais indices to 3 Miller indices of crystal direction [uvw] or plane normal (hkl). Transform 4 Miller–Bravais indices to 3 Miller indices of crystal direction [uvw] or plane normal (hkl).
Parameters Parameters
---------- ----------
uvtw|hkil : numpy.ndarray of shape (...,4) uvtw|hkil : numpy.ndarray, shape (...,4)
Miller–Bravais indices of crystallographic direction [uvtw] or plane normal (hkil). Miller–Bravais indices of crystallographic direction [uvtw] or plane normal (hkil).
Returns Returns
------- -------
uvw|hkl : numpy.ndarray of shape (...,3) uvw|hkl : numpy.ndarray, shape (...,3)
Miller indices of [uvw] direction or (hkl) plane normal. Miller indices of [uvw] direction or (hkl) plane normal.
""" """
@ -595,18 +620,18 @@ def Bravais_to_Miller(*,uvtw=None,hkil=None):
return np.einsum('il,...l',basis,axis) return np.einsum('il,...l',basis,axis)
def Miller_to_Bravais(*,uvw=None,hkl=None): def Miller_to_Bravais(*, uvw: np.ndarray = None, hkl: np.ndarray = None) -> np.ndarray:
""" """
Transform 3 Miller indices to 4 Miller–Bravais indices of crystal direction [uvtw] or plane normal (hkil). Transform 3 Miller indices to 4 Miller–Bravais indices of crystal direction [uvtw] or plane normal (hkil).
Parameters Parameters
---------- ----------
uvw|hkl : numpy.ndarray of shape (...,3) uvw|hkl : numpy.ndarray, shape (...,3)
Miller indices of crystallographic direction [uvw] or plane normal (hkl). Miller indices of crystallographic direction [uvw] or plane normal (hkl).
Returns Returns
------- -------
uvtw|hkil : numpy.ndarray of shape (...,4) uvtw|hkil : numpy.ndarray, shape (...,4)
Miller–Bravais indices of [uvtw] direction or (hkil) plane normal. Miller–Bravais indices of [uvtw] direction or (hkil) plane normal.
""" """
@ -624,7 +649,7 @@ def Miller_to_Bravais(*,uvw=None,hkl=None):
return np.einsum('il,...l',basis,axis) return np.einsum('il,...l',basis,axis)
def dict_prune(d): def dict_prune(d: Dict) -> Dict:
""" """
Recursively remove empty dictionaries. Recursively remove empty dictionaries.
@ -650,7 +675,7 @@ def dict_prune(d):
return new return new
def dict_flatten(d): def dict_flatten(d: Dict) -> Dict:
""" """
Recursively remove keys of single-entry dictionaries. Recursively remove keys of single-entry dictionaries.
@ -678,14 +703,14 @@ def dict_flatten(d):
#################################################################################################### ####################################################################################################
# Classes # Classes
#################################################################################################### ####################################################################################################
class _ProgressBar: class ProgressBar:
""" """
Report progress of an interation as a status bar. Report progress of an interation as a status bar.
Works for 0-based loops, ETA is estimated by linear extrapolation. Works for 0-based loops, ETA is estimated by linear extrapolation.
""" """
def __init__(self,total,prefix,bar_length): def __init__(self, total: int, prefix: str, bar_length: int):
""" """
Set current time as basis for ETA estimation. Set current time as basis for ETA estimation.
@ -708,7 +733,7 @@ class _ProgressBar:
sys.stderr.write(f"{self.prefix} {'░'*self.bar_length} 0% ETA n/a") sys.stderr.write(f"{self.prefix} {'░'*self.bar_length} 0% ETA n/a")
sys.stderr.flush() sys.stderr.flush()
def update(self,iteration): def update(self, iteration: int) -> None:
fraction = (iteration+1) / self.total fraction = (iteration+1) / self.total
filled_length = int(self.bar_length * fraction) filled_length = int(self.bar_length * fraction)

4
python/setup.py
View File

@ -16,7 +16,7 @@ setuptools.setup(
url='https://damask.mpie.de', url='https://damask.mpie.de',
packages=setuptools.find_packages(), packages=setuptools.find_packages(),
include_package_data=True, include_package_data=True,
python_requires = '>=3.7', python_requires = '>=3.8',
install_requires = [ install_requires = [
'pandas>=0.24', # requires numpy 'pandas>=0.24', # requires numpy
'numpy>=1.17', # needed for default_rng 'numpy>=1.17', # needed for default_rng
@ -30,7 +30,7 @@ setuptools.setup(
'Intended Audience :: Science/Research', 'Intended Audience :: Science/Research',
'Topic :: Scientific/Engineering', 'Topic :: Scientific/Engineering',
'Programming Language :: Python :: 3', 'Programming Language :: Python :: 3',
'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)', 'License :: OSI Approved :: GNU Affero General Public License v3 or later (AGPLv3+)',
'Operating System :: OS Independent', 'Operating System :: OS Independent',
], ],
) )

6
python/tests/test_Colormap.py
View File

@ -81,6 +81,7 @@ class TestColormap:
assert Colormap.from_predefined('strain') == Colormap.from_predefined('strain') assert Colormap.from_predefined('strain') == Colormap.from_predefined('strain')
assert Colormap.from_predefined('strain') != Colormap.from_predefined('stress') assert Colormap.from_predefined('strain') != Colormap.from_predefined('stress')
assert Colormap.from_predefined('strain',N=128) != Colormap.from_predefined('strain',N=64) assert Colormap.from_predefined('strain',N=128) != Colormap.from_predefined('strain',N=64)
assert not Colormap.from_predefined('strain',N=128) == 1
@pytest.mark.parametrize('low,high',[((0,0,0),(1,1,1)), @pytest.mark.parametrize('low,high',[((0,0,0),(1,1,1)),
([0,0,0],[1,1,1])]) ([0,0,0],[1,1,1])])
@ -139,6 +140,11 @@ class TestColormap:
c += c c += c
assert (np.allclose(c.colors[:len(c.colors)//2],c.colors[len(c.colors)//2:])) assert (np.allclose(c.colors[:len(c.colors)//2],c.colors[len(c.colors)//2:]))
def test_mul(self):
c = o = Colormap.from_predefined('jet')
o *= 2
assert c+c == o
@pytest.mark.parametrize('N,cmap,at,result',[ @pytest.mark.parametrize('N,cmap,at,result',[
(8,'gray',0.5,[0.5,0.5,0.5]), (8,'gray',0.5,[0.5,0.5,0.5]),
(17,'gray',0.5,[0.5,0.5,0.5]), (17,'gray',0.5,[0.5,0.5,0.5]),

3
python/tests/test_Crystal.py
View File

@ -40,6 +40,9 @@ class TestCrystal:
alpha=alpha,beta=beta,gamma=gamma) alpha=alpha,beta=beta,gamma=gamma)
assert np.allclose(np.eye(3),np.einsum('ik,jk',c.basis_real,c.basis_reciprocal)) assert np.allclose(np.eye(3),np.einsum('ik,jk',c.basis_real,c.basis_reciprocal))
def test_basis_invalid(self):
with pytest.raises(KeyError):
Crystal(family='cubic').basis_real
@pytest.mark.parametrize('keyFrame,keyLattice',[('uvw','direction'),('hkl','plane'),]) @pytest.mark.parametrize('keyFrame,keyLattice',[('uvw','direction'),('hkl','plane'),])
@pytest.mark.parametrize('vector',np.array([ @pytest.mark.parametrize('vector',np.array([

1
python/tests/test_Grid.py
View File

@ -44,6 +44,7 @@ class TestGrid:
def test_equal(self,default): def test_equal(self,default):
assert default == default assert default == default
assert not default == 42
def test_repr(self,default): def test_repr(self,default):
print(default) print(default)

9
python/tests/test_Orientation.py
View File

@ -364,6 +364,11 @@ class TestOrientation:
table.save(reference) table.save(reference)
assert np.allclose(P,Table.load(reference).get('Schmid')) assert np.allclose(P,Table.load(reference).get('Schmid'))
def test_Schmid_invalid(self):
with pytest.raises(KeyError):
Orientation(lattice='fcc').Schmid()
### vectorization tests ### ### vectorization tests ###
@pytest.mark.parametrize('lattice',['hP','cI','cF']) # tI not included yet @pytest.mark.parametrize('lattice',['hP','cI','cF']) # tI not included yet
@ -505,3 +510,7 @@ class TestOrientation:
for loc in np.random.randint(0,blend,(10,len(blend))): for loc in np.random.randint(0,blend,(10,len(blend))):
assert np.allclose(o[tuple(loc[:len(o.shape)])].to_pole(uvw=v[tuple(loc[-len(v.shape[:-1]):])]), assert np.allclose(o[tuple(loc[:len(o.shape)])].to_pole(uvw=v[tuple(loc[-len(v.shape[:-1]):])]),
o.to_pole(uvw=v)[tuple(loc)]) o.to_pole(uvw=v)[tuple(loc)])
def test_mul_invalid(self):
with pytest.raises(TypeError):
Orientation.from_random(lattice='cF')*np.ones(3)

3
python/tests/test_Result.py
View File

@ -102,6 +102,9 @@ class TestResult:
with pytest.raises(AttributeError): with pytest.raises(AttributeError):
default.view('invalid',True) default.view('invalid',True)
def test_add_invalid(self,default):
default.add_absolute('xxxx')
def test_add_absolute(self,default): def test_add_absolute(self,default):
default.add_absolute('F_e') default.add_absolute('F_e')
in_memory = np.abs(default.place('F_e')) in_memory = np.abs(default.place('F_e'))

9
python/tests/test_Rotation.py
View File

@ -792,6 +792,11 @@ class TestRotation:
R = Rotation.from_random(shape,rng_seed=1) R = Rotation.from_random(shape,rng_seed=1)
assert R == R if shape is None else (R == R).all() assert R == R if shape is None else (R == R).all()
@pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_allclose(self,shape):
R = Rotation.from_random(shape,rng_seed=1)
assert R.allclose(R)
@pytest.mark.parametrize('shape',[None,5,(4,6)]) @pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_unequal(self,shape): def test_unequal(self,shape):
R = Rotation.from_random(shape,rng_seed=1) R = Rotation.from_random(shape,rng_seed=1)
@ -1124,3 +1129,7 @@ class TestRotation:
weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights) weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5 assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5
def test_mul_invalid(self):
with pytest.raises(TypeError):
Rotation.from_random()*np.ones(3)

4
python/tests/test_VTK.py
View File

@ -28,6 +28,10 @@ class TestVTK:
def _patch_execution_stamp(self, patch_execution_stamp): def _patch_execution_stamp(self, patch_execution_stamp):
print('patched damask.util.execution_stamp') print('patched damask.util.execution_stamp')
def test_show(sef,default,monkeypatch):
monkeypatch.delenv('DISPLAY',raising=False)
default.show()
def test_rectilinearGrid(self,tmp_path): def test_rectilinearGrid(self,tmp_path):
cells = np.random.randint(5,10,3)*2 cells = np.random.randint(5,10,3)*2
size = np.random.random(3) + 1.0 size = np.random.random(3) + 1.0

250
src/YAML_types.f90
View File

@ -62,6 +62,8 @@ module YAML_types
tNode_get_byKey_as1dString => tNode_get_byKey_as1dString tNode_get_byKey_as1dString => tNode_get_byKey_as1dString
procedure :: & procedure :: &
getKey => tNode_get_byIndex_asKey getKey => tNode_get_byIndex_asKey
procedure :: &
Keys => tNode_getKeys
procedure :: & procedure :: &
getIndex => tNode_get_byKey_asIndex getIndex => tNode_get_byKey_asIndex
procedure :: & procedure :: &
@ -117,7 +119,7 @@ module YAML_types
type, extends(tNode), public :: tList type, extends(tNode), public :: tList
class(tItem), pointer :: first => null() class(tItem), pointer :: first => NULL()
contains contains
procedure :: asFormattedString => tList_asFormattedString procedure :: asFormattedString => tList_asFormattedString
@ -144,8 +146,8 @@ module YAML_types
type :: tItem type :: tItem
character(len=:), allocatable :: key character(len=:), allocatable :: key
class(tNode), pointer :: node => null() class(tNode), pointer :: node => NULL()
class(tItem), pointer :: next => null() class(tItem), pointer :: next => NULL()
contains contains
final :: tItem_finalize final :: tItem_finalize
@ -221,22 +223,22 @@ subroutine selfTest
select type(s1) select type(s1)
class is(tScalar) class is(tScalar)
s1 = '2' s1 = '2'
endselect end select
select type(s2) select type(s2)
class is(tScalar) class is(tScalar)
s2 = '3' s2 = '3'
endselect end select
select type(s3) select type(s3)
class is(tScalar) class is(tScalar)
s3 = '4' s3 = '4'
endselect end select
select type(s4) select type(s4)
class is(tScalar) class is(tScalar)
s4 = '5' s4 = '5'
endselect end select
allocate(tList::l1) allocate(tList::l1)
@ -249,14 +251,14 @@ subroutine selfTest
if (any(dNeq(l1%as1dFloat(),[2.0_pReal,3.0_pReal]))) error stop 'tList_as1dFloat' if (any(dNeq(l1%as1dFloat(),[2.0_pReal,3.0_pReal]))) error stop 'tList_as1dFloat'
if (n%get_asInt(1) /= 2) error stop 'byIndex_asInt' if (n%get_asInt(1) /= 2) error stop 'byIndex_asInt'
if (dNeq(n%get_asFloat(2),3.0_pReal)) error stop 'byIndex_asFloat' if (dNeq(n%get_asFloat(2),3.0_pReal)) error stop 'byIndex_asFloat'
endselect end select
allocate(tList::l3) allocate(tList::l3)
select type(l3) select type(l3)
class is(tList) class is(tList)
call l3%append(s3) call l3%append(s3)
call l3%append(s4) call l3%append(s4)
endselect end select
allocate(tList::l2) allocate(tList::l2)
select type(l2) select type(l2)
@ -332,9 +334,12 @@ function tNode_asScalar(self) result(scalar)
class(tNode), intent(in), target :: self class(tNode), intent(in), target :: self
class(tScalar), pointer :: scalar class(tScalar), pointer :: scalar
select type(self) select type(self)
class is(tScalar) class is(tScalar)
scalar => self scalar => self
class default
nullify(scalar)
end select end select
end function tNode_asScalar end function tNode_asScalar
@ -348,9 +353,12 @@ function tNode_asList(self) result(list)
class(tNode), intent(in), target :: self class(tNode), intent(in), target :: self
class(tList), pointer :: list class(tList), pointer :: list
select type(self) select type(self)
class is(tList) class is(tList)
list => self list => self
class default
nullify(list)
end select end select
end function tNode_asList end function tNode_asList
@ -364,9 +372,12 @@ function tNode_asDict(self) result(dict)
class(tNode), intent(in), target :: self class(tNode), intent(in), target :: self
class(tDict), pointer :: dict class(tDict), pointer :: dict
select type(self) select type(self)
class is(tDict) class is(tDict)
dict => self dict => self
class default
nullify(dict)
end select end select
end function tNode_asDict end function tNode_asDict
@ -385,12 +396,13 @@ function tNode_get_byIndex(self,i) result(node)
class(tItem), pointer :: item class(tItem), pointer :: item
integer :: j integer :: j
select type(self) select type(self)
class is(tList) class is(tList)
self_ => self%asList() self_ => self%asList()
class default class default
call IO_error(706,ext_msg='Expected list') call IO_error(706,ext_msg='Expected list')
endselect end select
item => self_%first item => self_%first
@ -409,15 +421,14 @@ end function tNode_get_byIndex
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_asFloat(self,i) result(nodeAsFloat) function tNode_get_byIndex_asFloat(self,i) result(nodeAsFloat)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
integer, intent(in) :: i integer, intent(in) :: i
real(pReal) :: nodeAsFloat real(pReal) :: nodeAsFloat
class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
node => self%get(i)
select type(node) select type(node => self%get(i))
class is(tScalar) class is(tScalar)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsFloat = scalar%asFloat() nodeAsFloat = scalar%asFloat()
@ -433,15 +444,15 @@ end function tNode_get_byIndex_asFloat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_asInt(self,i) result(nodeAsInt) function tNode_get_byIndex_asInt(self,i) result(nodeAsInt)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
integer, intent(in) :: i integer, intent(in) :: i
integer :: nodeAsInt integer :: nodeAsInt
class(tNode), pointer :: node class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
node => self%get(i)
select type(node) select type(node => self%get(i))
class is(tScalar) class is(tScalar)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsInt = scalar%asInt() nodeAsInt = scalar%asInt()
@ -457,21 +468,20 @@ end function tNode_get_byIndex_asInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_asBool(self,i) result(nodeAsBool) function tNode_get_byIndex_asBool(self,i) result(nodeAsBool)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
integer, intent(in) :: i integer, intent(in) :: i
logical :: nodeAsBool logical :: nodeAsBool
class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
node => self%get(i)
select type(node) select type(node => self%get(i))
class is(tScalar) class is(tScalar)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsBool = scalar%asBool() nodeAsBool = scalar%asBool()
class default class default
call IO_error(706,ext_msg='Expected scalar Boolean') call IO_error(706,ext_msg='Expected scalar Boolean')
endselect end select
end function tNode_get_byIndex_asBool end function tNode_get_byIndex_asBool
@ -481,21 +491,20 @@ end function tNode_get_byIndex_asBool
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_asString(self,i) result(nodeAsString) function tNode_get_byIndex_asString(self,i) result(nodeAsString)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
integer, intent(in) :: i integer, intent(in) :: i
character(len=:), allocatable :: nodeAsString character(len=:), allocatable :: nodeAsString
class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
node => self%get(i)
select type(node) select type(node => self%get(i))
class is(tScalar) class is(tScalar)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsString = scalar%asString() nodeAsString = scalar%asString()
class default class default
call IO_error(706,ext_msg='Expected scalar string') call IO_error(706,ext_msg='Expected scalar string')
endselect end select
end function tNode_get_byIndex_asString end function tNode_get_byIndex_asString
@ -505,21 +514,20 @@ end function tNode_get_byIndex_asString
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_as1dFloat(self,i) result(nodeAs1dFloat) function tNode_get_byIndex_as1dFloat(self,i) result(nodeAs1dFloat)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
integer, intent(in) :: i integer, intent(in) :: i
real(pReal), dimension(:), allocatable :: nodeAs1dFloat real(pReal), dimension(:), allocatable :: nodeAs1dFloat
class(tNode), pointer :: node
class(tList), pointer :: list class(tList), pointer :: list
node => self%get(i)
select type(node) select type(node => self%get(i))
class is(tList) class is(tList)
list => node%asList() list => node%asList()
nodeAs1dFloat = list%as1dFloat() nodeAs1dFloat = list%as1dFloat()
class default class default
call IO_error(706,ext_msg='Expected list of floats') call IO_error(706,ext_msg='Expected list of floats')
endselect end select
end function tNode_get_byIndex_as1dFloat end function tNode_get_byIndex_as1dFloat
@ -529,21 +537,20 @@ end function tNode_get_byIndex_as1dFloat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_as1dInt(self,i) result(nodeAs1dInt) function tNode_get_byIndex_as1dInt(self,i) result(nodeAs1dInt)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
integer, intent(in) :: i integer, intent(in) :: i
integer, dimension(:), allocatable :: nodeAs1dInt integer, dimension(:), allocatable :: nodeAs1dInt
class(tNode), pointer :: node
class(tList), pointer :: list class(tList), pointer :: list
node => self%get(i)
select type(node) select type(node => self%get(i))
class is(tList) class is(tList)
list => node%asList() list => node%asList()
nodeAs1dInt = list%as1dInt() nodeAs1dInt = list%as1dInt()
class default class default
call IO_error(706,ext_msg='Expected list of integers') call IO_error(706,ext_msg='Expected list of integers')
endselect end select
end function tNode_get_byIndex_as1dInt end function tNode_get_byIndex_as1dInt
@ -553,21 +560,20 @@ end function tNode_get_byIndex_as1dInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_as1dBool(self,i) result(nodeAs1dBool) function tNode_get_byIndex_as1dBool(self,i) result(nodeAs1dBool)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
integer, intent(in) :: i integer, intent(in) :: i
logical, dimension(:), allocatable :: nodeAs1dBool logical, dimension(:), allocatable :: nodeAs1dBool
class(tNode), pointer :: node
class(tList), pointer :: list class(tList), pointer :: list
node => self%get(i)
select type(node) select type(node => self%get(i))
class is(tList) class is(tList)
list => node%asList() list => node%asList()
nodeAs1dBool = list%as1dBool() nodeAs1dBool = list%as1dBool()
class default class default
call IO_error(706,ext_msg='Expected list of Booleans') call IO_error(706,ext_msg='Expected list of Booleans')
endselect end select
end function tNode_get_byIndex_as1dBool end function tNode_get_byIndex_as1dBool
@ -577,21 +583,20 @@ end function tNode_get_byIndex_as1dBool
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byIndex_as1dString(self,i) result(nodeAs1dString) function tNode_get_byIndex_as1dString(self,i) result(nodeAs1dString)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
integer, intent(in) :: i integer, intent(in) :: i
character(len=:), allocatable, dimension(:) :: nodeAs1dString character(len=:), allocatable, dimension(:) :: nodeAs1dString
class(tNode), pointer :: node
type(tList), pointer :: list type(tList), pointer :: list
node => self%get(i)
select type(node) select type(node => self%get(i))
class is(tList) class is(tList)
list => node%asList() list => node%asList()
nodeAs1dString = list%as1dString() nodeAs1dString = list%as1dString()
class default class default
call IO_error(706,ext_msg='Expected list of strings') call IO_error(706,ext_msg='Expected list of strings')
endselect end select
end function tNode_get_byIndex_as1dString end function tNode_get_byIndex_as1dString
@ -609,22 +614,50 @@ function tNode_get_byIndex_asKey(self,i) result(key)
type(tDict), pointer :: dict type(tDict), pointer :: dict
type(tItem), pointer :: item type(tItem), pointer :: item
select type(self) select type(self)
class is(tDict) class is(tDict)
dict => self%asDict() dict => self%asDict()
item => dict%first item => dict%first
do j = 1, min(i,dict%length)-1 do j = 1, min(i,dict%length)-1
item => item%next item => item%next
enddo end do
class default class default
call IO_error(706,ext_msg='Expected dict') call IO_error(706,ext_msg='Expected dict')
endselect end select
key = item%key key = item%key
end function tNode_get_byIndex_asKey end function tNode_get_byIndex_asKey
!--------------------------------------------------------------------------------------------------
!> @brief Get all keys from a dictionary
!--------------------------------------------------------------------------------------------------
function tNode_getKeys(self) result(keys)
class(tNode), intent(in) :: self
character(len=:), dimension(:), allocatable :: keys
character(len=pStringLen), dimension(:), allocatable :: temp
integer :: j, l
allocate(temp(self%length))
l = 0
do j = 1, self%length
temp(j) = self%getKey(j)
l = max(len_trim(temp(j)),l)
end do
allocate(character(l)::keys(self%length))
do j = 1, self%length
keys(j) = trim(temp(j))
end do
end function tNode_getKeys
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
!> @brief Checks if a given key/item is present in the dict/list !> @brief Checks if a given key/item is present in the dict/list
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
@ -658,7 +691,7 @@ function tNode_contains(self,k) result(exists)
enddo enddo
class default class default
call IO_error(706,ext_msg='Expected list or dict') call IO_error(706,ext_msg='Expected list or dict')
endselect end select
end function tNode_contains end function tNode_contains
@ -686,7 +719,7 @@ function tNode_get_byKey(self,k,defaultVal) result(node)
self_ => self%asDict() self_ => self%asDict()
class default class default
call IO_error(706,ext_msg='Expected dict for key '//k) call IO_error(706,ext_msg='Expected dict for key '//k)
endselect end select
j = 1 j = 1
item => self_%first item => self_%first
@ -713,23 +746,22 @@ end function tNode_get_byKey
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_asFloat(self,k,defaultVal) result(nodeAsFloat) function tNode_get_byKey_asFloat(self,k,defaultVal) result(nodeAsFloat)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
real(pReal), intent(in),optional :: defaultVal real(pReal), intent(in), optional :: defaultVal
real(pReal) :: nodeAsFloat real(pReal) :: nodeAsFloat
class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
if (self%contains(k)) then if (self%contains(k)) then
node => self%get(k) select type(node => self%get(k))
select type(node)
class is(tScalar) class is(tScalar)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsFloat = scalar%asFloat() nodeAsFloat = scalar%asFloat()
class default class default
call IO_error(706,ext_msg='Expected scalar float for key '//k) call IO_error(706,ext_msg='Expected scalar float for key '//k)
endselect end select
elseif (present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsFloat = defaultVal nodeAsFloat = defaultVal
else else
@ -744,23 +776,22 @@ end function tNode_get_byKey_asFloat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_asInt(self,k,defaultVal) result(nodeAsInt) function tNode_get_byKey_asInt(self,k,defaultVal) result(nodeAsInt)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
integer, intent(in),optional :: defaultVal integer, intent(in), optional :: defaultVal
integer :: nodeAsInt integer :: nodeAsInt
class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
if (self%contains(k)) then if (self%contains(k)) then
node => self%get(k) select type(node => self%get(k))
select type(node)
class is(tScalar) class is(tScalar)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsInt = scalar%asInt() nodeAsInt = scalar%asInt()
class default class default
call IO_error(706,ext_msg='Expected scalar integer for key '//k) call IO_error(706,ext_msg='Expected scalar integer for key '//k)
endselect end select
elseif (present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsInt = defaultVal nodeAsInt = defaultVal
else else
@ -775,23 +806,22 @@ end function tNode_get_byKey_asInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_asBool(self,k,defaultVal) result(nodeAsBool) function tNode_get_byKey_asBool(self,k,defaultVal) result(nodeAsBool)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
logical, intent(in),optional :: defaultVal logical, intent(in), optional :: defaultVal
logical :: nodeAsBool logical :: nodeAsBool
class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
if (self%contains(k)) then if (self%contains(k)) then
node => self%get(k) select type(node => self%get(k))
select type(node)
class is(tScalar) class is(tScalar)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsBool = scalar%asBool() nodeAsBool = scalar%asBool()
class default class default
call IO_error(706,ext_msg='Expected scalar Boolean for key '//k) call IO_error(706,ext_msg='Expected scalar Boolean for key '//k)
endselect end select
elseif (present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsBool = defaultVal nodeAsBool = defaultVal
else else
@ -806,23 +836,22 @@ end function tNode_get_byKey_asBool
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_asString(self,k,defaultVal) result(nodeAsString) function tNode_get_byKey_asString(self,k,defaultVal) result(nodeAsString)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
character(len=*), intent(in),optional :: defaultVal character(len=*), intent(in), optional :: defaultVal
character(len=:), allocatable :: nodeAsString character(len=:), allocatable :: nodeAsString
class(tNode), pointer :: node
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
if (self%contains(k)) then if (self%contains(k)) then
node => self%get(k) select type(node => self%get(k))
select type(node)
class is(tScalar) class is(tScalar)
scalar => node%asScalar() scalar => node%asScalar()
nodeAsString = scalar%asString() nodeAsString = scalar%asString()
class default class default
call IO_error(706,ext_msg='Expected scalar string for key '//k) call IO_error(706,ext_msg='Expected scalar string for key '//k)
endselect end select
elseif (present(defaultVal)) then elseif (present(defaultVal)) then
nodeAsString = defaultVal nodeAsString = defaultVal
else else
@ -837,25 +866,24 @@ end function tNode_get_byKey_asString
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as1dFloat(self,k,defaultVal,requiredSize) result(nodeAs1dFloat) function tNode_get_byKey_as1dFloat(self,k,defaultVal,requiredSize) result(nodeAs1dFloat)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
real(pReal), intent(in), dimension(:), optional :: defaultVal real(pReal), intent(in), dimension(:), optional :: defaultVal
integer, intent(in), optional :: requiredSize integer, intent(in), optional :: requiredSize
real(pReal), dimension(:), allocatable :: nodeAs1dFloat real(pReal), dimension(:), allocatable :: nodeAs1dFloat
class(tNode), pointer :: node
type(tList), pointer :: list type(tList), pointer :: list
if (self%contains(k)) then if (self%contains(k)) then
node => self%get(k) select type(node => self%get(k))
select type(node)
class is(tList) class is(tList)
list => node%asList() list => node%asList()
nodeAs1dFloat = list%as1dFloat() nodeAs1dFloat = list%as1dFloat()
class default class default
call IO_error(706,ext_msg='Expected 1D float array for key '//k) call IO_error(706,ext_msg='Expected 1D float array for key '//k)
endselect end select
elseif (present(defaultVal)) then elseif (present(defaultVal)) then
nodeAs1dFloat = defaultVal nodeAs1dFloat = defaultVal
else else
@ -874,25 +902,24 @@ end function tNode_get_byKey_as1dFloat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as2dFloat(self,k,defaultVal,requiredShape) result(nodeAs2dFloat) function tNode_get_byKey_as2dFloat(self,k,defaultVal,requiredShape) result(nodeAs2dFloat)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
real(pReal), intent(in), dimension(:,:), optional :: defaultVal real(pReal), intent(in), dimension(:,:), optional :: defaultVal
integer, intent(in), dimension(2), optional :: requiredShape integer, intent(in), dimension(2), optional :: requiredShape
real(pReal), dimension(:,:), allocatable :: nodeAs2dFloat real(pReal), dimension(:,:), allocatable :: nodeAs2dFloat
class(tNode), pointer :: node
type(tList), pointer :: rows type(tList), pointer :: rows
if(self%contains(k)) then if(self%contains(k)) then
node => self%get(k) select type(node => self%get(k))
select type(node)
class is(tList) class is(tList)
rows => node%asList() rows => node%asList()
nodeAs2dFloat = rows%as2dFloat() nodeAs2dFloat = rows%as2dFloat()
class default class default
call IO_error(706,ext_msg='Expected 2D float array for key '//k) call IO_error(706,ext_msg='Expected 2D float array for key '//k)
endselect end select
elseif(present(defaultVal)) then elseif(present(defaultVal)) then
nodeAs2dFloat = defaultVal nodeAs2dFloat = defaultVal
else else
@ -911,24 +938,22 @@ end function tNode_get_byKey_as2dFloat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as1dInt(self,k,defaultVal,requiredSize) result(nodeAs1dInt) function tNode_get_byKey_as1dInt(self,k,defaultVal,requiredSize) result(nodeAs1dInt)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
integer, dimension(:), intent(in), optional :: defaultVal integer, dimension(:), intent(in), optional :: defaultVal
integer, intent(in), optional :: requiredSize integer, intent(in), optional :: requiredSize
integer, dimension(:), allocatable :: nodeAs1dInt integer, dimension(:), allocatable :: nodeAs1dInt
class(tNode), pointer :: node
type(tList), pointer :: list type(tList), pointer :: list
if (self%contains(k)) then if (self%contains(k)) then
node => self%get(k) select type(node => self%get(k))
select type(node)
class is(tList) class is(tList)
list => node%asList() list => node%asList()
nodeAs1dInt = list%as1dInt() nodeAs1dInt = list%as1dInt()
class default class default
call IO_error(706,ext_msg='Expected 1D integer array for key '//k) call IO_error(706,ext_msg='Expected 1D integer array for key '//k)
endselect end select
elseif (present(defaultVal)) then elseif (present(defaultVal)) then
nodeAs1dInt = defaultVal nodeAs1dInt = defaultVal
else else
@ -947,23 +972,22 @@ end function tNode_get_byKey_as1dInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as1dBool(self,k,defaultVal) result(nodeAs1dBool) function tNode_get_byKey_as1dBool(self,k,defaultVal) result(nodeAs1dBool)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
logical, dimension(:), intent(in), optional :: defaultVal logical, dimension(:), intent(in), optional :: defaultVal
logical, dimension(:), allocatable :: nodeAs1dBool logical, dimension(:), allocatable :: nodeAs1dBool
class(tNode), pointer :: node
type(tList), pointer :: list type(tList), pointer :: list
if (self%contains(k)) then if (self%contains(k)) then
node => self%get(k) select type(node => self%get(k))
select type(node)
class is(tList) class is(tList)
list => node%asList() list => node%asList()
nodeAs1dBool = list%as1dBool() nodeAs1dBool = list%as1dBool()
class default class default
call IO_error(706,ext_msg='Expected 1D Boolean array for key '//k) call IO_error(706,ext_msg='Expected 1D Boolean array for key '//k)
endselect end select
elseif (present(defaultVal)) then elseif (present(defaultVal)) then
nodeAs1dBool = defaultVal nodeAs1dBool = defaultVal
else else
@ -978,23 +1002,22 @@ end function tNode_get_byKey_as1dBool
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as1dString(self,k,defaultVal) result(nodeAs1dString) function tNode_get_byKey_as1dString(self,k,defaultVal) result(nodeAs1dString)
class(tNode), intent(in), target :: self class(tNode), intent(in) :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
character(len=*), intent(in), dimension(:), optional :: defaultVal character(len=*), intent(in), dimension(:), optional :: defaultVal
character(len=:), allocatable, dimension(:) :: nodeAs1dString character(len=:), allocatable, dimension(:) :: nodeAs1dString
class(tNode), pointer :: node
type(tList), pointer :: list type(tList), pointer :: list
if (self%contains(k)) then if (self%contains(k)) then
node => self%get(k) select type(node => self%get(k))
select type(node)
class is(tList) class is(tList)
list => node%asList() list => node%asList()
nodeAs1dString = list%as1dString() nodeAs1dString = list%as1dString()
class default class default
call IO_error(706,ext_msg='Expected 1D string array for key '//k) call IO_error(706,ext_msg='Expected 1D string array for key '//k)
endselect end select
elseif (present(defaultVal)) then elseif (present(defaultVal)) then
nodeAs1dString = defaultVal nodeAs1dString = defaultVal
else else
@ -1019,7 +1042,7 @@ function output_as1dString(self) result(output) !ToDo: SR: Re
allocate(output(output_list%length)) allocate(output(output_list%length))
do o = 1, output_list%length do o = 1, output_list%length
output(o) = output_list%get_asString(o) output(o) = output_list%get_asString(o)
enddo end do
end function output_as1dString end function output_as1dString
@ -1042,7 +1065,7 @@ function tNode_get_byKey_asIndex(self,key) result(keyIndex)
do while (associated(item%next) .and. item%key /= key) do while (associated(item%next) .and. item%key /= key)
item => item%next item => item%next
keyIndex = keyIndex+1 keyIndex = keyIndex+1
enddo end do
if (item%key /= key) call IO_error(140,ext_msg=key) if (item%key /= key) call IO_error(140,ext_msg=key)
@ -1087,7 +1110,7 @@ recursive function tList_asFormattedString(self,indent) result(str)
if (i /= 1) str = str//repeat(' ',indent_) if (i /= 1) str = str//repeat(' ',indent_)
str = str//'- '//item%node%asFormattedString(indent_+2) str = str//'- '//item%node%asFormattedString(indent_+2)
item => item%next item => item%next
enddo end do
end function tList_asFormattedString end function tList_asFormattedString
@ -1119,9 +1142,9 @@ recursive function tDict_asFormattedString(self,indent) result(str)
str = str//trim(item%key)//': '//item%node%asFormattedString(indent_+len_trim(item%key)+2) str = str//trim(item%key)//': '//item%node%asFormattedString(indent_+len_trim(item%key)+2)
class default class default
str = str//trim(item%key)//':'//IO_EOL//repeat(' ',indent_+2)//item%node%asFormattedString(indent_+2) str = str//trim(item%key)//':'//IO_EOL//repeat(' ',indent_+2)//item%node%asFormattedString(indent_+2)
endselect end select
item => item%next item => item%next
enddo end do
end function tDict_asFormattedString end function tDict_asFormattedString
@ -1190,13 +1213,14 @@ function tList_as1dFloat(self)
type(tItem), pointer :: item type(tItem), pointer :: item
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
allocate(tList_as1dFloat(self%length)) allocate(tList_as1dFloat(self%length))
item => self%first item => self%first
do i = 1, self%length do i = 1, self%length
scalar => item%node%asScalar() scalar => item%node%asScalar()
tList_as1dFloat(i) = scalar%asFloat() tList_as1dFloat(i) = scalar%asFloat()
item => item%next item => item%next
enddo end do
end function tList_as1dFloat end function tList_as1dFloat
@ -1213,6 +1237,7 @@ function tList_as2dFloat(self)
class(tNode), pointer :: row class(tNode), pointer :: row
type(tList), pointer :: row_data type(tList), pointer :: row_data
row => self%get(1) row => self%get(1)
row_data => row%asList() row_data => row%asList()
allocate(tList_as2dFloat(self%length,row_data%length)) allocate(tList_as2dFloat(self%length,row_data%length))
@ -1220,9 +1245,9 @@ function tList_as2dFloat(self)
do i=1,self%length do i=1,self%length
row => self%get(i) row => self%get(i)
row_data => row%asList() row_data => row%asList()
if(row_data%length /= size(tList_as2dFloat,2)) call IO_error(709,ext_msg='Varying number of columns') if (row_data%length /= size(tList_as2dFloat,2)) call IO_error(709,ext_msg='Varying number of columns')
tList_as2dFloat(i,:) = self%get_as1dFloat(i) tList_as2dFloat(i,:) = self%get_as1dFloat(i)
enddo end do
end function tList_as2dFloat end function tList_as2dFloat
@ -1239,13 +1264,14 @@ function tList_as1dInt(self)
type(tItem), pointer :: item type(tItem), pointer :: item
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
allocate(tList_as1dInt(self%length)) allocate(tList_as1dInt(self%length))
item => self%first item => self%first
do i = 1, self%length do i = 1, self%length
scalar => item%node%asScalar() scalar => item%node%asScalar()
tList_as1dInt(i) = scalar%asInt() tList_as1dInt(i) = scalar%asInt()
item => item%next item => item%next
enddo end do
end function tList_as1dInt end function tList_as1dInt
@ -1262,13 +1288,14 @@ function tList_as1dBool(self)
type(tItem), pointer :: item type(tItem), pointer :: item
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
allocate(tList_as1dBool(self%length)) allocate(tList_as1dBool(self%length))
item => self%first item => self%first
do i = 1, self%length do i = 1, self%length
scalar => item%node%asScalar() scalar => item%node%asScalar()
tList_as1dBool(i) = scalar%asBool() tList_as1dBool(i) = scalar%asBool()
item => item%next item => item%next
enddo end do
end function tList_as1dBool end function tList_as1dBool
@ -1285,13 +1312,14 @@ function tList_as1dString(self)
type(tItem), pointer :: item type(tItem), pointer :: item
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
len_max = 0 len_max = 0
item => self%first item => self%first
do i = 1, self%length do i = 1, self%length
scalar => item%node%asScalar() scalar => item%node%asScalar()
len_max = max(len_max, len_trim(scalar%asString())) len_max = max(len_max, len_trim(scalar%asString()))
item => item%next item => item%next
enddo end do
allocate(character(len=len_max) :: tList_as1dString(self%length)) allocate(character(len=len_max) :: tList_as1dString(self%length))
item => self%first item => self%first

101
src/grid/grid_damage_spectral.f90
View File

@ -39,9 +39,8 @@ module grid_damage_spectral
type(tSolutionParams) :: params type(tSolutionParams) :: params
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
SNES :: damage_snes SNES :: SNES_damage
Vec :: solution_vec Vec :: solution_vec
PetscInt :: xstart, xend, ystart, yend, zstart, zend
real(pReal), dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
phi_current, & !< field of current damage phi_current, & !< field of current damage
phi_lastInc, & !< field of previous damage phi_lastInc, & !< field of previous damage
@ -105,10 +104,18 @@ subroutine grid_damage_spectral_init()
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! init fields
allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,damage_snes,err_PETSc); CHKERRQ(err_PETSc) call SNESCreate(PETSC_COMM_WORLD,SNES_damage,err_PETSc)
call SNESSetOptionsPrefix(damage_snes,'damage_',err_PETSc);CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(SNES_damage,'damage_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt) localK(worldrank) = int(grid3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@ -122,39 +129,41 @@ subroutine grid_damage_spectral_init()
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
damage_grid,err_PETSc) ! handle, error damage_grid,err_PETSc) ! handle, error
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetDM(damage_snes,damage_grid,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da call DMsetFromOptions(damage_grid,err_PETSc)
call DMsetFromOptions(damage_grid,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMsetUp(damage_grid,err_PETSc); CHKERRQ(err_PETSc) call DMsetUp(damage_grid,err_PETSc)
call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor) CHKERRQ(err_PETSc)
call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetFromOptions(damage_snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments call SNESSetDM(SNES_damage,damage_grid,err_PETSc)
call SNESGetType(damage_snes,snes_type,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetFromOptions(SNES_damage,err_PETSc) ! pull it all together with additional CLI arguments
CHKERRQ(err_PETSc)
call SNESGetType(SNES_damage,snes_type,err_PETSc)
CHKERRQ(err_PETSc)
if (trim(snes_type) == 'vinewtonrsls' .or. & if (trim(snes_type) == 'vinewtonrsls' .or. &
trim(snes_type) == 'vinewtonssls') then trim(snes_type) == 'vinewtonssls') then
call DMGetGlobalVector(damage_grid,lBound,err_PETSc); CHKERRQ(err_PETSc) call DMGetGlobalVector(damage_grid,lBound,err_PETSc)
call DMGetGlobalVector(damage_grid,uBound,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecSet(lBound,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc) call DMGetGlobalVector(damage_grid,uBound,err_PETSc)
call VecSet(uBound,1.0_pReal,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESVISetVariableBounds(damage_snes,lBound,uBound,err_PETSc) ! variable bounds for variational inequalities like contact mechanics, damage etc. call VecSet(lBound,0.0_pReal,err_PETSc)
call DMRestoreGlobalVector(damage_grid,lBound,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMRestoreGlobalVector(damage_grid,uBound,err_PETSc); CHKERRQ(err_PETSc) call VecSet(uBound,1.0_pReal,err_PETSc)
end if CHKERRQ(err_PETSc)
call SNESVISetVariableBounds(SNES_damage,lBound,uBound,err_PETSc) ! variable bounds for variational inequalities
!-------------------------------------------------------------------------------------------------- CHKERRQ(err_PETSc)
! init fields call DMRestoreGlobalVector(damage_grid,lBound,err_PETSc)
call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,err_PETSc) CHKERRQ(err_PETSc)
call DMRestoreGlobalVector(damage_grid,uBound,err_PETSc)
CHKERRQ(err_PETSc)
end if
call VecSet(solution_vec,1.0_pReal,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
xend = xstart + xend - 1
yend = ystart + yend - 1
zend = zstart + zend - 1
allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
call VecSet(solution_vec,1.0_pReal,err_PETSc); CHKERRQ(err_PETSc) call updateReference()
call updateReference
end subroutine grid_damage_spectral_init end subroutine grid_damage_spectral_init
@ -181,9 +190,9 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
! set module wide availabe data ! set module wide availabe data
params%Delta_t = Delta_t params%Delta_t = Delta_t
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,err_PETSc) call SNESSolve(SNES_damage,PETSC_NULL_VEC,solution_vec,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESGetConvergedReason(damage_snes,reason,err_PETSc) call SNESGetConvergedReason(SNES_damage,reason,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
if (reason < 1) then if (reason < 1) then
@ -209,8 +218,10 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
call homogenization_set_phi(phi_current(i,j,k),ce) call homogenization_set_phi(phi_current(i,j,k),ce)
end do; end do; end do end do; end do; end do
call VecMin(solution_vec,devNull,phi_min,err_PETSc); CHKERRQ(err_PETSc) call VecMin(solution_vec,devNull,phi_min,err_PETSc)
call VecMax(solution_vec,devNull,phi_max,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecMax(solution_vec,devNull,phi_max,err_PETSc)
CHKERRQ(err_PETSc)
if (solution%converged) & if (solution%converged) &
print'(/,1x,a)', '... nonlocal damage converged .....................................' print'(/,1x,a)', '... nonlocal damage converged .....................................'
print'(/,1x,a,f8.6,2x,f8.6,2x,e11.4)', 'Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm print'(/,1x,a,f8.6,2x,f8.6,2x,e11.4)', 'Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm
@ -228,7 +239,7 @@ subroutine grid_damage_spectral_forward(cutBack)
logical, intent(in) :: cutBack logical, intent(in) :: cutBack
integer :: i, j, k, ce integer :: i, j, k, ce
DM :: dm_local DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: phi_PETSc
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
if (cutBack) then if (cutBack) then
@ -236,11 +247,12 @@ subroutine grid_damage_spectral_forward(cutBack)
phi_stagInc = phi_lastInc phi_stagInc = phi_lastInc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reverting damage field state ! reverting damage field state
call SNESGetDM(damage_snes,dm_local,err_PETSc); CHKERRQ(err_PETSc) call SNESGetDM(SNES_damage,dm_local,err_PETSc)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,err_PETSc) !< get the data out of PETSc to work with
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current call DMDAVecGetArrayF90(dm_local,solution_vec,phi_PETSc,err_PETSc) !< get the data out of PETSc to work with
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,err_PETSc) CHKERRQ(err_PETSc)
phi_PETSc = phi_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,phi_PETSc,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
@ -258,7 +270,7 @@ end subroutine grid_damage_spectral_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forms the spectral damage residual vector !> @brief forms the spectral damage residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err) subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in in
@ -267,9 +279,9 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
x_scal x_scal
PetscScalar, dimension( & PetscScalar, dimension( &
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
f_scal r
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: dummy_err PetscErrorCode :: err_PETSc
integer :: i, j, k, ce integer :: i, j, k, ce
@ -310,7 +322,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current r = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current
err_PETSc = 0
end subroutine formResidual end subroutine formResidual

182
src/grid/grid_mech_FEM.f90
View File

@ -50,7 +50,7 @@ module grid_mechanical_FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM :: mechanical_grid DM :: mechanical_grid
SNES :: mechanical_snes SNES :: SNES_mechanical
Vec :: solution_current, solution_lastInc, solution_rate Vec :: solution_current, solution_lastInc, solution_rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -60,7 +60,6 @@ module grid_mechanical_FEM
real(pReal), dimension(3) :: delta real(pReal), dimension(3) :: delta
real(pReal), dimension(3,8) :: BMat real(pReal), dimension(3,8) :: BMat
real(pReal), dimension(8,8) :: HGMat real(pReal), dimension(8,8) :: HGMat
PetscInt :: xstart,ystart,zstart,xend,yend,zend
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
@ -146,7 +145,7 @@ subroutine grid_mechanical_FEM_init
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS, & call PetscOptionsInsertString(PETSC_NULL_OPTIONS, &
'-mechanical_snes_type newtonls -mechanical_ksp_type fgmres & '-mechanical_snes_type newtonls -mechanical_ksp_type fgmres &
&-mechanical_ksp_max_it 25 -mechanical_mg_levels_ksp_type chebyshev', & &-mechanical_ksp_max_it 25', &
err_PETSc) err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc)
@ -160,9 +159,9 @@ subroutine grid_mechanical_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,err_PETSc) call SNESCreate(PETSC_COMM_WORLD,SNES_mechanical,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',err_PETSc) call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
localK = 0_pPetscInt localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt) localK(worldrank) = int(grid3,pPetscInt)
@ -177,8 +176,6 @@ subroutine grid_mechanical_FEM_init
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
mechanical_grid,err_PETSc) mechanical_grid,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetDM(mechanical_snes,mechanical_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetFromOptions(mechanical_grid,err_PETSc) call DMsetFromOptions(mechanical_grid,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMsetUp(mechanical_grid,err_PETSc) call DMsetUp(mechanical_grid,err_PETSc)
@ -195,28 +192,28 @@ subroutine grid_mechanical_FEM_init
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMSNESSetJacobianLocal(mechanical_grid,formJacobian,PETSC_NULL_SNES,err_PETSc) call DMSNESSetJacobianLocal(mechanical_grid,formJacobian,PETSC_NULL_SNES,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetConvergenceTest(mechanical_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "_converged" call SNESSetConvergenceTest(SNES_mechanical,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "_converged"
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1_pPetscInt), err_PETSc) ! ignore linear solve failures call SNESSetMaxLinearSolveFailures(SNES_mechanical, huge(1_pPetscInt), err_PETSc) ! ignore linear solve failures
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetFromOptions(mechanical_snes,err_PETSc) ! pull it all together with additional cli arguments call SNESSetDM(SNES_mechanical,mechanical_grid,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetFromOptions(SNES_mechanical,err_PETSc) ! pull it all together with additional cli arguments
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call VecSet(solution_current,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc) call VecSet(solution_current,0.0_pReal,err_PETSc)
call VecSet(solution_lastInc,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecSet(solution_rate ,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc) call VecSet(solution_lastInc,0.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
call VecSet(solution_rate ,0.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,err_PETSc) call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,err_PETSc) call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMDAGetCorners(mechanical_grid,xstart,ystart,zstart,xend,yend,zend,err_PETSc) ! local grid extent
CHKERRQ(err_PETSc)
xend = xstart+xend-1
yend = ystart+yend-1
zend = zstart+zend-1
delta = geomSize/real(grid,pReal) ! grid spacing delta = geomSize/real(grid,pReal) ! grid spacing
detJ = product(delta) ! cell volume detJ = product(delta) ! cell volume
@ -311,14 +308,9 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution)
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
!-------------------------------------------------------------------------------------------------- call SNESsolve(SNES_mechanical,PETSC_NULL_VEC,solution_current,err_PETSc)
! solve BVP
call SNESsolve(mechanical_snes,PETSC_NULL_VEC,solution_current,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESGetConvergedReason(SNES_mechanical,reason,err_PETSc)
!--------------------------------------------------------------------------------------------------
! check convergence
call SNESGetConvergedReason(mechanical_snes,reason,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
solution%converged = reason > 0 solution%converged = reason > 0
@ -386,9 +378,11 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
call VecScale(solution_rate,1.0_pReal/Delta_t_old,err_PETSc) call VecScale(solution_rate,1.0_pReal/Delta_t_old,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
else else
call VecSet(solution_rate,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc) call VecSet(solution_rate,0.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
endif endif
call VecCopy(solution_current,solution_lastInc,err_PETSc); CHKERRQ(err_PETSc) call VecCopy(solution_current,solution_lastInc,err_PETSc)
CHKERRQ(err_PETSc)
F_lastInc = F F_lastInc = F
@ -515,6 +509,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,e
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
print'(/,1x,a)', '===========================================================================' print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
err_PETSc = 0
end subroutine converged end subroutine converged
@ -527,7 +522,7 @@ subroutine formResidual(da_local,x_local, &
DM :: da_local DM :: da_local
Vec :: x_local, f_local Vec :: x_local, f_local
PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, r
PetscScalar, dimension(8,3) :: x_elem, f_elem PetscScalar, dimension(8,3) :: x_elem, f_elem
PetscInt :: i, ii, j, jj, k, kk, ctr, ele PetscInt :: i, ii, j, jj, k, kk, ctr, ele
PetscInt :: & PetscInt :: &
@ -538,9 +533,9 @@ subroutine formResidual(da_local,x_local, &
integer(MPI_INTEGER_KIND) :: err_MPI integer(MPI_INTEGER_KIND) :: err_MPI
real(pReal), dimension(3,3,3,3) :: devNull real(pReal), dimension(3,3,3,3) :: devNull
call SNESGetNumberFunctionEvals(mechanical_snes,nfuncs,err_PETSc) call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESGetIterationNumber(mechanical_snes,PETScIter,err_PETSc) call SNESGetIterationNumber(SNES_mechanical,PETScIter,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
@ -559,17 +554,18 @@ subroutine formResidual(da_local,x_local, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! get deformation gradient ! get deformation gradient
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc);CHKERRQ(err_PETSc) call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc)
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend CHKERRQ(err_PETSc)
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
ctr = 0 ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1 ctr = ctr + 1
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
enddo; enddo; enddo enddo; enddo; enddo
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 F(1:3,1:3,i,j,k-grid3offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
enddo; enddo; enddo enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc);CHKERRQ(err_PETSc) call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
@ -586,47 +582,53 @@ subroutine formResidual(da_local,x_local, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
call VecSet(f_local,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc) call VecSet(f_local,0.0_pReal,err_PETSc)
call DMDAVecGetArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc);CHKERRQ(err_PETSc) call DMDAVecGetArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
ele = 0 ele = 0
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
ctr = 0 ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1 ctr = ctr + 1
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
enddo; enddo; enddo enddo; enddo; enddo
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
ele = ele + 1 ele = ele + 1
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-grid3offset)))*detJ + &
matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + & matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,ele) + & homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
ctr = 0 ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1 ctr = ctr + 1
f_scal(0:2,i+ii,j+jj,k+kk) = f_scal(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3) r(0:2,i+ii,j+jj,k+kk) = r(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3)
enddo; enddo; enddo enddo; enddo; enddo
enddo; enddo; enddo enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc);CHKERRQ(err_PETSc) call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc)
call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! applying boundary conditions ! applying boundary conditions
call DMDAVecGetArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc) call DMDAVecGetArrayF90(da_local,f_local,r,err_PETSc)
if (zstart == 0) then CHKERRQ(err_PETSc)
f_scal(0:2,xstart,ystart,zstart) = 0.0 if (grid3offset == 0) then
f_scal(0:2,xend+1,ystart,zstart) = 0.0 r(0:2,0, 0, 0) = 0.0_pReal
f_scal(0:2,xstart,yend+1,zstart) = 0.0 r(0:2,grid(1),0, 0) = 0.0_pReal
f_scal(0:2,xend+1,yend+1,zstart) = 0.0 r(0:2,0, grid(2),0) = 0.0_pReal
endif r(0:2,grid(1),grid(2),0) = 0.0_pReal
if (zend + 1 == grid(3)) then end if
f_scal(0:2,xstart,ystart,zend+1) = 0.0 if (grid3+grid3offset == grid(3)) then
f_scal(0:2,xend+1,ystart,zend+1) = 0.0 r(0:2,0, 0, grid(3)) = 0.0_pReal
f_scal(0:2,xstart,yend+1,zend+1) = 0.0 r(0:2,grid(1),0, grid(3)) = 0.0_pReal
f_scal(0:2,xend+1,yend+1,zend+1) = 0.0 r(0:2,0, grid(2),grid(3)) = 0.0_pReal
endif r(0:2,grid(1),grid(2),grid(3)) = 0.0_pReal
call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc) end if
call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
end subroutine formResidual end subroutine formResidual
@ -643,7 +645,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
PetscScalar,pointer,dimension(:,:,:,:) :: x_scal PetscScalar,pointer,dimension(:,:,:,:) :: x_scal
PetscScalar,dimension(24,24) :: K_ele PetscScalar,dimension(24,24) :: K_ele
PetscScalar,dimension(9,24) :: BMatFull PetscScalar,dimension(9,24) :: BMatFull
PetscInt :: i, ii, j, jj, k, kk, ctr, ele PetscInt :: i, ii, j, jj, k, kk, ctr, ce
PetscInt,dimension(3),parameter :: rows = [0, 1, 2] PetscInt,dimension(3),parameter :: rows = [0, 1, 2]
PetscScalar :: diag PetscScalar :: diag
PetscObject :: dummy PetscObject :: dummy
@ -658,11 +660,12 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,err_PETSc) call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call MatZeroEntries(Jac,err_PETSc); CHKERRQ(err_PETSc) call MatZeroEntries(Jac,err_PETSc)
ele = 0 CHKERRQ(err_PETSc)
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend ce = 0
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
ctr = 0 ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1 ctr = ctr + 1
col(MatStencil_i,ctr ) = i+ii col(MatStencil_i,ctr ) = i+ii
col(MatStencil_j,ctr ) = j+jj col(MatStencil_j,ctr ) = j+jj
@ -678,44 +681,47 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
col(MatStencil_c,ctr+16) = 2 col(MatStencil_c,ctr+16) = 2
enddo; enddo; enddo enddo; enddo; enddo
row = col row = col
ele = ele + 1 ce = ce + 1
K_ele = 0.0 K_ele = 0.0
K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + & K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
homogenization_dPdF(2,2,2,2,ele) + & homogenization_dPdF(2,2,2,2,ce) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,ce))/3.0_pReal
K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + & K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
homogenization_dPdF(2,2,2,2,ele) + & homogenization_dPdF(2,2,2,2,ce) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,ce))/3.0_pReal
K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + & K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
homogenization_dPdF(2,2,2,2,ele) + & homogenization_dPdF(2,2,2,2,ce) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,ce))/3.0_pReal
K_ele = K_ele + & K_ele = K_ele + &
matmul(transpose(BMatFull), & matmul(transpose(BMatFull), &
matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,ele), & matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,ce), &
shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ
call MatSetValuesStencil(Jac,24_pPETScInt,row,24_pPetscInt,col,K_ele,ADD_VALUES,err_PETSc) call MatSetValuesStencil(Jac,24_pPETScInt,row,24_pPetscInt,col,K_ele,ADD_VALUES,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
enddo; enddo; enddo enddo; enddo; enddo
call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc) call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc) call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc)
CHKERRQ(err_PETSc)
call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc)
CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! applying boundary conditions ! applying boundary conditions
diag = (C_volAvg(1,1,1,1)/delta(1)**2 + & diag = (C_volAvg(1,1,1,1)/delta(1)**2 + C_volAvg(2,2,2,2)/delta(2)**2 + C_volAvg(3,3,3,3)/delta(3)**2) &
C_volAvg(2,2,2,2)/delta(2)**2 + & * detJ
C_volAvg(3,3,3,3)/delta(3)**2)*detJ
call MatZeroRowsColumns(Jac,size(rows,kind=pPetscInt),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,err_PETSc) call MatZeroRowsColumns(Jac,size(rows,kind=pPetscInt),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMGetGlobalVector(da_local,coordinates,err_PETSc) call DMGetGlobalVector(da_local,coordinates,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da_local,coordinates,x_scal,err_PETSc) call DMDAVecGetArrayF90(da_local,coordinates,x_scal,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
ele = 0 ce = 0
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
ele = ele + 1 ce = ce + 1
x_scal(0:2,i,j,k) = discretization_IPcoords(1:3,ele) x_scal(0:2,i-1,j-1,k-1) = discretization_IPcoords(1:3,ce)
enddo; enddo; enddo enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,err_PETSc) call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,err_PETSc)
CHKERRQ(err_PETSc) ! initialize to undeformed coordinates (ToDo: use ip coordinates) CHKERRQ(err_PETSc) ! initialize to undeformed coordinates (ToDo: use ip coordinates)

77
src/grid/grid_mech_spectral_basic.f90
View File

@ -50,7 +50,7 @@ module grid_mechanical_spectral_basic
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM :: da DM :: da
SNES :: snes SNES :: SNES_mechanical
Vec :: solution_vec Vec :: solution_vec
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -158,8 +158,10 @@ subroutine grid_mechanical_spectral_basic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,err_PETSc); CHKERRQ(err_PETSc) call SNESCreate(PETSC_COMM_WORLD,SNES_mechanical,err_PETSc)
call SNESSetOptionsPrefix(snes,'mechanical_',err_PETSc);CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt) localK(worldrank) = int(grid3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@ -173,19 +175,25 @@ subroutine grid_mechanical_spectral_basic_init
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
da,err_PETSc) ! handle, error da,err_PETSc) ! handle, error
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetDM(snes,da,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da call DMsetFromOptions(da,err_PETSc)
call DMsetFromOptions(da,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMsetUp(da,err_PETSc); CHKERRQ(err_PETSc) call DMsetUp(da,err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 9, i.e. every def grad tensor) CHKERRQ(err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 9, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "converged" call SNESsetConvergenceTest(SNES_mechanical,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "converged"
CHKERRQ(err_PETSc)
call SNESSetDM(SNES_mechanical,da,err_PETSc)
CHKERRQ(err_PETSc)
call SNESsetFromOptions(SNES_mechanical,err_PETSc) ! pull it all together with additional CLI arguments
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESsetFromOptions(snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc) ! places pointer on PETSc data call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc) ! places pointer on PETSc data
CHKERRQ(err_PETSc)
restartRead: if (interface_restartInc > 0) then restartRead: if (interface_restartInc > 0) then
print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file' print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
@ -218,7 +226,8 @@ subroutine grid_mechanical_spectral_basic_init
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal) ! time increment 0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc) ! deassociate pointer call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc) ! deassociate pointer
CHKERRQ(err_PETSc)
restartRead2: if (interface_restartInc > 0) then restartRead2: if (interface_restartInc > 0) then
print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file' print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
@ -270,13 +279,10 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg) if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
!-------------------------------------------------------------------------------------------------- call SNESsolve(SNES_mechanical,PETSC_NULL_VEC,solution_vec,err_PETSc)
! solve BVP CHKERRQ(err_PETSc)
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,err_PETSc); CHKERRQ(err_PETSc) call SNESGetConvergedReason(SNES_mechanical,reason,err_PETSc)
CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! check convergence
call SNESGetConvergedReason(snes,reason,err_PETSc); CHKERRQ(err_PETSc)
solution%converged = reason > 0 solution%converged = reason > 0
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
@ -310,7 +316,8 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
PetscScalar, pointer, dimension(:,:,:,:) :: F PetscScalar, pointer, dimension(:,:,:,:) :: F
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc) call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
if (cutBack) then if (cutBack) then
C_volAvg = C_volAvgLastInc C_volAvg = C_volAvgLastInc
@ -353,7 +360,8 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc) call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide available data ! set module wide available data
@ -372,9 +380,11 @@ subroutine grid_mechanical_spectral_basic_updateCoords
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
PetscScalar, dimension(:,:,:,:), pointer :: F PetscScalar, dimension(:,:,:,:), pointer :: F
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc) call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
call utilities_updateCoords(F) call utilities_updateCoords(F)
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc) call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
end subroutine grid_mechanical_spectral_basic_updateCoords end subroutine grid_mechanical_spectral_basic_updateCoords
@ -388,7 +398,8 @@ subroutine grid_mechanical_spectral_basic_restartWrite
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
PetscScalar, dimension(:,:,:,:), pointer :: F PetscScalar, dimension(:,:,:,:), pointer :: F
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc) call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT) print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
@ -415,7 +426,8 @@ subroutine grid_mechanical_spectral_basic_restartWrite
if (num%update_gamma) call utilities_saveReferenceStiffness if (num%update_gamma) call utilities_saveReferenceStiffness
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc) call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
end subroutine grid_mechanical_spectral_basic_restartWrite end subroutine grid_mechanical_spectral_basic_restartWrite
@ -457,6 +469,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
print'(/,1x,a)', '===========================================================================' print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
err_PETSc = 0
end subroutine converged end subroutine converged
@ -465,13 +478,13 @@ end subroutine converged
!> @brief forms the residual vector !> @brief forms the residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine formResidual(in, F, & subroutine formResidual(in, F, &
residuum, dummy, err_PETSc) r, dummy, err_PETSc)
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work)
PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), & PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), &
intent(in) :: F !< deformation gradient field intent(in) :: F !< deformation gradient field
PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), & PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), &
intent(out) :: residuum !< residuum field intent(out) :: r !< residuum field
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
deltaF_aim deltaF_aim
PetscInt :: & PetscInt :: &
@ -481,8 +494,10 @@ subroutine formResidual(in, F, &
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
integer(MPI_INTEGER_KIND) :: err_MPI integer(MPI_INTEGER_KIND) :: err_MPI
call SNESGetNumberFunctionEvals(snes,nfuncs,err_PETSc); CHKERRQ(err_PETSc) call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc)
call SNESGetIterationNumber(snes,PETScIter,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESGetIterationNumber(SNES_mechanical,PETScIter,err_PETSc)
CHKERRQ(err_PETSc)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
@ -500,7 +515,7 @@ subroutine formResidual(in, F, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) call utilities_constitutiveResponse(r, & ! residuum gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, & P_av,C_volAvg,C_minMaxAvg, &
F,params%Delta_t,params%rotation_BC) F,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
@ -515,7 +530,7 @@ subroutine formResidual(in, F, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme ! updated deformation gradient using fix point algorithm of basic scheme
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r ! store fPK field for subsequent FFT forward transform
call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" call utilities_FFTtensorForward ! FFT forward of global "tensorField_real"
err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier
@ -523,7 +538,7 @@ subroutine formResidual(in, F, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too r = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
end subroutine formResidual end subroutine formResidual

78
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -55,7 +55,7 @@ module grid_mechanical_spectral_polarisation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM :: da DM :: da
SNES :: snes SNES :: SNES_mechanical
Vec :: solution_vec Vec :: solution_vec
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -178,8 +178,10 @@ subroutine grid_mechanical_spectral_polarisation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,err_PETSc); CHKERRQ(err_PETSc) call SNESCreate(PETSC_COMM_WORLD,SNES_mechanical,err_PETSc)
call SNESSetOptionsPrefix(snes,'mechanical_',err_PETSc);CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt) localK(worldrank) = int(grid3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@ -193,19 +195,25 @@ subroutine grid_mechanical_spectral_polarisation_init
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
da,err_PETSc) ! handle, error da,err_PETSc) ! handle, error
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetDM(snes,da,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da call DMsetFromOptions(da,err_PETSc)
call DMsetFromOptions(da,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMsetUp(da,err_PETSc); CHKERRQ(err_PETSc) call DMsetUp(da,err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 18, i.e. every def grad tensor) CHKERRQ(err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 18, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "converged" call SNESsetConvergenceTest(SNES_mechanical,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "converged"
CHKERRQ(err_PETSc)
call SNESSetDM(SNES_mechanical,da,err_PETSc)
CHKERRQ(err_PETSc)
call SNESsetFromOptions(SNES_mechanical,err_PETSc) ! pull it all together with additional CLI arguments
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESsetFromOptions(snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc); CHKERRQ(err_PETSc) ! places pointer on PETSc data call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc) ! places pointer on PETSc data
CHKERRQ(err_PETSc)
F => FandF_tau(0: 8,:,:,:) F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:) F_tau => FandF_tau(9:17,:,:,:)
@ -303,13 +311,10 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
S_scale = math_invSym3333(C_minMaxAvg) S_scale = math_invSym3333(C_minMaxAvg)
end if end if
!-------------------------------------------------------------------------------------------------- call SNESSolve(SNES_mechanical,PETSC_NULL_VEC,solution_vec,err_PETSc)
! solve BVP CHKERRQ(err_PETSc)
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,err_PETSc); CHKERRQ(err_PETSc) call SNESGetConvergedReason(SNES_mechanical,reason,err_PETSc)
CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! check convergence
call SNESGetConvergedReason(snes,reason,err_PETSc); CHKERRQ(err_PETSc)
solution%converged = reason > 0 solution%converged = reason > 0
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
@ -345,7 +350,8 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
real(pReal), dimension(3,3) :: F_lambda33 real(pReal), dimension(3,3) :: F_lambda33
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc); CHKERRQ(err_PETSc) call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc)
CHKERRQ(err_PETSc)
F => FandF_tau(0: 8,:,:,:) F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:) F_tau => FandF_tau(9:17,:,:,:)
@ -446,7 +452,8 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc); CHKERRQ(err_PETSc) call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc)
CHKERRQ(err_PETSc)
F => FandF_tau(0: 8,:,:,:) F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:) F_tau => FandF_tau(9:17,:,:,:)
@ -523,6 +530,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
print'(/,1x,a)', '===========================================================================' print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
err_PETSc = 0
end subroutine converged end subroutine converged
@ -531,18 +539,18 @@ end subroutine converged
!> @brief forms the residual vector !> @brief forms the residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine formResidual(in, FandF_tau, & subroutine formResidual(in, FandF_tau, &
residuum, dummy,err_PETSc) r, dummy,err_PETSc)
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work)
PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), & PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), &
target, intent(in) :: FandF_tau target, intent(in) :: FandF_tau
PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),& PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),&
target, intent(out) :: residuum !< residuum field target, intent(out) :: r !< residuum field
PetscScalar, pointer, dimension(:,:,:,:,:) :: & PetscScalar, pointer, dimension(:,:,:,:,:) :: &
F, & F, &
F_tau, & F_tau, &
residual_F, & r_F, &
residual_F_tau r_F_tau
PetscInt :: & PetscInt :: &
PETScIter, & PETScIter, &
nfuncs nfuncs
@ -558,17 +566,19 @@ subroutine formResidual(in, FandF_tau, &
XG_RANGE,YG_RANGE,ZG_RANGE) XG_RANGE,YG_RANGE,ZG_RANGE)
F_tau => FandF_tau(1:3,1:3,2,& F_tau => FandF_tau(1:3,1:3,2,&
XG_RANGE,YG_RANGE,ZG_RANGE) XG_RANGE,YG_RANGE,ZG_RANGE)
residual_F => residuum(1:3,1:3,1,& r_F => r(1:3,1:3,1,&
X_RANGE, Y_RANGE, Z_RANGE) X_RANGE, Y_RANGE, Z_RANGE)
residual_F_tau => residuum(1:3,1:3,2,& r_F_tau => r(1:3,1:3,2,&
X_RANGE, Y_RANGE, Z_RANGE) X_RANGE, Y_RANGE, Z_RANGE)
F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,F_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) call MPI_Allreduce(MPI_IN_PLACE,F_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call SNESGetNumberFunctionEvals(snes,nfuncs,err_PETSc); CHKERRQ(err_PETSc) call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc)
call SNESGetIterationNumber(snes,PETScIter,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESGetIterationNumber(SNES_mechanical,PETScIter,err_PETSc)
CHKERRQ(err_PETSc)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
@ -602,13 +612,13 @@ subroutine formResidual(in, FandF_tau, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
residual_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) call utilities_constitutiveResponse(r_F, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, & P_av,C_volAvg,C_minMaxAvg, &
F - residual_F_tau/num%beta,params%Delta_t,params%rotation_BC) F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -619,7 +629,7 @@ subroutine formResidual(in, FandF_tau, &
params%stress_mask))) params%stress_mask)))
! calculate divergence ! calculate divergence
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r_F !< stress field in disguise
call utilities_FFTtensorForward call utilities_FFTtensorForward
err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress
@ -628,11 +638,11 @@ subroutine formResidual(in, FandF_tau, &
e = 0 e = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
e = e + 1 e = e + 1
residual_F(1:3,1:3,i,j,k) = & r_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), & math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), &
residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & r_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
+ residual_F_tau(1:3,1:3,i,j,k) + r_F_tau(1:3,1:3,i,j,k)
end do; end do; end do end do; end do; end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

85
src/grid/grid_thermal_spectral.f90
View File

@ -38,9 +38,8 @@ module grid_thermal_spectral
type(tSolutionParams) :: params type(tSolutionParams) :: params
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
SNES :: thermal_snes SNES :: SNES_thermal
Vec :: solution_vec Vec :: solution_vec
PetscInt :: xstart, xend, ystart, yend, zstart, zend
real(pReal), dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
T_current, & !< field of current temperature T_current, & !< field of current temperature
T_lastInc, & !< field of previous temperature T_lastInc, & !< field of previous temperature
@ -100,10 +99,24 @@ subroutine grid_thermal_spectral_init(T_0)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! init fields
allocate(T_current(grid(1),grid(2),grid3), source=T_0)
allocate(T_lastInc(grid(1),grid(2),grid3), source=T_0)
allocate(T_stagInc(grid(1),grid(2),grid3), source=T_0)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1
call homogenization_thermal_setField(T_0,0.0_pReal,ce)
end do; end do; end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,thermal_snes,err_PETSc); CHKERRQ(err_PETSc) call SNESCreate(PETSC_COMM_WORLD,SNES_thermal,err_PETSc)
call SNESSetOptionsPrefix(thermal_snes,'thermal_',err_PETSc);CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(SNES_thermal,'thermal_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt) localK(worldrank) = int(grid3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@ -117,42 +130,25 @@ subroutine grid_thermal_spectral_init(T_0)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
thermal_grid,err_PETSc) ! handle, error thermal_grid,err_PETSc) ! handle, error
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetDM(thermal_snes,thermal_grid,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da call DMsetFromOptions(thermal_grid,err_PETSc)
call DMsetFromOptions(thermal_grid,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMsetUp(thermal_grid,err_PETSc); CHKERRQ(err_PETSc) call DMsetUp(thermal_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor) call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESSetFromOptions(thermal_snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments call SNESSetDM(SNES_thermal,thermal_grid,err_PETSc)
CHKERRQ(err_PETSc)
!-------------------------------------------------------------------------------------------------- call SNESSetFromOptions(SNES_thermal,err_PETSc) ! pull it all together with additional CLI arguments
! init fields
call DMDAGetCorners(thermal_grid,xstart,ystart,zstart,xend,yend,zend,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
xend = xstart + xend - 1
yend = ystart + yend - 1
zend = zstart + zend - 1
allocate(T_current(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1
T_current(i,j,k) = T_0
T_lastInc(i,j,k) = T_current(i,j,k)
T_stagInc(i,j,k) = T_current(i,j,k)
call homogenization_thermal_setField(T_0,0.0_pReal,ce)
end do; end do; end do
call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc) call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
T_PETSc(xstart:xend,ystart:yend,zstart:zend) = T_current T_PETSc = T_current
call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc) call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call updateReference call updateReference()
end subroutine grid_thermal_spectral_init end subroutine grid_thermal_spectral_init
@ -179,9 +175,9 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
! set module wide availabe data ! set module wide availabe data
params%Delta_t = Delta_t params%Delta_t = Delta_t
call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,err_PETSc) call SNESSolve(SNES_thermal,PETSC_NULL_VEC,solution_vec,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESGetConvergedReason(thermal_snes,reason,err_PETSc) call SNESGetConvergedReason(SNES_thermal,reason,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
if (reason < 1) then if (reason < 1) then
@ -207,8 +203,10 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
end do; end do; end do end do; end do; end do
call VecMin(solution_vec,devNull,T_min,err_PETSc); CHKERRQ(err_PETSc) call VecMin(solution_vec,devNull,T_min,err_PETSc)
call VecMax(solution_vec,devNull,T_max,err_PETSc); CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecMax(solution_vec,devNull,T_max,err_PETSc)
CHKERRQ(err_PETSc)
if (solution%converged) & if (solution%converged) &
print'(/,1x,a)', '... thermal conduction converged ..................................' print'(/,1x,a)', '... thermal conduction converged ..................................'
print'(/,1x,a,f8.4,2x,f8.4,2x,f8.4)', 'Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm print'(/,1x,a,f8.4,2x,f8.4,2x,f8.4)', 'Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm
@ -226,7 +224,7 @@ subroutine grid_thermal_spectral_forward(cutBack)
logical, intent(in) :: cutBack logical, intent(in) :: cutBack
integer :: i, j, k, ce integer :: i, j, k, ce
DM :: dm_local DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: T_PETSc
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
if (cutBack) then if (cutBack) then
@ -235,12 +233,12 @@ subroutine grid_thermal_spectral_forward(cutBack)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reverting thermal field state ! reverting thermal field state
call SNESGetDM(thermal_snes,dm_local,err_PETSc) call SNESGetDM(SNES_thermal,dm_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,err_PETSc) !< get the data out of PETSc to work with call DMDAVecGetArrayF90(dm_local,solution_vec,T_PETSc,err_PETSc) !< get the data out of PETSc to work with
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current T_PETSc = T_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,err_PETSc) call DMDAVecRestoreArrayF90(dm_local,solution_vec,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
@ -258,7 +256,7 @@ end subroutine grid_thermal_spectral_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forms the spectral thermal residual vector !> @brief forms the spectral thermal residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err) subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in in
@ -267,9 +265,9 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
x_scal x_scal
PetscScalar, dimension( & PetscScalar, dimension( &
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
f_scal r
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: dummy_err PetscErrorCode :: err_PETSc
integer :: i, j, k, ce integer :: i, j, k, ce
T_current = x_scal T_current = x_scal
@ -304,7 +302,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
f_scal = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) r = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3)
err_PETSc = 0
end subroutine formResidual end subroutine formResidual

60
src/grid/spectral_utilities.f90
View File

@ -156,7 +156,7 @@ subroutine spectral_utilities_init
integer(C_INTPTR_T) :: alloc_local, local_K, local_K_offset integer(C_INTPTR_T) :: alloc_local, local_K, local_K_offset
integer(C_INTPTR_T), parameter :: & integer(C_INTPTR_T), parameter :: &
scalarSize = 1_C_INTPTR_T, & scalarSize = 1_C_INTPTR_T, &
vecSize = 3_C_INTPTR_T, & vectorSize = 3_C_INTPTR_T, &
tensorSize = 9_C_INTPTR_T tensorSize = 9_C_INTPTR_T
character(len=*), parameter :: & character(len=*), parameter :: &
PETSCDEBUG = ' -snes_view -snes_monitor ' PETSCDEBUG = ' -snes_view -snes_monitor '
@ -274,7 +274,7 @@ subroutine spectral_utilities_init
call c_f_pointer(tensorField, tensorField_fourier, [3_C_INTPTR_T,3_C_INTPTR_T, & call c_f_pointer(tensorField, tensorField_fourier, [3_C_INTPTR_T,3_C_INTPTR_T, &
gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T , gridFFTW(2),local_K]) ! place a pointer for a fourier tensor representation gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T , gridFFTW(2),local_K]) ! place a pointer for a fourier tensor representation
vectorField = fftw_alloc_complex(vecSize*alloc_local) vectorField = fftw_alloc_complex(vectorSize*alloc_local)
call c_f_pointer(vectorField, vectorField_real, [3_C_INTPTR_T,& call c_f_pointer(vectorField, vectorField_real, [3_C_INTPTR_T,&
2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K]) ! place a pointer for a real vector representation 2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K]) ! place a pointer for a real vector representation
call c_f_pointer(vectorField, vectorField_fourier,[3_C_INTPTR_T,& call c_f_pointer(vectorField, vectorField_fourier,[3_C_INTPTR_T,&
@ -288,42 +288,42 @@ subroutine spectral_utilities_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! tensor MPI fftw plans ! tensor MPI fftw plans
planTensorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order planTensorForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),tensorSize, &
tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
tensorField_real, tensorField_fourier, & ! input data, output data tensorField_real,tensorField_fourier, &
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision PETSC_COMM_WORLD,FFTW_planner_flag)
if (.not. c_associated(planTensorForth)) error stop 'FFTW error' if (.not. c_associated(planTensorForth)) error stop 'FFTW error'
planTensorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order planTensorBack = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),tensorSize, &
tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
tensorField_fourier,tensorField_real, & ! input data, output data tensorField_fourier,tensorField_real, &
PETSC_COMM_WORLD, FFTW_planner_flag) ! all processors, planer precision PETSC_COMM_WORLD, FFTW_planner_flag)
if (.not. c_associated(planTensorBack)) error stop 'FFTW error' if (.not. c_associated(planTensorBack)) error stop 'FFTW error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! vector MPI fftw plans ! vector MPI fftw plans
planVectorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order planVectorForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),vectorSize, &
vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,&! no. of transforms, default iblock and oblock FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
vectorField_real, vectorField_fourier, & ! input data, output data vectorField_real,vectorField_fourier, &
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision PETSC_COMM_WORLD,FFTW_planner_flag)
if (.not. C_ASSOCIATED(planVectorForth)) error stop 'FFTW error' if (.not. c_associated(planVectorForth)) error stop 'FFTW error'
planVectorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order planVectorBack = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),vectorSize, &
vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
vectorField_fourier,vectorField_real, & ! input data, output data vectorField_fourier,vectorField_real, &
PETSC_COMM_WORLD, FFTW_planner_flag) ! all processors, planer precision PETSC_COMM_WORLD, FFTW_planner_flag)
if (.not. C_ASSOCIATED(planVectorBack)) error stop 'FFTW error' if (.not. c_associated(planVectorBack)) error stop 'FFTW error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! scalar MPI fftw plans ! scalar MPI fftw plans
planScalarForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order planScalarForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),scalarSize, &
scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
scalarField_real, scalarField_fourier, & ! input data, output data scalarField_real,scalarField_fourier, &
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision PETSC_COMM_WORLD,FFTW_planner_flag)
if (.not. C_ASSOCIATED(planScalarForth)) error stop 'FFTW error' if (.not. c_associated(planScalarForth)) error stop 'FFTW error'
planScalarBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order, no. of transforms planScalarBack = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),scalarSize, &
scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
scalarField_fourier,scalarField_real, & ! input data, output data scalarField_fourier,scalarField_real, &
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision PETSC_COMM_WORLD, FFTW_planner_flag)
if (.not. C_ASSOCIATED(planScalarBack)) error stop 'FFTW error' if (.not. c_associated(planScalarBack)) error stop 'FFTW error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)

142
src/lattice.f90
View File

@ -587,8 +587,8 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Non-schmid projections for bcc with up to 6 coefficients !> @brief Non-schmid projections for bcc with up to 6 coefficients
! Koester et al. 2012, Acta Materialia 60 (2012) 3894–3901, eq. (17) ! https://doi.org/10.1016/j.actamat.2012.03.053, eq. (17)
! Gröger et al. 2008, Acta Materialia 56 (2008) 5412–5425, table 1 ! https://doi.org/10.1016/j.actamat.2008.07.037, table 1
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix) function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix)
@ -602,6 +602,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
type(rotation) :: R type(rotation) :: R
integer :: i integer :: i
if (abs(sense) /= 1) error stop 'Sense in lattice_nonSchmidMatrix' if (abs(sense) /= 1) error stop 'Sense in lattice_nonSchmidMatrix'
coordinateSystem = buildCoordinateSystem(Nslip,BCC_NSLIPSYSTEM,BCC_SYSTEMSLIP,'cI',0.0_pReal) coordinateSystem = buildCoordinateSystem(Nslip,BCC_NSLIPSYSTEM,BCC_SYSTEMSLIP,'cI',0.0_pReal)
@ -634,7 +635,9 @@ end function lattice_nonSchmidMatrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Slip-slip interaction matrix !> @brief Slip-slip interaction matrix
!> details only active slip systems are considered !> @details only active slip systems are considered
!> @details https://doi.org/10.1016/j.actamat.2016.12.040 (fcc: Tab S4-1, bcc: Tab S5-1)
!> @details https://doi.org/10.1016/j.ijplas.2014.06.010 (hex: Tab 3b)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix) function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix)
@ -646,6 +649,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
integer, dimension(:), allocatable :: NslipMax integer, dimension(:), allocatable :: NslipMax
integer, dimension(:,:), allocatable :: interactionTypes integer, dimension(:,:), allocatable :: interactionTypes
integer, dimension(FCC_NSLIP,FCC_NSLIP), parameter :: & integer, dimension(FCC_NSLIP,FCC_NSLIP), parameter :: &
FCC_INTERACTIONSLIPSLIP = reshape( [& FCC_INTERACTIONSLIPSLIP = reshape( [&
1, 2, 2, 4, 7, 5, 3, 5, 5, 4, 6, 7, 10,11,10,11,12,13, & ! -----> acting (forest) 1, 2, 2, 4, 7, 5, 3, 5, 5, 4, 6, 7, 10,11,10,11,12,13, & ! -----> acting (forest)
@ -750,41 +754,113 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: & integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: &
HEX_INTERACTIONSLIPSLIP = reshape( [& HEX_INTERACTIONSLIPSLIP = reshape( [&
! basal prism 1. pyr<a> 1. pyr<c+a> 2. pyr<c+a> ! basal prism 1. pyr<a> 1. pyr<c+a> 2. pyr<c+a>
1, 2, 2, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13,13,13,13,13,13,13, 21,21,21,21,21,21, & ! -----> acting (forest) 1, 2, 2, 3, 4, 4, 9,10, 9, 9,10, 9, 20,21,22,22,21,20,20,21,22,22,21,20, 47,47,48,47,47,48, & ! -----> acting (forest)
2, 1, 2, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13,13,13,13,13,13,13, 21,21,21,21,21,21, & ! | basal 2, 1, 2, 4, 3, 4, 10, 9, 9,10, 9, 9, 22,22,21,20,20,21,22,22,21,20,20,21, 47,48,47,47,48,47, & ! | basal
2, 2, 1, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13,13,13,13,13,13,13, 21,21,21,21,21,21, & ! | 2, 2, 1, 4, 4, 3, 9, 9,10, 9, 9,10, 21,20,20,21,22,22,21,20,20,21,22,22, 48,47,47,48,47,47, & ! |
! v ! v
6, 6, 6, 4, 5, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14,14,14,14,14,14,14, 22,22,22,22,22,22, & ! reacting (primary) 7, 8, 8, 5, 6, 6, 11,12,11,11,12,11, 23,24,25,25,24,23,23,24,25,25,24,23, 49,49,50,49,49,50, & ! reacting (primary)
6, 6, 6, 5, 4, 5, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14,14,14,14,14,14,14, 22,22,22,22,22,22, & ! prism 8, 7, 8, 6, 5, 6, 12,11,11,12,11,11, 25,25,24,23,23,24,25,25,24,23,23,24, 49,50,49,49,50,49, & ! prism
6, 6, 6, 5, 5, 4, 8, 8, 8, 8, 8, 8, 14,14,14,14,14,14,14,14,14,14,14,14, 22,22,22,22,22,22, & 8, 8, 7, 6, 6, 5, 11,11,12,11,11,12, 24,23,23,24,25,25,24,23,23,24,25,25, 50,49,49,50,49,49, &
12,12,12, 11,11,11, 9,10,10,10,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & 18,19,18, 16,17,16, 13,14,14,15,14,14, 26,26,27,28,28,27,29,29,27,28,28,27, 51,52,51,51,52,51, &
12,12,12, 11,11,11, 10, 9,10,10,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & 19,18,18, 17,16,16, 14,13,14,14,15,14, 28,27,26,26,27,28,28,27,29,29,27,28, 51,51,52,51,51,52, &
12,12,12, 11,11,11, 10,10, 9,10,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & 18,18,19, 16,16,17, 14,14,13,14,14,15, 27,28,28,27,26,26,27,28,28,27,29,29, 52,51,51,52,51,51, &
12,12,12, 11,11,11, 10,10,10, 9,10,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & ! 1. pyr<a> 18,19,18, 16,17,16, 15,14,14,13,14,14, 29,29,27,28,28,27,26,26,27,28,28,27, 51,52,51,51,52,51, & ! 1. pyr<a>
12,12,12, 11,11,11, 10,10,10,10, 9,10, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & 19,18,18, 17,16,16, 14,15,14,14,13,14, 28,27,29,29,27,28,28,27,26,26,27,28, 51,51,52,51,51,52, &
12,12,12, 11,11,11, 10,10,10,10,10, 9, 15,15,15,15,15,15,15,15,15,15,15,15, 23,23,23,23,23,23, & 18,18,19, 16,16,17, 14,14,15,14,14,13, 27,28,28,27,29,29,27,28,28,27,26,26, 52,51,51,52,51,51, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 16,17,17,17,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, & 44,45,46, 41,42,43, 37,38,39,40,38,39, 30,31,32,32,32,33,34,35,32,32,32,36, 53,54,55,53,54,56, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,16,17,17,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, & 46,45,44, 43,42,41, 37,39,38,40,39,38, 31,30,36,32,32,32,35,34,33,32,32,32, 56,54,53,55,54,53, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,16,17,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, & 45,46,44, 42,43,41, 39,37,38,39,40,38, 32,36,30,31,32,32,32,33,34,35,32,32, 56,53,54,55,53,54, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,16,17,17,17,17,17,17,17,17, 24,24,24,24,24,24, & 45,44,46, 42,41,43, 38,37,39,38,40,39, 32,32,31,30,36,32,32,32,35,34,33,32, 53,56,54,53,55,54, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,16,17,17,17,17,17,17,17, 24,24,24,24,24,24, & 46,44,45, 43,41,42, 38,39,37,38,39,40, 32,32,32,36,30,31,32,32,32,33,34,35, 54,56,53,54,55,53, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16,17,17,17,17,17,17, 24,24,24,24,24,24, & 44,46,45, 41,43,42, 39,38,37,39,38,40, 33,32,32,32,31,30,36,32,32,32,35,34, 54,53,56,54,53,55, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,16,17,17,17,17,17, 24,24,24,24,24,24, & ! 1. pyr<c+a> 44,45,46, 41,42,43, 40,38,39,37,38,39, 34,35,32,32,32,36,30,31,32,32,32,33, 53,54,56,53,54,55, & ! 1. pyr<c+a>
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,16,17,17,17,17, 24,24,24,24,24,24, & 46,45,44, 43,42,41, 40,39,38,37,39,38, 35,34,33,32,32,32,31,30,36,32,32,32, 55,54,53,56,54,53, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,16,17,17,17, 24,24,24,24,24,24, & 45,46,44, 42,43,41, 39,40,38,39,37,38, 32,33,34,35,32,32,32,36,30,31,32,32, 55,53,54,56,53,54, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,17,16,17,17, 24,24,24,24,24,24, & 45,44,46, 42,41,43, 38,40,39,38,37,39, 32,32,35,34,33,32,32,32,31,30,36,32, 53,55,54,53,56,54, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,17,17,16,17, 24,24,24,24,24,24, & 46,44,45, 43,41,42, 38,39,40,38,39,37, 32,32,32,33,34,35,32,32,32,36,30,31, 54,55,53,54,56,53, &
20,20,20, 19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,17,17,17,17,17,17,16, 24,24,24,24,24,24, & 44,46,45, 41,43,42, 39,38,40,39,38,37, 36,32,32,32,35,34,33,32,32,32,31,30, 54,53,55,54,53,56, &
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 25,26,26,26,26,26, & 68,68,69, 66,66,67, 64,64,65,64,65,65, 60,61,61,60,62,62,60,63,63,60,62,62, 57,58,58,59,58,58, &
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,25,26,26,26,26, & 68,69,68, 66,67,66, 65,64,64,65,64,64, 62,62,60,61,61,60,62,62,60,63,63,60, 58,57,58,58,59,58, &
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,25,26,26,26, & 69,68,68, 67,66,66, 64,65,64,64,65,64, 63,60,62,62,60,61,61,60,62,62,60,63, 58,58,57,58,58,59, &
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,26,25,26,26, & ! 2. pyr<c+a> 68,68,69, 66,66,67, 64,64,65,64,64,65, 60,63,63,60,62,62,60,61,61,60,62,62, 59,58,58,57,58,58, & ! 2. pyr<c+a>
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,26,26,25,26, & 68,69,68, 66,67,66, 65,64,64,65,64,64, 62,62,60,63,63,60,62,62,60,61,61,60, 58,59,58,58,57,58, &
30,30,30, 29,29,29, 28,28,28,28,28,28, 27,27,27,27,27,27,27,27,27,27,27,27, 26,26,26,26,26,25 & 69,68,68, 67,66,66, 64,65,64,64,65,64, 61,60,62,62,60,63,63,60,62,62,60,61, 58,58,59,58,58,57 &
],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for hex (onion peel naming scheme) ],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for hex (onion peel naming scheme)
!< 10.1016/j.ijplas.2014.06.010 table 3
!< 10.1080/14786435.2012.699689 table 2 and 3
!< index & label & description
!< 1 & S1 & basal self-interaction
!< 2 & 1 & basal/basal coplanar
!< 3 & 3 & basal/prismatic collinear
!< 4 & 4 & basal/prismatic non-collinear
!< 5 & S2 & prismatic self-interaction
!< 6 & 2 & prismatic/prismatic
!< 7 & 5 & prismatic/basal collinear
!< 8 & 6 & prismatic/basal non-collinear
!< 9 & - & basal/pyramidal <a> non-collinear
!< 10 & - & basal/pyramidal <a> collinear
!< 11 & - & prismatic/pyramidal <a> non-collinear
!< 12 & - & prismatic/pyramidal <a> collinear
!< 13 & - & pyramidal <a> self-interaction
!< 14 & - & pyramidal <a> non-collinear
!< 15 & - & pyramidal <a> collinear
!< 16 & - & pyramidal <a>/prismatic non-collinear
!< 17 & - & pyramidal <a>/prismatic collinear
!< 18 & - & pyramidal <a>/basal non-collinear
!< 19 & - & pyramidal <a>/basal collinear
!< 20 & - & basal/1. order pyramidal <c+a> semi-collinear
!< 21 & - & basal/1. order pyramidal <c+a>
!< 22 & - & basal/1. order pyramidal <c+a>
!< 23 & - & prismatic/1. order pyramidal <c+a> semi-collinear
!< 24 & - & prismatic/1. order pyramidal <c+a>
!< 25 & - & prismatic/1. order pyramidal <c+a> semi-coplanar?
!< 26 & - & pyramidal <a>/1. order pyramidal <c+a> coplanar
!< 27 & - & pyramidal <a>/1. order pyramidal <c+a>
!< 28 & - & pyramidal <a>/1. order pyramidal <c+a> semi-collinear
!< 29 & - & pyramidal <a>/1. order pyramidal <c+a> semi-coplanar
!< 30 & - & 1. order pyramidal <c+a> self-interaction
!< 31 & - & 1. order pyramidal <c+a> coplanar
!< 32 & - & 1. order pyramidal <c+a>
!< 33 & - & 1. order pyramidal <c+a>
!< 34 & - & 1. order pyramidal <c+a> semi-coplanar
!< 35 & - & 1. order pyramidal <c+a> semi-coplanar
!< 36 & - & 1. order pyramidal <c+a> collinear
!< 37 & - & 1. order pyramidal <c+a>/pyramidal <a> coplanar
!< 38 & - & 1. order pyramidal <c+a>/pyramidal <a> semi-collinear
!< 39 & - & 1. order pyramidal <c+a>/pyramidal <a>
!< 40 & - & 1. order pyramidal <c+a>/pyramidal <a> semi-coplanar
!< 41 & - & 1. order pyramidal <c+a>/prismatic semi-collinear
!< 42 & - & 1. order pyramidal <c+a>/prismatic semi-coplanar
!< 43 & - & 1. order pyramidal <c+a>/prismatic
!< 44 & - & 1. order pyramidal <c+a>/basal semi-collinear
!< 45 & - & 1. order pyramidal <c+a>/basal
!< 46 & - & 1. order pyramidal <c+a>/basal
!< 47 & 8 & basal/2. order pyramidal <c+a> non-collinear
!< 48 & 7 & basal/2. order pyramidal <c+a> semi-collinear
!< 49 & 10 & prismatic/2. order pyramidal <c+a>
!< 50 & 9 & prismatic/2. order pyramidal <c+a> semi-collinear
!< 51 & - & pyramidal <a>/2. order pyramidal <c+a>
!< 52 & - & pyramidal <a>/2. order pyramidal <c+a> semi collinear
!< 53 & - & 1. order pyramidal <c+a>/2. order pyramidal <c+a>
!< 54 & - & 1. order pyramidal <c+a>/2. order pyramidal <c+a>
!< 55 & - & 1. order pyramidal <c+a>/2. order pyramidal <c+a>
!< 56 & - & 1. order pyramidal <c+a>/2. order pyramidal <c+a> collinear
!< 57 & S3 & 2. order pyramidal <c+a> self-interaction
!< 58 & 16 & 2. order pyramidal <c+a> non-collinear
!< 59 & 15 & 2. order pyramidal <c+a> semi-collinear
!< 60 & - & 2. order pyramidal <c+a>/1. order pyramidal <c+a>
!< 61 & - & 2. order pyramidal <c+a>/1. order pyramidal <c+a> collinear
!< 62 & - & 2. order pyramidal <c+a>/1. order pyramidal <c+a>
!< 63 & - & 2. order pyramidal <c+a>/1. order pyramidal <c+a>
!< 64 & - & 2. order pyramidal <c+a>/pyramidal <a> non-collinear
!< 65 & - & 2. order pyramidal <c+a>/pyramidal <a> semi-collinear
!< 66 & 14 & 2. order pyramidal <c+a>/prismatic non-collinear
!< 67 & 13 & 2. order pyramidal <c+a>/prismatic semi-collinear
!< 68 & 12 & 2. order pyramidal <c+a>/basal non-collinear
!< 69 & 11 & 2. order pyramidal <c+a>/basal semi-collinear
integer, dimension(BCT_NSLIP,BCT_NSLIP), parameter :: & integer, dimension(BCT_NSLIP,BCT_NSLIP), parameter :: &
BCT_INTERACTIONSLIPSLIP = reshape( [& BCT_INTERACTIONSLIPSLIP = reshape( [&

11
src/material.f90
View File

@ -108,8 +108,14 @@ subroutine parse()
homogenizations => config_material%get('homogenization') homogenizations => config_material%get('homogenization')
call sanityCheck(materials, homogenizations) call sanityCheck(materials, homogenizations)
#if defined (__GFORTRAN__)
material_name_phase = getKeys(phases) material_name_phase = getKeys(phases)
material_name_homogenization = getKeys(homogenizations) material_name_homogenization = getKeys(homogenizations)
#else
material_name_phase = phases%Keys()
material_name_homogenization = homogenizations%Keys()
#endif
allocate(homogenization_Nconstituents(homogenizations%length)) allocate(homogenization_Nconstituents(homogenizations%length))
do h=1, homogenizations%length do h=1, homogenizations%length
@ -203,9 +209,9 @@ subroutine sanityCheck(materials,homogenizations)
end subroutine sanityCheck end subroutine sanityCheck
#if defined (__GFORTRAN__)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Get all keys from a dictionary !> @brief %keys() is broken on gfortran
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function getKeys(dict) function getKeys(dict)
@ -228,5 +234,6 @@ function getKeys(dict)
end do end do
end function getKeys end function getKeys
#endif
end module material end module material

30
src/math.f90
View File

@ -84,38 +84,34 @@ contains
subroutine math_init subroutine math_init
real(pReal), dimension(4) :: randTest real(pReal), dimension(4) :: randTest
integer :: & integer :: randSize
randSize, & integer, dimension(:), allocatable :: seed
randomSeed !< fixed seeding for pseudo-random number generator, Default 0: use random seed
integer, dimension(:), allocatable :: randInit
class(tNode), pointer :: & class(tNode), pointer :: &
num_generic num_generic
print'(/,1x,a)', '<<<+- math init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- math init -+>>>'; flush(IO_STDOUT)
num_generic => config_numerics%get('generic',defaultVal=emptyDict) num_generic => config_numerics%get('generic',defaultVal=emptyDict)
randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0)
call random_seed(size=randSize) call random_seed(size=randSize)
allocate(randInit(randSize)) allocate(seed(randSize))
if (randomSeed > 0) then
randInit = randomSeed if (num_generic%contains('random_seed')) then
seed = num_generic%get_as1dInt('random_seed',requiredSize=randSize)
else else
call random_seed() call random_seed()
call random_seed(get = randInit) call random_seed(get = seed)
randInit(2:randSize) = randInit(1) end if
endif
call random_seed(put = randInit) call random_seed(put = seed)
call random_number(randTest) call random_number(randTest)
print'(/,a,i2)', ' size of random seed: ', randSize print'(/,a,i2)', ' size of random seed: ', randSize
print'( a,i0)', ' value of random seed: ', randInit(1) print*, 'value of random seed: ', seed
print'( a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest print'( a,4(/,26x,f17.14))', ' start of random sequence: ', randTest
call random_seed(put = randInit) call selfTest()
call selfTest
end subroutine math_init end subroutine math_init

2
src/parallelization.f90
View File

@ -53,6 +53,7 @@ contains
subroutine parallelization_init subroutine parallelization_init
integer(MPI_INTEGER_KIND) :: err_MPI, typeSize integer(MPI_INTEGER_KIND) :: err_MPI, typeSize
character(len=4) :: rank_str
!$ integer :: got_env, threadLevel !$ integer :: got_env, threadLevel
!$ integer(pI32) :: OMP_NUM_THREADS !$ integer(pI32) :: OMP_NUM_THREADS
!$ character(len=6) NumThreadsString !$ character(len=6) NumThreadsString
@ -112,6 +113,7 @@ subroutine parallelization_init
if (worldrank /= 0) then if (worldrank /= 0) then
close(OUTPUT_UNIT) ! disable output close(OUTPUT_UNIT) ! disable output
write(rank_str,'(i4.4)') worldrank ! use for MPI debug filenames
open(OUTPUT_UNIT,file='/dev/null',status='replace') ! close() alone will leave some temp files in cwd open(OUTPUT_UNIT,file='/dev/null',status='replace') ! close() alone will leave some temp files in cwd
endif endif

4
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -99,11 +99,11 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
associate(prm => param(kinematics_thermal_expansion_instance(ph))) associate(prm => param(kinematics_thermal_expansion_instance(ph)))
Li = dot_T * ( & Li = dot_T * ( &
prm%A(1:3,1:3,1) & ! constant coefficient prm%A(1:3,1:3,1) & ! constant coefficient
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient + prm%A(1:3,1:3,2)*(T - prm%T_ref) & ! linear coefficient
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
) / & ) / &
(1.0_pReal & (1.0_pReal &
+ prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1.0_pReal & + prm%A(1:3,1:3,1)*(T - prm%T_ref) / 1.0_pReal &
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2.0_pReal & + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2.0_pReal &
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3.0_pReal & + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3.0_pReal &
) )

12
src/phase_mechanical_elastic.f90
View File

@ -101,32 +101,32 @@ pure module function elastic_C66(ph,en) result(C66)
T = thermal_T(ph,en) T = thermal_T(ph,en)
C66(1,1) = prm%C_11(1) & C66(1,1) = prm%C_11(1) &
+ prm%C_11(2)*(T - prm%T_ref)**1 & + prm%C_11(2)*(T - prm%T_ref) &
+ prm%C_11(3)*(T - prm%T_ref)**2 + prm%C_11(3)*(T - prm%T_ref)**2
C66(1,2) = prm%C_12(1) & C66(1,2) = prm%C_12(1) &
+ prm%C_12(2)*(T - prm%T_ref)**1 & + prm%C_12(2)*(T - prm%T_ref) &
+ prm%C_12(3)*(T - prm%T_ref)**2 + prm%C_12(3)*(T - prm%T_ref)**2
C66(4,4) = prm%C_44(1) & C66(4,4) = prm%C_44(1) &
+ prm%C_44(2)*(T - prm%T_ref)**1 & + prm%C_44(2)*(T - prm%T_ref) &
+ prm%C_44(3)*(T - prm%T_ref)**2 + prm%C_44(3)*(T - prm%T_ref)**2
if (any(phase_lattice(ph) == ['hP','tI'])) then if (any(phase_lattice(ph) == ['hP','tI'])) then
C66(1,3) = prm%C_13(1) & C66(1,3) = prm%C_13(1) &
+ prm%C_13(2)*(T - prm%T_ref)**1 & + prm%C_13(2)*(T - prm%T_ref) &
+ prm%C_13(3)*(T - prm%T_ref)**2 + prm%C_13(3)*(T - prm%T_ref)**2
C66(3,3) = prm%C_33(1) & C66(3,3) = prm%C_33(1) &
+ prm%C_33(2)*(T - prm%T_ref)**1 & + prm%C_33(2)*(T - prm%T_ref) &
+ prm%C_33(3)*(T - prm%T_ref)**2 + prm%C_33(3)*(T - prm%T_ref)**2
end if end if
if (phase_lattice(ph) == 'tI') then if (phase_lattice(ph) == 'tI') then
C66(6,6) = prm%C_66(1) & C66(6,6) = prm%C_66(1) &
+ prm%C_66(2)*(T - prm%T_ref)**1 & + prm%C_66(2)*(T - prm%T_ref) &
+ prm%C_66(3)*(T - prm%T_ref)**2 + prm%C_66(3)*(T - prm%T_ref)**2
end if end if

14
src/rotations.f90
View File

@ -270,7 +270,7 @@ pure elemental subroutine standardize(self)
class(rotation), intent(inout) :: self class(rotation), intent(inout) :: self
if (self%q(1) < 0.0_pReal) self%q = - self%q if (sign(1.0_pReal,self%q(1)) < 0.0_pReal) self%q = - self%q
end subroutine standardize end subroutine standardize
@ -450,7 +450,7 @@ pure function qu2om(qu) result(om)
om(3,2) = 2.0_pReal*(qu(4)*qu(3)+qu(1)*qu(2)) om(3,2) = 2.0_pReal*(qu(4)*qu(3)+qu(1)*qu(2))
om(1,3) = 2.0_pReal*(qu(2)*qu(4)+qu(1)*qu(3)) om(1,3) = 2.0_pReal*(qu(2)*qu(4)+qu(1)*qu(3))
if (P < 0.0_pReal) om = transpose(om) if (sign(1.0_pReal,P) < 0.0_pReal) om = transpose(om)
end function qu2om end function qu2om
@ -480,7 +480,7 @@ pure function qu2eu(qu) result(eu)
atan2( 2.0_pReal*chi, q03-q12 ), & atan2( 2.0_pReal*chi, q03-q12 ), &
atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )] atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
endif degenerated endif degenerated
where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
end function qu2eu end function qu2eu
@ -602,7 +602,7 @@ pure function om2qu(om) result(qu)
qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s] qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s]
endif endif
endif endif
if(qu(1)<0._pReal) qu =-1.0_pReal * qu if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
qu = qu*[1.0_pReal,P,P,P] qu = qu*[1.0_pReal,P,P,P]
end function om2qu end function om2qu
@ -628,7 +628,7 @@ pure function om2eu(om) result(eu)
eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ] eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ]
end if end if
where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal
where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
end function om2eu end function om2eu
@ -735,7 +735,7 @@ pure function eu2qu(eu) result(qu)
-P*sPhi*cos(ee(1)-ee(3)), & -P*sPhi*cos(ee(1)-ee(3)), &
-P*sPhi*sin(ee(1)-ee(3)), & -P*sPhi*sin(ee(1)-ee(3)), &
-P*cPhi*sin(ee(1)+ee(3))] -P*cPhi*sin(ee(1)+ee(3))]
if(qu(1) < 0.0_pReal) qu = qu * (-1.0_pReal) if(sign(1.0_pReal,qu(1)) < 0.0_pReal) qu = qu * (-1.0_pReal)
end function eu2qu end function eu2qu
@ -792,7 +792,7 @@ pure function eu2ax(eu) result(ax)
else else
ax(1:3) = -P/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ] ! passive axis-angle pair so a minus sign in front ax(1:3) = -P/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ] ! passive axis-angle pair so a minus sign in front
ax(4) = alpha ax(4) = alpha
if (alpha < 0.0_pReal) ax = -ax ! ensure alpha is positive if (sign(1.0_pReal,alpha) < 0.0_pReal) ax = -ax ! ensure alpha is positive
end if end if
end function eu2ax end function eu2ax