diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..293b66b --- /dev/null +++ b/.clang-format @@ -0,0 +1,134 @@ +# SPDX-FileCopyrightText: 2021 SeisSol Group +# +# SPDX-License-Identifier: BSD-3-Clause +# SPDX-LicenseComments: Full text under /LICENSE and /LICENSES/ +# +# SPDX-FileContributor: Author lists in /AUTHORS and /CITATION.cff + +Language: Cpp +AccessModifierOffset: -4 +AlignAfterOpenBracket: Align +AlignConsecutiveMacros: false +AlignConsecutiveAssignments: false +AlignConsecutiveDeclarations: false +AlignEscapedNewlines: Right +AlignOperands: true +AlignTrailingComments: true +AllowAllArgumentsOnNextLine: true +AllowAllConstructorInitializersOnNextLine: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortBlocksOnASingleLine: false +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortLambdasOnASingleLine: All +AllowShortIfStatementsOnASingleLine: Never +AllowShortLoopsOnASingleLine: false +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: Yes +BinPackArguments: false +BinPackParameters: false +BraceWrapping: + AfterCaseLabel: false + AfterClass: false + AfterControlStatement: false + AfterEnum: false + AfterFunction: false + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + AfterExternBlock: false + BeforeCatch: false + BeforeElse: false + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeBraces: Attach +BreakBeforeInheritanceComma: false +BreakInheritanceList: BeforeColon +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +BreakConstructorInitializers: BeforeColon +BreakAfterJavaFieldAnnotations: false +BreakStringLiterals: true +ColumnLimit: 100 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerAllOnOneLineOrOnePerLine: false +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DerivePointerAlignment: false +DisableFormat: false +ExperimentalAutoDetectBinPacking: false +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IncludeBlocks: Regroup +IncludeCategories: + # keep the doctest headers in front + - Regex: '^(<|")doctest' + Priority: 1 + - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + Priority: 3 + - Regex: '^(<|"(gtest|gmock|isl|json)/)' + Priority: 4 + - Regex: '.*' + Priority: 2 +IncludeIsMainRegex: '(Test)?$' +IndentCaseLabels: false +IndentPPDirectives: None +IndentWidth: 2 +IndentWrappedFunctionNames: true +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: true +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBinPackProtocolList: Auto +ObjCBlockIndentWidth: 2 +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyBreakTemplateDeclaration: 10 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 60 +PointerAlignment: Left +QualifierAlignment: Left +ReflowComments: true +SortIncludes: true +SortUsingDeclarations: true +SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: true +SpaceBeforeAssignmentOperators: true +SpaceBeforeCpp11BracedList: false +SpaceBeforeCtorInitializerColon: true +SpaceBeforeInheritanceColon: true +SpaceBeforeParens: ControlStatements +SpaceBeforeRangeBasedForLoopColon: true +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInParentheses: false +SpacesInSquareBrackets: false +Standard: Cpp11 +StatementMacros: + - Q_UNUSED + - QT_REQUIRE_VERSION +TabWidth: 8 +UseTab: Never diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs new file mode 100644 index 0000000..c2ae4d6 --- /dev/null +++ b/.git-blame-ignore-revs @@ -0,0 +1,11 @@ +# SPDX-FileCopyrightText: 2026 SeisSol Group +# +# SPDX-License-Identifier: BSD-3-Clause +# SPDX-LicenseComments: Full text under /LICENSE and /LICENSES/ +# +# SPDX-FileContributor: Author lists in /AUTHORS and /CITATION.cff + +# apply pre-commit (whitespaces etc.) +412e0d94440dd4371f6d269f16ef6b50396e2907 +# apply clang-format +c007c326b54aa0f43deea770ba17ad465eb3d107 diff --git a/.github/workflows/pre-commit.yml b/.github/workflows/pre-commit.yml new file mode 100644 index 0000000..6f56a26 --- /dev/null +++ b/.github/workflows/pre-commit.yml @@ -0,0 +1,21 @@ +# SPDX-FileCopyrightText: 2025 SeisSol Group +# +# SPDX-License-Identifier: BSD-3-Clause +# SPDX-LicenseComments: Full text under /LICENSE and /LICENSES/ +# +# SPDX-FileContributor: Author lists in /AUTHORS and /CITATION.cff + +name: pre-commit +on: + - push + +jobs: + pre-commit: + name: pre-commit + runs-on: ubuntu-24.04 + steps: + - uses: actions/checkout@v6 + with: + submodules: recursive + + - uses: pre-commit/action@v3.0.1 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..fc7c532 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,92 @@ +# SPDX-FileCopyrightText: 2025 SeisSol Group +# +# SPDX-License-Identifier: BSD-3-Clause +# SPDX-LicenseComments: Full text under /LICENSE and /LICENSES/ +# +# SPDX-FileContributor: Author lists in /AUTHORS and /CITATION.cff + +--- + +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v6.0.0 + hooks: + - id: check-merge-conflict + name: '[GENERIC] merge conflict check' + - id: check-symlinks + name: '[GENERIC] symlink check' + - id: destroyed-symlinks + name: '[GENERIC] detect broken symlinks' + - id: detect-private-key + name: '[GENERIC] detect private keys uploaded by accident' + - id: check-case-conflict + name: '[GENERIC] detect OS file naming case conflicts' + - id: check-executables-have-shebangs + name: '[GENERIC] check for shebangs in executable files' + - id: check-illegal-windows-names + name: '[GENERIC] detect illegal Windows file names' + - id: check-json + name: '[JSON] check' + - id: check-xml + name: '[XML] check' + - id: check-shebang-scripts-are-executable + name: '[GENERIC] check that shebang-containing files are executable' + + +- repo: https://github.com/DavidAnson/markdownlint-cli2 + rev: v0.18.1 + hooks: + - id: markdownlint-cli2 + name: '[MARKDOWN] lint' + +#- repo: https://github.com/fsfe/reuse-tool +# rev: v6.0.0 +# hooks: +# - id: reuse +# name: '[GENERIC] REUSE compatibiltiy' + +#- repo: https://github.com/psf/black-pre-commit-mirror +# rev: 25.1.0 +# hooks: +# - id: black +# files: ^(?!preprocessing|postprocessing) +# name: '[PYTHON] black' +#- repo: https://github.com/pycqa/isort +# rev: 6.0.1 +# hooks: +# - id: isort +# files: ^(?!preprocessing|postprocessing) +# args: ["--profile", "black"] +# name: '[PYTHON] isort' +- repo: https://github.com/pycqa/bandit + rev: 1.8.6 + hooks: + - id: bandit + args: ["--confidence-level", "high", "--severity-level", "high"] + name: '[PYTHON] bandit' +#- repo: https://github.com/pycqa/flake8 +# rev: '7.3.0' +# hooks: +# - id: flake8 +# files: ^(?!preprocessing|postprocessing) +# name: '[PYTHON] Flake8' + +- repo: https://github.com/sphinx-contrib/sphinx-lint + rev: 'v1.0.0' + hooks: + - id: sphinx-lint + name: '[SPHINX/RST] sphinx lint' + +- repo: https://github.com/pre-commit/mirrors-clang-format + rev: 'v22.1.0' + hooks: + - id: clang-format + name: '[C++] clang-format' + +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v6.0.0 + hooks: + - id: end-of-file-fixer + name: '[GENERIC] newline eof' + - id: trailing-whitespace + name: '[GENERIC] remove trailing whitespace' diff --git a/README.md b/README.md index 8f2656b..8ff0fb9 100644 --- a/README.md +++ b/README.md @@ -1 +1,3 @@ +# Preprocessing Tools for SeisSol + A collection of useful tools for mesh generation and preprocessing for [SeisSol](http://www.seissol.org). diff --git a/SimModelerDownloadingBuilding/README.md b/SimModelerDownloadingBuilding/README.md index 50bc965..bc2f74d 100644 --- a/SimModelerDownloadingBuilding/README.md +++ b/SimModelerDownloadingBuilding/README.md @@ -1,20 +1,31 @@ # Downloading SimModeler -Go to [the support area of Simmetrix web site](http://www.simmetrix.com/index.php/support/support-downloads) and download the binaries of the code. +Go to [the support area of Simmetrix web site](http://www.simmetrix.com/index.php/support/support-downloads) +and download the binaries of the code. -# Customizing SimModeler for SeisSol -Add the files of [this folder](https://github.com/SeisSol/Meshing/tree/master/SimModelerDownloadingBuilding/SimModelerCustomization) to the main directory of SimModeler. These files allow defining properly boundary conditions and exporting the mesh in the proper format. +## Customizing SimModeler for SeisSol +Add the files of [this folder](https://github.com/SeisSol/Meshing/tree/master/SimModelerDownloadingBuilding/SimModelerCustomization) +to the main directory of SimModeler. These files allow +defining properly boundary conditions and exporting the +mesh in the proper format. -# Downloading SimModeler modeling suite +## Downloading SimModeler modeling suite -Go to [the support area of Simmetrix web site](http://www.simmetrix.com/index.php/support/support-downloads) to see the available versions of SimModeler modeling suite. Then use this [script](https://github.com/SeisSol/Meshing/tree/master/SimModelerDownloadingBuilding/downloadSimLib.sh) to download all components of the library. Please make sure that all components that you download are included in your licence agreeement. +Go to [the support area of Simmetrix web site](http://www.simmetrix.com/index.php/support/support-downloads) +to see the available versions of SimModeler modeling suite. +Then use this [script](https://github.com/SeisSol/Meshing/tree/master/SimModelerDownloadingBuilding/downloadSimLib.sh) +to download all components of the library. Please +make sure that all components that you download +are included in your licence agreeement. Finally, untar all the files, using for instance: + ```bash for filename in *.tgz do tar zxf $filename done ``` + The libraries are precompiled, so no need to build them! diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/README.md b/SimModelerDownloadingBuilding/SimModelerCustomization/README.md index 780544b..8940b3e 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/README.md +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/README.md @@ -1 +1,6 @@ -These files need to be added to the main directory of SimModeler to be able to properly define the boundary conditions and to be able to export the mesh in the proper format. +# Note + +These files need to be added to the main directory +of SimModeler to be able to properly define the +boundary conditions and to be able to export +the mesh in the proper format. diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/seisSol.adi b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/seisSol.adi index 1854a38..051a8cf 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/seisSol.adi +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/seisSol.adi @@ -4,4 +4,3 @@ SolverName=seisSol [0] AttDefs=simSeisSol/attDefs/v000/ PatternInfo=simSeisSol/export/case/v000/gambit_seisSol.sxi - diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/66/6633a98134d729503662fdb7d18fb7fceeb22283.svn-base b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/66/6633a98134d729503662fdb7d18fb7fceeb22283.svn-base index 709b35a..10208e1 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/66/6633a98134d729503662fdb7d18fb7fceeb22283.svn-base +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/66/6633a98134d729503662fdb7d18fb7fceeb22283.svn-base @@ -4,4 +4,3 @@ SolverName=seisSol [0] AttDefs=simSeisSol/attDefs/v000/ PatternInfo=case/simSeisSol/export/case/v000/gambit_seisSol.sxi - diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/6d/6d97471d1a600f14e76448d9c0150455786aa03d.svn-base b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/6d/6d97471d1a600f14e76448d9c0150455786aa03d.svn-base index 4afa7aa..8f58704 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/6d/6d97471d1a600f14e76448d9c0150455786aa03d.svn-base +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/6d/6d97471d1a600f14e76448d9c0150455786aa03d.svn-base @@ -17,11 +17,11 @@ define-procedure hasBoundaryCondition(attName) { gfTag = $tag attributes["boundaryCondition"]/* { match = function:strcmp($image,$(inputAttName)) - hasBcAtt = function:notI($match) + hasBcAtt = function:notI($match) } } - hasBoundaryCondition = $hasBcAtt + hasBoundaryCondition = $hasBcAtt } define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version) { @@ -55,11 +55,11 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version numFaces = function:addI($count,$numFaces) info:header = " num classified faces on gFace $gfTag = $count ... total classified = $numFaces" } - #get number of regions with a face classified on gface + #get number of regions with a face classified on gface mfaces/*/mregions/* { regionExists = function:mapI2I($id,"<-1 0>", 1) skip = function:notI($regionExists) - numRegions = function:addI(1,$numRegions) + numRegions = function:addI(1,$numRegions) } info:header = " num regions with a face classied on gFace $gfTag = $numRegions" } @@ -68,11 +68,11 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version info:header = " writing $numFaces faces .... " info:header = " writing $numRegions regions .... " - #write the header string + #write the header string header = "BOUNDARY CONDITIONS $version \n" + " $bcID 1 $numRegions 0 6" - #loop over the mesh faces classified on the geometric model faces with the bcAttribute + #loop over the mesh faces classified on the geometric model faces with the bcAttribute # and write the mesh elements and the face number gmodel/gfaces/* { gfTag = $tag @@ -90,38 +90,38 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version } skip = function:notI($hasBcAtt) - #need to determine the vertex numbering so the face numbering can be determined - # see www.stanford.edu/class/me469b/handouts/gambit_write.pdf pg C-33 - # for the tet ordering + #need to determine the vertex numbering so the face numbering can be determined + # see www.stanford.edu/class/me469b/handouts/gambit_write.pdf pg C-33 + # for the tet ordering # gambit and MeshSim have the same ordering for tets # see the R_face documentation for the MeshSim ordering mfaces/* { mRgnId_A = $(mregions[0]/id) mRgnId_B = $(mregions[1]/id) - . { + . { regionExists = function:mapI2I($mRgnId_A,"<-1 0>", 1) skip = function:notI($regionExists) # Need the raw mesh region id to call mfaceNum() fnum = function:mfaceNum($mRgnId_A) snum = function:addI($fnum, 1) # for each face write: - # - # <6 for tet> + # + # <6 for tet> # - header = " $mRgnId_A 6 $snum" + header = " $mRgnId_A 6 $snum" } # end region A - . { + . { regionExists = function:mapI2I($mRgnId_B,"<-1 0>", 1) skip = function:notI($regionExists) # Need the raw mesh region id to call mfaceNum() fnum = function:mfaceNum($mRgnId_B) snum = function:addI($fnum, 1) # for each face write: - # - # <6 for tet> + # + # <6 for tet> # - header = " $mRgnId_B 6 $snum" + header = " $mRgnId_B 6 $snum" } # end region B } # end faces } # end gfaces @@ -244,7 +244,7 @@ mesh/mregions { header = "$(id,I8) $(ntype,I2) $(ndp,I2) $(mvertices[*]/id,I8)" } . { - # Pyrmaid elements + # Pyrmaid elements skip = function:compI($ndp,5) header = "$(id,I8) $(ntype,I2) $(ndp,I2) " + "$(mvertices[0]/id,I8)$(mvertices[1]/id,I8)$(mvertices[3]/id,I8)$(mvertices[2]/id,I8)" + diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/e8/e8a0186ae71e0dac8f12bb52f2669169e15522d8.svn-base b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/e8/e8a0186ae71e0dac8f12bb52f2669169e15522d8.svn-base index b5cae27..9b454e5 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/e8/e8a0186ae71e0dac8f12bb52f2669169e15522d8.svn-base +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/.svn/pristine/e8/e8a0186ae71e0dac8f12bb52f2669169e15522d8.svn-base @@ -3,4 +3,3 @@ atd 13 image "seisSol" : "analysis" root { C "problem definition":"seisSol" 0 1 1; }; - diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/attDefs/v000/analysis.atd b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/attDefs/v000/analysis.atd index b5cae27..9b454e5 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/attDefs/v000/analysis.atd +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/attDefs/v000/analysis.atd @@ -3,4 +3,3 @@ atd 13 image "seisSol" : "analysis" root { C "problem definition":"seisSol" 0 1 1; }; - diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/export/case/v000/gambit_seisSol.sxp b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/export/case/v000/gambit_seisSol.sxp index 4afa7aa..8f58704 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/export/case/v000/gambit_seisSol.sxp +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/export/case/v000/gambit_seisSol.sxp @@ -17,11 +17,11 @@ define-procedure hasBoundaryCondition(attName) { gfTag = $tag attributes["boundaryCondition"]/* { match = function:strcmp($image,$(inputAttName)) - hasBcAtt = function:notI($match) + hasBcAtt = function:notI($match) } } - hasBoundaryCondition = $hasBcAtt + hasBoundaryCondition = $hasBcAtt } define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version) { @@ -55,11 +55,11 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version numFaces = function:addI($count,$numFaces) info:header = " num classified faces on gFace $gfTag = $count ... total classified = $numFaces" } - #get number of regions with a face classified on gface + #get number of regions with a face classified on gface mfaces/*/mregions/* { regionExists = function:mapI2I($id,"<-1 0>", 1) skip = function:notI($regionExists) - numRegions = function:addI(1,$numRegions) + numRegions = function:addI(1,$numRegions) } info:header = " num regions with a face classied on gFace $gfTag = $numRegions" } @@ -68,11 +68,11 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version info:header = " writing $numFaces faces .... " info:header = " writing $numRegions regions .... " - #write the header string + #write the header string header = "BOUNDARY CONDITIONS $version \n" + " $bcID 1 $numRegions 0 6" - #loop over the mesh faces classified on the geometric model faces with the bcAttribute + #loop over the mesh faces classified on the geometric model faces with the bcAttribute # and write the mesh elements and the face number gmodel/gfaces/* { gfTag = $tag @@ -90,38 +90,38 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version } skip = function:notI($hasBcAtt) - #need to determine the vertex numbering so the face numbering can be determined - # see www.stanford.edu/class/me469b/handouts/gambit_write.pdf pg C-33 - # for the tet ordering + #need to determine the vertex numbering so the face numbering can be determined + # see www.stanford.edu/class/me469b/handouts/gambit_write.pdf pg C-33 + # for the tet ordering # gambit and MeshSim have the same ordering for tets # see the R_face documentation for the MeshSim ordering mfaces/* { mRgnId_A = $(mregions[0]/id) mRgnId_B = $(mregions[1]/id) - . { + . { regionExists = function:mapI2I($mRgnId_A,"<-1 0>", 1) skip = function:notI($regionExists) # Need the raw mesh region id to call mfaceNum() fnum = function:mfaceNum($mRgnId_A) snum = function:addI($fnum, 1) # for each face write: - # - # <6 for tet> + # + # <6 for tet> # - header = " $mRgnId_A 6 $snum" + header = " $mRgnId_A 6 $snum" } # end region A - . { + . { regionExists = function:mapI2I($mRgnId_B,"<-1 0>", 1) skip = function:notI($regionExists) # Need the raw mesh region id to call mfaceNum() fnum = function:mfaceNum($mRgnId_B) snum = function:addI($fnum, 1) # for each face write: - # - # <6 for tet> + # + # <6 for tet> # - header = " $mRgnId_B 6 $snum" + header = " $mRgnId_B 6 $snum" } # end region B } # end faces } # end gfaces @@ -244,7 +244,7 @@ mesh/mregions { header = "$(id,I8) $(ntype,I2) $(ndp,I2) $(mvertices[*]/id,I8)" } . { - # Pyrmaid elements + # Pyrmaid elements skip = function:compI($ndp,5) header = "$(id,I8) $(ntype,I2) $(ndp,I2) " + "$(mvertices[0]/id,I8)$(mvertices[1]/id,I8)$(mvertices[3]/id,I8)$(mvertices[2]/id,I8)" + diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/seisSol.adi b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/seisSol.adi index 709b35a..10208e1 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/seisSol.adi +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/attDefs/simSeisSol/seisSol.adi @@ -4,4 +4,3 @@ SolverName=seisSol [0] AttDefs=simSeisSol/attDefs/v000/ PatternInfo=case/simSeisSol/export/case/v000/gambit_seisSol.sxi - diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/66/6633a98134d729503662fdb7d18fb7fceeb22283.svn-base b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/66/6633a98134d729503662fdb7d18fb7fceeb22283.svn-base index 709b35a..10208e1 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/66/6633a98134d729503662fdb7d18fb7fceeb22283.svn-base +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/66/6633a98134d729503662fdb7d18fb7fceeb22283.svn-base @@ -4,4 +4,3 @@ SolverName=seisSol [0] AttDefs=simSeisSol/attDefs/v000/ PatternInfo=case/simSeisSol/export/case/v000/gambit_seisSol.sxi - diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/6d/6d97471d1a600f14e76448d9c0150455786aa03d.svn-base b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/6d/6d97471d1a600f14e76448d9c0150455786aa03d.svn-base index 4afa7aa..8f58704 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/6d/6d97471d1a600f14e76448d9c0150455786aa03d.svn-base +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/6d/6d97471d1a600f14e76448d9c0150455786aa03d.svn-base @@ -17,11 +17,11 @@ define-procedure hasBoundaryCondition(attName) { gfTag = $tag attributes["boundaryCondition"]/* { match = function:strcmp($image,$(inputAttName)) - hasBcAtt = function:notI($match) + hasBcAtt = function:notI($match) } } - hasBoundaryCondition = $hasBcAtt + hasBoundaryCondition = $hasBcAtt } define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version) { @@ -55,11 +55,11 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version numFaces = function:addI($count,$numFaces) info:header = " num classified faces on gFace $gfTag = $count ... total classified = $numFaces" } - #get number of regions with a face classified on gface + #get number of regions with a face classified on gface mfaces/*/mregions/* { regionExists = function:mapI2I($id,"<-1 0>", 1) skip = function:notI($regionExists) - numRegions = function:addI(1,$numRegions) + numRegions = function:addI(1,$numRegions) } info:header = " num regions with a face classied on gFace $gfTag = $numRegions" } @@ -68,11 +68,11 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version info:header = " writing $numFaces faces .... " info:header = " writing $numRegions regions .... " - #write the header string + #write the header string header = "BOUNDARY CONDITIONS $version \n" + " $bcID 1 $numRegions 0 6" - #loop over the mesh faces classified on the geometric model faces with the bcAttribute + #loop over the mesh faces classified on the geometric model faces with the bcAttribute # and write the mesh elements and the face number gmodel/gfaces/* { gfTag = $tag @@ -90,38 +90,38 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version } skip = function:notI($hasBcAtt) - #need to determine the vertex numbering so the face numbering can be determined - # see www.stanford.edu/class/me469b/handouts/gambit_write.pdf pg C-33 - # for the tet ordering + #need to determine the vertex numbering so the face numbering can be determined + # see www.stanford.edu/class/me469b/handouts/gambit_write.pdf pg C-33 + # for the tet ordering # gambit and MeshSim have the same ordering for tets # see the R_face documentation for the MeshSim ordering mfaces/* { mRgnId_A = $(mregions[0]/id) mRgnId_B = $(mregions[1]/id) - . { + . { regionExists = function:mapI2I($mRgnId_A,"<-1 0>", 1) skip = function:notI($regionExists) # Need the raw mesh region id to call mfaceNum() fnum = function:mfaceNum($mRgnId_A) snum = function:addI($fnum, 1) # for each face write: - # - # <6 for tet> + # + # <6 for tet> # - header = " $mRgnId_A 6 $snum" + header = " $mRgnId_A 6 $snum" } # end region A - . { + . { regionExists = function:mapI2I($mRgnId_B,"<-1 0>", 1) skip = function:notI($regionExists) # Need the raw mesh region id to call mfaceNum() fnum = function:mfaceNum($mRgnId_B) snum = function:addI($fnum, 1) # for each face write: - # - # <6 for tet> + # + # <6 for tet> # - header = " $mRgnId_B 6 $snum" + header = " $mRgnId_B 6 $snum" } # end region B } # end faces } # end gfaces @@ -244,7 +244,7 @@ mesh/mregions { header = "$(id,I8) $(ntype,I2) $(ndp,I2) $(mvertices[*]/id,I8)" } . { - # Pyrmaid elements + # Pyrmaid elements skip = function:compI($ndp,5) header = "$(id,I8) $(ntype,I2) $(ndp,I2) " + "$(mvertices[0]/id,I8)$(mvertices[1]/id,I8)$(mvertices[3]/id,I8)$(mvertices[2]/id,I8)" + diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/e8/e8a0186ae71e0dac8f12bb52f2669169e15522d8.svn-base b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/e8/e8a0186ae71e0dac8f12bb52f2669169e15522d8.svn-base index b5cae27..9b454e5 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/e8/e8a0186ae71e0dac8f12bb52f2669169e15522d8.svn-base +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/.svn/pristine/e8/e8a0186ae71e0dac8f12bb52f2669169e15522d8.svn-base @@ -3,4 +3,3 @@ atd 13 image "seisSol" : "analysis" root { C "problem definition":"seisSol" 0 1 1; }; - diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/attDefs/v000/analysis.atd b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/attDefs/v000/analysis.atd index b5cae27..9b454e5 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/attDefs/v000/analysis.atd +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/attDefs/v000/analysis.atd @@ -3,4 +3,3 @@ atd 13 image "seisSol" : "analysis" root { C "problem definition":"seisSol" 0 1 1; }; - diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/export/case/v000/gambit_seisSol.sxp b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/export/case/v000/gambit_seisSol.sxp index 4afa7aa..8f58704 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/export/case/v000/gambit_seisSol.sxp +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/export/case/v000/gambit_seisSol.sxp @@ -17,11 +17,11 @@ define-procedure hasBoundaryCondition(attName) { gfTag = $tag attributes["boundaryCondition"]/* { match = function:strcmp($image,$(inputAttName)) - hasBcAtt = function:notI($match) + hasBcAtt = function:notI($match) } } - hasBoundaryCondition = $hasBcAtt + hasBoundaryCondition = $hasBcAtt } define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version) { @@ -55,11 +55,11 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version numFaces = function:addI($count,$numFaces) info:header = " num classified faces on gFace $gfTag = $count ... total classified = $numFaces" } - #get number of regions with a face classified on gface + #get number of regions with a face classified on gface mfaces/*/mregions/* { regionExists = function:mapI2I($id,"<-1 0>", 1) skip = function:notI($regionExists) - numRegions = function:addI(1,$numRegions) + numRegions = function:addI(1,$numRegions) } info:header = " num regions with a face classied on gFace $gfTag = $numRegions" } @@ -68,11 +68,11 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version info:header = " writing $numFaces faces .... " info:header = " writing $numRegions regions .... " - #write the header string + #write the header string header = "BOUNDARY CONDITIONS $version \n" + " $bcID 1 $numRegions 0 6" - #loop over the mesh faces classified on the geometric model faces with the bcAttribute + #loop over the mesh faces classified on the geometric model faces with the bcAttribute # and write the mesh elements and the face number gmodel/gfaces/* { gfTag = $tag @@ -90,38 +90,38 @@ define-procedure writeMeshFacesWithAttribute(attName, displayName, bcID, version } skip = function:notI($hasBcAtt) - #need to determine the vertex numbering so the face numbering can be determined - # see www.stanford.edu/class/me469b/handouts/gambit_write.pdf pg C-33 - # for the tet ordering + #need to determine the vertex numbering so the face numbering can be determined + # see www.stanford.edu/class/me469b/handouts/gambit_write.pdf pg C-33 + # for the tet ordering # gambit and MeshSim have the same ordering for tets # see the R_face documentation for the MeshSim ordering mfaces/* { mRgnId_A = $(mregions[0]/id) mRgnId_B = $(mregions[1]/id) - . { + . { regionExists = function:mapI2I($mRgnId_A,"<-1 0>", 1) skip = function:notI($regionExists) # Need the raw mesh region id to call mfaceNum() fnum = function:mfaceNum($mRgnId_A) snum = function:addI($fnum, 1) # for each face write: - # - # <6 for tet> + # + # <6 for tet> # - header = " $mRgnId_A 6 $snum" + header = " $mRgnId_A 6 $snum" } # end region A - . { + . { regionExists = function:mapI2I($mRgnId_B,"<-1 0>", 1) skip = function:notI($regionExists) # Need the raw mesh region id to call mfaceNum() fnum = function:mfaceNum($mRgnId_B) snum = function:addI($fnum, 1) # for each face write: - # - # <6 for tet> + # + # <6 for tet> # - header = " $mRgnId_B 6 $snum" + header = " $mRgnId_B 6 $snum" } # end region B } # end faces } # end gfaces @@ -244,7 +244,7 @@ mesh/mregions { header = "$(id,I8) $(ntype,I2) $(ndp,I2) $(mvertices[*]/id,I8)" } . { - # Pyrmaid elements + # Pyrmaid elements skip = function:compI($ndp,5) header = "$(id,I8) $(ntype,I2) $(ndp,I2) " + "$(mvertices[0]/id,I8)$(mvertices[1]/id,I8)$(mvertices[3]/id,I8)$(mvertices[2]/id,I8)" + diff --git a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/seisSol.adi b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/seisSol.adi index 709b35a..10208e1 100644 --- a/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/seisSol.adi +++ b/SimModelerDownloadingBuilding/SimModelerCustomization/export/simSeisSol/seisSol.adi @@ -4,4 +4,3 @@ SolverName=seisSol [0] AttDefs=simSeisSol/attDefs/v000/ PatternInfo=case/simSeisSol/export/case/v000/gambit_seisSol.sxi - diff --git a/SimModelerDownloadingBuilding/downloadSimLib.sh b/SimModelerDownloadingBuilding/downloadSimLib.sh old mode 100644 new mode 100755 index 9a4461c..655efea --- a/SimModelerDownloadingBuilding/downloadSimLib.sh +++ b/SimModelerDownloadingBuilding/downloadSimLib.sh @@ -18,27 +18,27 @@ str=$4 i=$((${#str}-3)) last3char="${str:$i:3}" -if [ "$last3char" = "dev" ]; then - fold='DM' +if [ "$last3char" = "dev" ]; then + fold='DM' else - fold='M' + fold='M' fi; # in vi add newlines # %s/<\/a>/<\/a>\r/g -# to source.php( source code of support download site ) then run the command +# to source.php( source code of support download site ) then run the command # grep documentation.*zip ~/Downloads/source.php | sed 's/.*documentation\/\([A-Za-z]*.zip\).*/\1/g' # to extract the zip file names -if [ "$downloadcode" = "1" ]; then +if [ "$downloadcode" = "1" ]; then documentation=( GeomSim.zip GeomSimAcis.zip GeomSimDiscrete.zip FieldSim.zip GeomSimAbstract.zip MeshSimCrack.zip MeshSimAdapt.zip MeshSimAdv.zip MeshSim.zip MeshSimCrack.zip ParallelMeshSimAdapt.zip ParallelMeshSim.zip GeomSimProe.zip GeomSimParasolid.zip GeomSimSolidWorks.zip ) echo $documentation else documentation=() fi; -for doc in ${documentation[@]}; do +for doc in ${documentation[@]}; do wget --user=${username} --password=${password} http://www.simmetrix.com/application/release/${fold}/${release}/documentation/${doc} done @@ -46,13 +46,12 @@ done # grep linux64.tgz ~/Downloads/source.php | sed 's/.*release\/\([a-z]*-linux64.tgz\).*/\1/g' # to extract the tarball names -if [ "$downloadcode" = "0" ]; then +if [ "$downloadcode" = "0" ]; then components=( gmcore-linux64.tgz aciskrnl-linux64.tgz discrete-linux64.tgz fdcore-linux64.tgz gmabstract-linux64.tgz gmadv-linux64.tgz msadapt-linux64.tgz msadv-linux64.tgz mscore-linux64.tgz msparalleladapt-linux64.tgz msparallelmesh-linux64.tgz pskrnl-linux64.tgz simlicense-linux64.tgz ) else components=() fi; -for comp in ${components[@]}; do +for comp in ${components[@]}; do wget --user=${username} --password=${password} http://www.simmetrix.com/application/release/${fold}/${release}/release/${comp} done - diff --git a/analyse_communication/SConstruct b/analyse_communication/SConstruct index 3127b7b..5588558 100644 --- a/analyse_communication/SConstruct +++ b/analyse_communication/SConstruct @@ -1,4 +1,3 @@ -#! /usr/bin/python ## # @file # This file is part of SeisSol. @@ -8,21 +7,21 @@ # @section LICENSE # Copyright (c) 2017, SeisSol Group # All rights reserved. -# +# # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: -# +# # 1. Redistributions of source code must retain the above copyright notice, # this list of conditions and the following disclaimer. -# +# # 2. Redistributions in binary form must reproduce the above copyright notice, # this list of conditions and the following disclaimer in the documentation # and/or other materials provided with the distribution. -# +# # 3. Neither the name of the copyright holder nor the names of its # contributors may be used to endorse or promote products derived from this # software without specific prior written permission. -# +# # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE diff --git a/analyse_communication/src/Graph.cpp b/analyse_communication/src/Graph.cpp index cc96bc1..fb01373 100644 --- a/analyse_communication/src/Graph.cpp +++ b/analyse_communication/src/Graph.cpp @@ -1,36 +1,30 @@ #include "Graph.h" +#include #include #include -#include -Graph::Graph(unsigned numberOfNodes) - : m_numberOfNodes(numberOfNodes) -{ +Graph::Graph(unsigned numberOfNodes) : m_numberOfNodes(numberOfNodes) { m_edges = new counter_t[numberOfNodes * numberOfNodes]; memset(m_edges, 0, numberOfNodes * numberOfNodes * sizeof(counter_t)); } -Graph::~Graph() -{ - delete[] m_edges; -} +Graph::~Graph() { delete[] m_edges; } -void Graph::addEdge(unsigned fromNode, unsigned toNode) -{ +void Graph::addEdge(unsigned fromNode, unsigned toNode) { ++m_edges[fromNode * m_numberOfNodes + toNode]; } -void Graph::printDOT(std::string const& filename) -{ +void Graph::printDOT(const std::string& filename) { FILE* fp = fopen(filename.c_str(), "w"); fprintf(fp, "graph communication {\n"); fprintf(fp, "overlap=scale;\n"); fprintf(fp, "splines=true;\n"); counter_t maxEdgeCut = 0; for (unsigned from = 0; from < m_numberOfNodes; ++from) { - for (unsigned to = from+1; to < m_numberOfNodes; ++to) { - maxEdgeCut = std::max(maxEdgeCut, m_edges[from * m_numberOfNodes + to] + m_edges[to * m_numberOfNodes + from]); + for (unsigned to = from + 1; to < m_numberOfNodes; ++to) { + maxEdgeCut = std::max( + maxEdgeCut, m_edges[from * m_numberOfNodes + to] + m_edges[to * m_numberOfNodes + from]); } } for (unsigned from = 0; from < m_numberOfNodes; ++from) { @@ -39,23 +33,24 @@ void Graph::printDOT(std::string const& filename) edgeCut += m_edges[from * m_numberOfNodes + to]; } fprintf(fp, "%u [label=\"P%u\n%u\"];\n", from, from, edgeCut); - - for (unsigned to = from+1; to < m_numberOfNodes; ++to) { - counter_t edgeCut = m_edges[from * m_numberOfNodes + to] + m_edges[to * m_numberOfNodes + from]; + + for (unsigned to = from + 1; to < m_numberOfNodes; ++to) { + counter_t edgeCut = + m_edges[from * m_numberOfNodes + to] + m_edges[to * m_numberOfNodes + from]; if (edgeCut > 0) { double edgeWeight = static_cast(edgeCut) / maxEdgeCut; - fprintf(fp, "%u -- %u [weight=%lf,penwidth=%lf];\n", from, to, edgeWeight, 7.0*edgeWeight); + fprintf( + fp, "%u -- %u [weight=%lf,penwidth=%lf];\n", from, to, edgeWeight, 7.0 * edgeWeight); } } } - + fprintf(fp, "}\n"); - + fclose(fp); } -void Graph::printMatrix(std::string const& filename) -{ +void Graph::printMatrix(const std::string& filename) { FILE* fp = fopen(filename.c_str(), "wb"); fwrite(m_edges, sizeof(counter_t), m_numberOfNodes * m_numberOfNodes, fp); fclose(fp); diff --git a/analyse_communication/src/Graph.h b/analyse_communication/src/Graph.h index 6533837..97856aa 100644 --- a/analyse_communication/src/Graph.h +++ b/analyse_communication/src/Graph.h @@ -5,18 +5,18 @@ typedef unsigned long long counter_t; -class Graph { -private: +class Graph { + private: counter_t* m_edges; unsigned m_numberOfNodes; -public: + public: explicit Graph(unsigned numberOfNodes); ~Graph(); - + void addEdge(unsigned fromNode, unsigned toNode); - void printDOT(std::string const& filename); - void printMatrix(std::string const& filename); + void printDOT(const std::string& filename); + void printMatrix(const std::string& filename); }; #endif // GRAPH_H_ diff --git a/analyse_communication/src/SConscript b/analyse_communication/src/SConscript index 920d55c..bbc09af 100644 --- a/analyse_communication/src/SConscript +++ b/analyse_communication/src/SConscript @@ -1,4 +1,3 @@ -#!/usr/bin/python ## # @file # This file is part of SeisSol. @@ -8,21 +7,21 @@ # @section LICENSE # Copyright (c) 2017, SeisSol Group # All rights reserved. -# +# # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: -# +# # 1. Redistributions of source code must retain the above copyright notice, # this list of conditions and the following disclaimer. -# +# # 2. Redistributions in binary form must reproduce the above copyright notice, # this list of conditions and the following disclaimer in the documentation # and/or other materials provided with the distribution. -# +# # 3. Neither the name of the copyright holder nor the names of its # contributors may be used to endorse or promote products derived from this # software without specific prior written permission. -# +# # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE diff --git a/analyse_communication/src/main.cpp b/analyse_communication/src/main.cpp index df69573..58cf048 100644 --- a/analyse_communication/src/main.cpp +++ b/analyse_communication/src/main.cpp @@ -1,13 +1,12 @@ +#include "Graph.h" +#include "common/PartitionReader.h" + +#include #include #include -#include #include -#include "common/PartitionReader.h" -#include "Graph.h" - -void printCounter(std::string const& name, counter_t* counter, unsigned numPartitions) -{ +void printCounter(const std::string& name, counter_t* counter, unsigned numPartitions) { std::cout << std::endl << name; for (unsigned p = 0; p < numPartitions; ++p) { if (p % 10 == 0) { @@ -18,8 +17,7 @@ void printCounter(std::string const& name, counter_t* counter, unsigned numParti std::cout << std::endl << std::endl; } -int main(int argc, char** argv) -{ +int main(int argc, char** argv) { utils::Args args; args.addOption("dot", 'd', "Write dot file for visualisation.", utils::Args::Required, false); args.addOption("matrix", 'm', "Write edge-cut matrix.", utils::Args::Required, false); @@ -34,64 +32,70 @@ int main(int argc, char** argv) std::string matrixFile = args.getArgument("matrix", ""); PartitionReader reader(meshFile); - + counter_t totalEdgeCut = 0; counter_t totalCommVolume = 0; - + counter_t* edgeCut = new counter_t[reader.partitions]; counter_t* commVolume = new counter_t[reader.partitions]; - - memset(edgeCut, 0, reader.partitions*sizeof(counter_t)); - memset(commVolume, 0, reader.partitions*sizeof(counter_t)); - + + memset(edgeCut, 0, reader.partitions * sizeof(counter_t)); + memset(commVolume, 0, reader.partitions * sizeof(counter_t)); + Graph edgeCutGraph(reader.partitions); for (unsigned p = 0; p < reader.partitions; ++p) { - if (p%10 == 0) { + if (p % 10 == 0) { std::cout << "Reading partition " << p << std::endl; } - + reader.readPartition(p); - + for (unsigned elem = 0; elem < reader.elementSize[p]; ++elem) { for (unsigned face = 0; face < 4; ++face) { - int neighborRank = reader.elementNeighborRanks[4*elem + face]; + int neighborRank = reader.elementNeighborRanks[4 * elem + face]; if (neighborRank != p) { ++edgeCut[p]; edgeCutGraph.addEdge(p, static_cast(neighborRank)); } } - std::set neighbors(reader.elementNeighborRanks + 4*elem, reader.elementNeighborRanks + 4*elem + 4); + std::set neighbors(reader.elementNeighborRanks + 4 * elem, + reader.elementNeighborRanks + 4 * elem + 4); neighbors.erase(p); commVolume[p] += neighbors.size(); } - + totalEdgeCut += edgeCut[p]; totalCommVolume += commVolume[p]; } - + if (dotFile.size() > 0) { edgeCutGraph.printDOT(dotFile); } if (matrixFile.size() > 0) { edgeCutGraph.printMatrix(matrixFile); } - + printCounter("Edge cut", edgeCut, reader.partitions); printCounter("Communication volume", commVolume, reader.partitions); - + std::cout << "Total edge cut: " << totalEdgeCut << std::endl; std::cout << "Total communication volume: " << totalCommVolume << std::endl; - std::cout << "Ratio: " << static_cast(totalEdgeCut) / totalCommVolume << std::endl << std::endl; - + std::cout << "Ratio: " << static_cast(totalEdgeCut) / totalCommVolume << std::endl + << std::endl; + std::cout << "Order\t\tEdge cut (MB)\tComm volume (MB)" << std::endl; for (unsigned order = 1; order <= 8; ++order) { - std::cout << order << "\t\t" << sizeof(double)*9*order*(order+1)/2.0 * totalEdgeCut / (1024.0*1024.0) - << "\t\t" << sizeof(double)*9*order*(order+1)*(order+2)/6.0 * totalCommVolume / (1024.0*1024.0) << std::endl; + std::cout << order << "\t\t" + << sizeof(double) * 9 * order * (order + 1) / 2.0 * totalEdgeCut / (1024.0 * 1024.0) + << "\t\t" + << sizeof(double) * 9 * order * (order + 1) * (order + 2) / 6.0 * totalCommVolume / + (1024.0 * 1024.0) + << std::endl; } - + delete[] edgeCut; delete[] commVolume; - + return 0; } diff --git a/check_mesh/CMakeLists.txt b/check_mesh/CMakeLists.txt index 260153e..dfadf5d 100644 --- a/check_mesh/CMakeLists.txt +++ b/check_mesh/CMakeLists.txt @@ -22,4 +22,3 @@ target_compile_definitions(check_mesh PUBLIC USE_MPI) find_package(HDF5 REQUIRED COMPONENTS C HL) target_include_directories(check_mesh PUBLIC ${HDF5_INCLUDE_DIRS}) target_link_libraries(check_mesh PUBLIC ${HDF5_C_HL_LIBRARIES} ${HDF5_C_LIBRARIES}) - diff --git a/check_mesh/README.md b/check_mesh/README.md index ae1d350..a3bb956 100644 --- a/check_mesh/README.md +++ b/check_mesh/README.md @@ -1,8 +1,8 @@ # Evaluate material model on mesh -## Compilation +## Building -``` +```bash mkdir build && cd build cmake .. make @@ -10,26 +10,31 @@ make ## Usage -``` +```bash mpirun -n 4 ./check_mesh -m mesh.h5 ``` ## Description 1. Read `mesh.h5` -2. For each element check that: - - if a face is an internal face, the neighbor *has to exist*. - - if a face is an external face, the neighbor *must not exist*. +2. For each element check that: the neighbor exists +if and _only_ if the face is an internal face. 3. If an element breaks above rules, the program writes an error message. -4. If all elements are correct, the programm exits with a success output. If at least one element is broken, the program returns with an error code. +4. If all elements are correct, the programm exits with +a success output. If at least one element is broken, +the program returns with an error code. Internal faces are: -* Regular -* Dynamic Rupture -* Dynamic Rupture with Fault Tagging + +- Regular (0) +- Dynamic Rupture (3) +- Dynamic Rupture with Fault Tagging (64 or higher) External faces are: -* Free Surface -* Absorbing -* Periodic -* Free Surface with Gravity + +- Free Surface (1) +- Absorbing (2) +- Periodic (6) +- Free Surface with Gravity (5) +- Dirichlet (4) +- Analytical (7) diff --git a/check_mesh/src/main.cpp b/check_mesh/src/main.cpp index 5e54929..1f1a3b3 100644 --- a/check_mesh/src/main.cpp +++ b/check_mesh/src/main.cpp @@ -1,12 +1,11 @@ -#include -#include - -#include - #include "mesh.h" #include "utils/args.h" #include "utils/logger.h" +#include +#include +#include + int main(int argc, char** argv) { MPI_Init(&argc, &argv); int rank; diff --git a/check_mesh/src/mesh.cpp b/check_mesh/src/mesh.cpp index 39b95cc..d8b5417 100644 --- a/check_mesh/src/mesh.cpp +++ b/check_mesh/src/mesh.cpp @@ -1,13 +1,12 @@ #include "mesh.h" -#include - #include +#include #pragma clang diagnostic push #pragma clang diagnostic ignored "-Weverything" -#include "PUML/Neighbor.h" #include "PUML/Downward.h" +#include "PUML/Neighbor.h" #pragma clang diagnostic pop Mesh::Mesh(const std::string& fileName) { @@ -18,7 +17,7 @@ Mesh::Mesh(const std::string& fileName) { std::vector cellIdsAsInFile(numTotalCells); std::iota(cellIdsAsInFile.begin(), cellIdsAsInFile.end(), 0); puml.addData(cellIdsAsInFile.data(), numTotalCells, PUML::CELL); - + // Keep all cells on the same rank int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); @@ -40,7 +39,7 @@ bool Mesh::checkNeighbors() const { const auto& cell = puml.cells().at(elemIdx); const auto bndInfo = puml.cellData(1)[elemIdx]; const auto cellIdAsInFile = puml.cellData(2)[elemIdx]; - auto decodeBC = [&bndInfo](std::size_t side) { return (bndInfo >> (side*8)) & 0xFF; }; + auto decodeBC = [&bndInfo](std::size_t side) { return (bndInfo >> (side * 8)) & 0xFF; }; PUML::Downward::faces(puml, cell, faceids); PUML::Neighbor::face(puml, elemIdx, neighbors); @@ -49,23 +48,31 @@ bool Mesh::checkNeighbors() const { for (size_t side = 0; side < 4; side++) { const auto& face = puml.faces()[faceids[side]]; const auto sideBC = decodeBC(side); - // if a face is an internal face, it has to have a neighbor on either this rank or somewhere else: + // if a face is an internal face, it has to have a neighbor on either this + // rank or somewhere else: if (bcToType(sideBC) == BCType::internal) { if (neighbors[side] < 0 && !face.isShared()) { - logInfo() << "Element" << cellIdAsInFile << ", side" << side << " has a" << bcToString(sideBC) << "boundary condition, but the neighboring element doesn't exist"; + logInfo() << "Element" << cellIdAsInFile << ", side" << side << " has a" + << bcToString(sideBC) + << "boundary condition, but the neighboring element " + "doesn't exist"; result = false; } } // external boundaries must not have neighboring elements: else if (bcToType(sideBC) == BCType::external) { if (neighbors[side] >= 0 || face.isShared()) { - logInfo() << "Element" << cellIdAsInFile << ", side" << side << " has a" << bcToString(sideBC) << "boundary condition, but the neighboring element is not flagged -1"; + logInfo() << "Element" << cellIdAsInFile << ", side" << side << " has a" + << bcToString(sideBC) + << "boundary condition, but the neighboring element is not " + "flagged -1"; result = false; } } // ignore unknown boundary conditions and warn else { - logWarning() << "Element" << cellIdAsInFile << ", side" << side << " has a boundary condition, which I don't understand" << sideBC; + logWarning() << "Element" << cellIdAsInFile << ", side" << side + << " has a boundary condition, which I don't understand" << sideBC; } } } diff --git a/check_mesh/src/mesh.h b/check_mesh/src/mesh.h index 8385ea0..4e184ec 100644 --- a/check_mesh/src/mesh.h +++ b/check_mesh/src/mesh.h @@ -1,18 +1,17 @@ #ifndef MESH_H_ #define MESH_H_ +#include "PUML/PUML.h" + #include #include #include #include #include -#include "PUML/PUML.h" - - class Mesh { public: - explicit Mesh(std::string const& fileName); + explicit Mesh(const std::string& fileName); ~Mesh() = default; [[nodiscard]] std::size_t numVertices() const { return puml.vertices().size(); } [[nodiscard]] std::size_t numElements() const { return puml.cells().size(); } @@ -21,9 +20,7 @@ class Mesh { private: PUML::TETPUML puml; - enum class BCType { - internal, external, unknown - }; + enum class BCType { internal, external, unknown }; BCType bcToType(int id) const { if (id == 0 || id == 3 || id > 64) { @@ -36,19 +33,26 @@ class Mesh { } std::string bcToString(int id) const { - if (id == 0) { return std::string("regular"); } - else if (id == 1) { return std::string("free surface"); } - else if (id == 2) { return std::string("free surface with gravity"); } - else if (id == 3) { return std::string("dynamic rupture"); } - else if (id == 5) { return std::string("absorbing"); } - else if (id == 6) { return std::string("periodic"); } - else if (id > 64) { + if (id == 0) { + return std::string("regular"); + } else if (id == 1) { + return std::string("free surface"); + } else if (id == 2) { + return std::string("free surface with gravity"); + } else if (id == 3) { + return std::string("dynamic rupture"); + } else if (id == 5) { + return std::string("absorbing"); + } else if (id == 6) { + return std::string("periodic"); + } else if (id > 64) { std::stringstream s; s << "fault-tagging (" << id << ")"; return s.str(); - } else {return std::string(""); } + } else { + return std::string(""); + } } - }; #endif diff --git a/common/PartitionReader.cpp b/common/PartitionReader.cpp index 782bbf6..a3bfd51 100644 --- a/common/PartitionReader.cpp +++ b/common/PartitionReader.cpp @@ -2,22 +2,23 @@ * @file * This file is part of SeisSol. * - * @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) + * @author Carsten Uphoff (c.uphoff AT tum.de, + * http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) * * @section LICENSE * Copyright (c) 2017, SeisSol Group * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, * this list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from this * software without specific prior written permission. @@ -38,86 +39,86 @@ **/ #include "PartitionReader.h" -#include #include +#include #include -static void check_err(const int stat, const int line, const char *file) { +static void check_err(const int stat, const int line, const char* file) { if (stat != NC_NOERR) { - fprintf(stderr,"line %d of %s: %s\n", line, file, nc_strerror(stat)); + fprintf(stderr, "line %d of %s: %s\n", line, file, nc_strerror(stat)); fflush(stderr); exit(-1); } } -PartitionReader::PartitionReader(std::string const& fileName){ +PartitionReader::PartitionReader(const std::string& fileName) { int stat; - size_t maxVertices; - size_t maxElements; + size_t maxVertices; + size_t maxElements; stat = nc_open(fileName.c_str(), 0, &ncid); - check_err(stat,__LINE__,__FILE__); - - int ncDimPart; - stat = nc_inq_dimid(ncid, "partitions", &ncDimPart); - check_err(stat,__LINE__,__FILE__); - stat = nc_inq_dimlen(ncid, ncDimPart, &partitions); - check_err(stat,__LINE__,__FILE__); - - int ncDimVrtx; - stat = nc_inq_dimid(ncid, "vertices", &ncDimVrtx); - check_err(stat,__LINE__,__FILE__); - stat = nc_inq_dimlen(ncid, ncDimVrtx, &maxVertices); - check_err(stat,__LINE__,__FILE__); - - int ncDimElem; - stat = nc_inq_dimid(ncid, "elements", &ncDimElem); - check_err(stat,__LINE__,__FILE__); - stat = nc_inq_dimlen(ncid, ncDimElem, &maxElements); - check_err(stat,__LINE__,__FILE__); + check_err(stat, __LINE__, __FILE__); + + int ncDimPart; + stat = nc_inq_dimid(ncid, "partitions", &ncDimPart); + check_err(stat, __LINE__, __FILE__); + stat = nc_inq_dimlen(ncid, ncDimPart, &partitions); + check_err(stat, __LINE__, __FILE__); + + int ncDimVrtx; + stat = nc_inq_dimid(ncid, "vertices", &ncDimVrtx); + check_err(stat, __LINE__, __FILE__); + stat = nc_inq_dimlen(ncid, ncDimVrtx, &maxVertices); + check_err(stat, __LINE__, __FILE__); + + int ncDimElem; + stat = nc_inq_dimid(ncid, "elements", &ncDimElem); + check_err(stat, __LINE__, __FILE__); + stat = nc_inq_dimlen(ncid, ncDimElem, &maxElements); + check_err(stat, __LINE__, __FILE__); int ncVarElemSize; - stat = nc_inq_varid(ncid, "element_size", &ncVarElemSize); - check_err(stat,__LINE__,__FILE__); + stat = nc_inq_varid(ncid, "element_size", &ncVarElemSize); + check_err(stat, __LINE__, __FILE__); - stat = nc_inq_varid(ncid, "element_vertices", &ncVarElemVertices); - check_err(stat,__LINE__,__FILE__); + stat = nc_inq_varid(ncid, "element_vertices", &ncVarElemVertices); + check_err(stat, __LINE__, __FILE__); - stat = nc_inq_varid(ncid, "element_boundaries", &ncVarElemBoundaries); - check_err(stat,__LINE__,__FILE__); + stat = nc_inq_varid(ncid, "element_boundaries", &ncVarElemBoundaries); + check_err(stat, __LINE__, __FILE__); - stat = nc_inq_varid(ncid, "element_neighbor_ranks", &ncVarElemNeighborRanks); - check_err(stat,__LINE__,__FILE__); + stat = nc_inq_varid(ncid, "element_neighbor_ranks", &ncVarElemNeighborRanks); + check_err(stat, __LINE__, __FILE__); - stat = nc_inq_varid(ncid, "element_group", &ncVarElemGroup); - check_err(stat,__LINE__,__FILE__); + stat = nc_inq_varid(ncid, "element_group", &ncVarElemGroup); + check_err(stat, __LINE__, __FILE__); int ncVarVrtxSize; - stat = nc_inq_varid(ncid, "vertex_size", &ncVarVrtxSize); - check_err(stat,__LINE__,__FILE__); - - stat = nc_inq_varid(ncid, "vertex_coordinates", &ncVarVrtxCoords); - check_err(stat,__LINE__,__FILE__); - - elementSize = new int[partitions]; - stat = nc_get_var_int(ncid, ncVarElemSize, elementSize); - check_err(stat,__LINE__,__FILE__); - - vertexSize = new int[partitions]; - stat = nc_get_var_int(ncid, ncVarVrtxSize, vertexSize); - check_err(stat,__LINE__,__FILE__); - - vertexCoordinates = new double[3*maxVertices]; - elementBoundaries = new int[4*maxElements]; - elementVertices = new int[4*maxElements]; - elementNeighborRanks = new int[4*maxElements]; + stat = nc_inq_varid(ncid, "vertex_size", &ncVarVrtxSize); + check_err(stat, __LINE__, __FILE__); + + stat = nc_inq_varid(ncid, "vertex_coordinates", &ncVarVrtxCoords); + check_err(stat, __LINE__, __FILE__); + + elementSize = new int[partitions]; + stat = nc_get_var_int(ncid, ncVarElemSize, elementSize); + check_err(stat, __LINE__, __FILE__); + + vertexSize = new int[partitions]; + stat = nc_get_var_int(ncid, ncVarVrtxSize, vertexSize); + check_err(stat, __LINE__, __FILE__); + + vertexCoordinates = new double[3 * maxVertices]; + elementBoundaries = new int[4 * maxElements]; + elementVertices = new int[4 * maxElements]; + elementNeighborRanks = new int[4 * maxElements]; elementGroup = new int[maxElements]; } PartitionReader::~PartitionReader() { int stat = nc_close(ncid); - check_err(stat,__LINE__,__FILE__); - + check_err(stat, __LINE__, __FILE__); + delete[] elementSize; delete[] vertexSize; delete[] elementBoundaries; @@ -131,23 +132,23 @@ void PartitionReader::readPartition(int partition) { int stat; size_t elementsStart[3] = {partition, 0, 0}; size_t elementsSize[3] = {1, elementSize[partition], 4}; - + stat = nc_get_vara_int(ncid, ncVarElemVertices, elementsStart, elementsSize, elementVertices); - check_err(stat,__LINE__,__FILE__); - + check_err(stat, __LINE__, __FILE__); + stat = nc_get_vara_int(ncid, ncVarElemBoundaries, elementsStart, elementsSize, elementBoundaries); - check_err(stat,__LINE__,__FILE__); - - stat = nc_get_vara_int(ncid, ncVarElemNeighborRanks, elementsStart, elementsSize, elementNeighborRanks); - check_err(stat,__LINE__,__FILE__); - + check_err(stat, __LINE__, __FILE__); + + stat = nc_get_vara_int( + ncid, ncVarElemNeighborRanks, elementsStart, elementsSize, elementNeighborRanks); + check_err(stat, __LINE__, __FILE__); + stat = nc_get_vara_int(ncid, ncVarElemGroup, elementsStart, elementsSize, elementGroup); - check_err(stat,__LINE__,__FILE__); - - + check_err(stat, __LINE__, __FILE__); + size_t verticesStart[3] = {partition, 0, 0}; size_t verticesSize[3] = {1, vertexSize[partition], 3}; - + stat = nc_get_vara_double(ncid, ncVarVrtxCoords, verticesStart, verticesSize, vertexCoordinates); - check_err(stat,__LINE__,__FILE__); + check_err(stat, __LINE__, __FILE__); } diff --git a/common/PartitionReader.h b/common/PartitionReader.h index abb4265..0fae5ee 100644 --- a/common/PartitionReader.h +++ b/common/PartitionReader.h @@ -2,22 +2,23 @@ * @file * This file is part of SeisSol. * - * @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) + * @author Carsten Uphoff (c.uphoff AT tum.de, + * http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) * * @section LICENSE * Copyright (c) 2017, SeisSol Group * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, * this list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from this * software without specific prior written permission. @@ -42,13 +43,13 @@ #include class PartitionReader { -public: - explicit PartitionReader(std::string const& fileName); + public: + explicit PartitionReader(const std::string& fileName); ~PartitionReader(); - + void readPartition(int partition); - size_t partitions; + size_t partitions; int* elementSize; int* vertexSize; int* elementBoundaries; @@ -56,14 +57,14 @@ class PartitionReader { int* elementNeighborRanks; int* elementGroup; double* vertexCoordinates; - -private: + + private: int ncid; - int ncVarElemVertices; - int ncVarElemBoundaries; - int ncVarElemNeighborRanks; + int ncVarElemVertices; + int ncVarElemBoundaries; + int ncVarElemNeighborRanks; int ncVarElemGroup; - int ncVarVrtxCoords; + int ncVarVrtxCoords; }; #endif // PARTITIONREADER_H_ diff --git a/creating_geometric_models/ExampleFiles/Gocad_workflow/ExampleCurve_createFaultFromCurve.pl b/creating_geometric_models/ExampleFiles/Gocad_workflow/ExampleCurve_createFaultFromCurve.pl index f38ea67..a47f41a 100644 --- a/creating_geometric_models/ExampleFiles/Gocad_workflow/ExampleCurve_createFaultFromCurve.pl +++ b/creating_geometric_models/ExampleFiles/Gocad_workflow/ExampleCurve_createFaultFromCurve.pl @@ -1,9 +1,9 @@ -GOCAD PLine 1 +GOCAD PLine 1 HEADER { name: test } GOCAD_ORIGINAL_COORDINATE_SYSTEM -NAME " gocad Local" +NAME " gocad Local" PROJECTION Unknown DATUM Unknown AXIS_NAME X Y Z @@ -24,57 +24,57 @@ is_z: on } ILINE -VRTX 1 6169194.75 -3927111.59375 3982.27001953125 -VRTX 2 6169012.765625 -3928064.78125 3982.27001953125 -VRTX 3 6168830.78125 -3929017.96875 3982.27001953125 -VRTX 4 6168648.96875 -3929970.5 3982.27001953125 -VRTX 5 6168467.140625 -3930923.03125 3982.27001953125 -VRTX 6 6168285.375 -3931875.875 3982.27001953125 -VRTX 7 6168103.59375 -3932828.71875 3982.27001953125 -VRTX 8 6167941.84375 -3933787.96875 3982.27001953125 -VRTX 9 6167780.09375 -3934747.21875 3982.27001953125 -VRTX 10 6167621.46875 -3935706 3982.27001953125 -VRTX 11 6167462.828125 -3936664.75 3982.27001953125 -VRTX 12 6167304.1875 -3937623.71875 3982.27001953125 -VRTX 13 6167145.546875 -3938582.6875 3982.27001953125 -VRTX 14 6166840.125 -3939505.28125 3982.27001953125 -VRTX 15 6166534.71875 -3940427.875 3982.27001953125 -VRTX 16 6166224.6875 -3941349.8125 3982.27001953125 -VRTX 17 6165914.65625 -3942271.75 3982.27001953125 -VRTX 18 6165604.625 -3943193.9375 3982.27001953125 -VRTX 19 6165294.578125 -3944116.09375 3982.27001953125 -VRTX 20 6164941.75 -3945022.90625 3982.27001953125 -VRTX 21 6164588.9375 -3945929.71875 3982.27001953125 -VRTX 22 6164224.796875 -3946832.15625 3982.27001953125 -VRTX 23 6163860.65625 -3947734.59375 3982.27001953125 -VRTX 24 6163503.25 -3948639.8125 3982.27001953125 -VRTX 25 6163145.84375 -3949545.0625 3982.27001953125 -VRTX 26 6162799.21875 -3950455.53125 3982.27001953125 -VRTX 27 6162452.578125 -3951366 3982.27001953125 -SEG 1 2 -SEG 2 3 -SEG 3 4 -SEG 4 5 -SEG 5 6 -SEG 6 7 -SEG 7 8 -SEG 8 9 -SEG 9 10 -SEG 10 11 -SEG 11 12 -SEG 12 13 -SEG 13 14 -SEG 14 15 -SEG 15 16 -SEG 16 17 -SEG 17 18 -SEG 18 19 -SEG 19 20 -SEG 20 21 -SEG 21 22 -SEG 22 23 -SEG 23 24 -SEG 24 25 -SEG 25 26 -SEG 26 27 +VRTX 1 6169194.75 -3927111.59375 3982.27001953125 +VRTX 2 6169012.765625 -3928064.78125 3982.27001953125 +VRTX 3 6168830.78125 -3929017.96875 3982.27001953125 +VRTX 4 6168648.96875 -3929970.5 3982.27001953125 +VRTX 5 6168467.140625 -3930923.03125 3982.27001953125 +VRTX 6 6168285.375 -3931875.875 3982.27001953125 +VRTX 7 6168103.59375 -3932828.71875 3982.27001953125 +VRTX 8 6167941.84375 -3933787.96875 3982.27001953125 +VRTX 9 6167780.09375 -3934747.21875 3982.27001953125 +VRTX 10 6167621.46875 -3935706 3982.27001953125 +VRTX 11 6167462.828125 -3936664.75 3982.27001953125 +VRTX 12 6167304.1875 -3937623.71875 3982.27001953125 +VRTX 13 6167145.546875 -3938582.6875 3982.27001953125 +VRTX 14 6166840.125 -3939505.28125 3982.27001953125 +VRTX 15 6166534.71875 -3940427.875 3982.27001953125 +VRTX 16 6166224.6875 -3941349.8125 3982.27001953125 +VRTX 17 6165914.65625 -3942271.75 3982.27001953125 +VRTX 18 6165604.625 -3943193.9375 3982.27001953125 +VRTX 19 6165294.578125 -3944116.09375 3982.27001953125 +VRTX 20 6164941.75 -3945022.90625 3982.27001953125 +VRTX 21 6164588.9375 -3945929.71875 3982.27001953125 +VRTX 22 6164224.796875 -3946832.15625 3982.27001953125 +VRTX 23 6163860.65625 -3947734.59375 3982.27001953125 +VRTX 24 6163503.25 -3948639.8125 3982.27001953125 +VRTX 25 6163145.84375 -3949545.0625 3982.27001953125 +VRTX 26 6162799.21875 -3950455.53125 3982.27001953125 +VRTX 27 6162452.578125 -3951366 3982.27001953125 +SEG 1 2 +SEG 2 3 +SEG 3 4 +SEG 4 5 +SEG 5 6 +SEG 6 7 +SEG 7 8 +SEG 8 9 +SEG 9 10 +SEG 10 11 +SEG 11 12 +SEG 12 13 +SEG 13 14 +SEG 14 15 +SEG 15 16 +SEG 16 17 +SEG 17 18 +SEG 18 19 +SEG 19 20 +SEG 20 21 +SEG 21 22 +SEG 22 23 +SEG 23 24 +SEG 24 25 +SEG 25 26 +SEG 26 27 END diff --git a/creating_geometric_models/ExampleFiles/Gocad_workflow/ExampleDipType1_createFaultFromCurve.dat b/creating_geometric_models/ExampleFiles/Gocad_workflow/ExampleDipType1_createFaultFromCurve.dat index 4831626..8a76f36 100644 --- a/creating_geometric_models/ExampleFiles/Gocad_workflow/ExampleDipType1_createFaultFromCurve.dat +++ b/creating_geometric_models/ExampleFiles/Gocad_workflow/ExampleDipType1_createFaultFromCurve.dat @@ -1,4 +1,4 @@ --40e3 30 +-40e3 30 -25e3 30 -15e3 40 -5e3 50 diff --git a/creating_geometric_models/ExampleFiles/Gocad_workflow/TopBottom.pl b/creating_geometric_models/ExampleFiles/Gocad_workflow/TopBottom.pl index a5135fb..1ec44cf 100644 --- a/creating_geometric_models/ExampleFiles/Gocad_workflow/TopBottom.pl +++ b/creating_geometric_models/ExampleFiles/Gocad_workflow/TopBottom.pl @@ -11,8 +11,8 @@ VRTX 2 -60e3 60e3 -60e3 VRTX 3 -60e3 -60e3 -60e3 VRTX 4 60e3 -60e3 -60e3 -SEG 1 2 -SEG 2 3 +SEG 1 2 +SEG 2 3 SEG 3 4 SEG 4 1 END @@ -25,8 +25,8 @@ END VRTX 2 -60e3 60e3 -0e3 VRTX 3 -60e3 -60e3 -0e3 VRTX 4 60e3 -60e3 -0e3 -SEG 1 2 -SEG 2 3 +SEG 1 2 +SEG 2 3 SEG 3 4 SEG 4 1 END diff --git a/creating_geometric_models/create_fault_from_trace.py b/creating_geometric_models/create_fault_from_trace.py index 83789a8..764faef 100755 --- a/creating_geometric_models/create_fault_from_trace.py +++ b/creating_geometric_models/create_fault_from_trace.py @@ -296,4 +296,3 @@ def generate_vertices(maxdepth, sign=1): fout.write("TRGL %d %d %d\n" % (i + j * NX, i + 1 + (j + 1) * NX, i + (j + 1) * NX)) fout.write("END") print(f"done writing {fname}") - diff --git a/estimate_runtime/SConstruct b/estimate_runtime/SConstruct index 7e2711c..dd3c428 100644 --- a/estimate_runtime/SConstruct +++ b/estimate_runtime/SConstruct @@ -1,4 +1,3 @@ -#! /usr/bin/python ## # @file # This file is part of SeisSol. @@ -8,21 +7,21 @@ # @section LICENSE # Copyright (c) 2017, SeisSol Group # All rights reserved. -# +# # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: -# +# # 1. Redistributions of source code must retain the above copyright notice, # this list of conditions and the following disclaimer. -# +# # 2. Redistributions in binary form must reproduce the above copyright notice, # this list of conditions and the following disclaimer in the documentation # and/or other materials provided with the distribution. -# +# # 3. Neither the name of the copyright holder nor the names of its # contributors may be used to endorse or promote products derived from this # software without specific prior written permission. -# +# # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE diff --git a/estimate_runtime/src/SConscript b/estimate_runtime/src/SConscript index 811a599..5afd646 100644 --- a/estimate_runtime/src/SConscript +++ b/estimate_runtime/src/SConscript @@ -1,4 +1,3 @@ -#!/usr/bin/python ## # @file # This file is part of SeisSol. @@ -8,21 +7,21 @@ # @section LICENSE # Copyright (c) 2017, SeisSol Group # All rights reserved. -# +# # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: -# +# # 1. Redistributions of source code must retain the above copyright notice, # this list of conditions and the following disclaimer. -# +# # 2. Redistributions in binary form must reproduce the above copyright notice, # this list of conditions and the following disclaimer in the documentation # and/or other materials provided with the distribution. -# +# # 3. Neither the name of the copyright holder nor the names of its # contributors may be used to endorse or promote products derived from this # software without specific prior written permission. -# +# # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE diff --git a/estimate_runtime/src/main.cpp b/estimate_runtime/src/main.cpp index 6dc940f..27221fb 100644 --- a/estimate_runtime/src/main.cpp +++ b/estimate_runtime/src/main.cpp @@ -1,26 +1,24 @@ #include #include #define GLM_FORCE_SWIZZLE -#include +#include "common/PartitionReader.h" + #include +#include #include -#include "common/PartitionReader.h" - -double calculateTimestep(double cfl, double insphereRadius, double pVelocity, unsigned order) -{ - return cfl * 2.0 * insphereRadius / (pVelocity) / (2*(order-1)+1); +double calculateTimestep(double cfl, double insphereRadius, double pVelocity, unsigned order) { + return cfl * 2.0 * insphereRadius / (pVelocity) / (2 * (order - 1) + 1); } -unsigned getCluster(double timestep, double globalMinTimestep, unsigned rate) -{ +unsigned getCluster(double timestep, double globalMinTimestep, unsigned rate) { if (rate == 1) { return 0; } double upper; upper = rate * globalMinTimestep; - + unsigned cluster = 0; while (upper <= timestep) { upper *= rate; @@ -29,8 +27,7 @@ unsigned getCluster(double timestep, double globalMinTimestep, unsigned rate) return cluster; } -int main(int argc, char** argv) -{ +int main(int argc, char** argv) { utils::Args args; args.addOption("mesh", 'm', "Netcdf mesh file"); args.addOption("average-flops-per-cell", 'f', "Hardware flops per cell"); @@ -41,11 +38,11 @@ int main(int argc, char** argv) args.addOption("rate", 'r', "Clusterd lts rate", utils::Args::Required, false); args.addOption("final-time", 't', "Final simulation time", utils::Args::Required, false); - if (args.parse(argc, argv) != utils::Args::Success) { + if (args.parse(argc, argv) != utils::Args::Success) { return -1; } - - std::string meshFile = args.getArgument("mesh"); + + std::string meshFile = args.getArgument("mesh"); double flopsPerCell = args.getArgument("average-flops-per-cell"); double nodePerformance = args.getArgument("node-performance"); double pVelocity = args.getArgument("p-wave-velocity", 6000.0); @@ -53,7 +50,7 @@ int main(int argc, char** argv) unsigned order = args.getArgument("order"); unsigned rate = args.getArgument("rate", 2); double finalTime = args.getArgument("final-time", 1.0); - + PartitionReader reader(meshFile); std::vector timesteps; @@ -62,7 +59,7 @@ int main(int argc, char** argv) double minTimestep = std::numeric_limits::max(); double maxTimestep = std::numeric_limits::min(); for (unsigned p = 0; p < reader.partitions; ++p) { - std::cout << "Reading partition " << p+1 << "..." << std::endl; + std::cout << "Reading partition " << p + 1 << "..." << std::endl; reader.readPartition(p); numberOfElements += reader.elementSize[p]; @@ -70,53 +67,62 @@ int main(int argc, char** argv) glm::dvec3 x[4]; for (unsigned vtx = 0; vtx < 4; ++vtx) { for (unsigned d = 0; d < 3; ++d) { - x[vtx][d] = reader.vertexCoordinates[3 * reader.elementVertices[4*elem + vtx] + d]; + x[vtx][d] = reader.vertexCoordinates[3 * reader.elementVertices[4 * elem + vtx] + d]; } } - double alpha = determinant(glm::dmat4(glm::dvec4(x[0], 1.0), glm::dvec4(x[1], 1.0), glm::dvec4(x[2], 1.0), glm::dvec4(x[3], 1.0))); - double Nabc = length(cross(x[1]-x[0], x[2]-x[0])); - double Nabd = length(cross(x[1]-x[0], x[3]-x[0])); - double Nacd = length(cross(x[2]-x[0], x[3]-x[0])); - double Nbcd = length(cross(x[2]-x[1], x[3]-x[1])); + double alpha = determinant(glm::dmat4(glm::dvec4(x[0], 1.0), + glm::dvec4(x[1], 1.0), + glm::dvec4(x[2], 1.0), + glm::dvec4(x[3], 1.0))); + double Nabc = length(cross(x[1] - x[0], x[2] - x[0])); + double Nabd = length(cross(x[1] - x[0], x[3] - x[0])); + double Nacd = length(cross(x[2] - x[0], x[3] - x[0])); + double Nbcd = length(cross(x[2] - x[1], x[3] - x[1])); double insphere = std::fabs(alpha) / (Nabc + Nabd + Nacd + Nbcd); minInsphereRadius = std::min(minInsphereRadius, insphere); - + double timestep = calculateTimestep(cfl, insphere, pVelocity, order); minTimestep = std::min(minTimestep, timestep); maxTimestep = std::max(maxTimestep, timestep); - + timesteps.push_back(timestep); } } - + unsigned numClusters = getCluster(maxTimestep, minTimestep, rate); - std::vector clusterHistogram(numClusters+1, 0); + std::vector clusterHistogram(numClusters + 1, 0); for (std::vector::const_iterator it = timesteps.begin(); it != timesteps.end(); ++it) { unsigned cluster = getCluster(*it, minTimestep, rate); ++clusterHistogram[cluster]; } - + double numberOfUpdates = 0; double clusterTimestep = minTimestep; unsigned cluster = 0; - for (std::vector::const_iterator it = clusterHistogram.begin(); it != clusterHistogram.end(); ++it) { - std::cout << "Cluster " << cluster++ << " updates (" << *it << " elements): " << std::ceil(finalTime / clusterTimestep) << std::endl; + for (std::vector::const_iterator it = clusterHistogram.begin(); + it != clusterHistogram.end(); + ++it) { + std::cout << "Cluster " << cluster++ << " updates (" << *it + << " elements): " << std::ceil(finalTime / clusterTimestep) << std::endl; numberOfUpdates += (*it) * std::ceil(finalTime / clusterTimestep); clusterTimestep *= rate; } - + std::cout << "Minimum insphere radius: " << minInsphereRadius << std::endl; std::cout << "Elements: " << numberOfElements << std::endl; std::cout << "Minimum timestep: " << minTimestep << std::endl; std::cout << "Maximum timestep: " << maxTimestep << std::endl; std::cout << "Elements in time clusters: "; - for (std::vector::const_iterator it = clusterHistogram.begin(); it != clusterHistogram.end(); ++it) { + for (std::vector::const_iterator it = clusterHistogram.begin(); + it != clusterHistogram.end(); + ++it) { std::cout << *it << " "; } std::cout << std::endl; std::cout << "Number of cell updates (@ " << finalTime << "s): " << numberOfUpdates << std::endl; - std::cout << "Estimated time on one node: " << numberOfUpdates * flopsPerCell / (nodePerformance * 1.0e9) << std::endl; - + std::cout << "Estimated time on one node: " + << numberOfUpdates * flopsPerCell / (nodePerformance * 1.0e9) << std::endl; + return 0; } diff --git a/evaluate_material/CMakeLists.txt b/evaluate_material/CMakeLists.txt index 1fbda51..105e45f 100644 --- a/evaluate_material/CMakeLists.txt +++ b/evaluate_material/CMakeLists.txt @@ -25,5 +25,3 @@ target_link_libraries(evaluate_easi PUBLIC ${HDF5_C_HL_LIBRARIES} ${HDF5_C_LIBRA find_package(easi 1.0.0 REQUIRED) target_link_libraries(evaluate_easi PUBLIC easi::easi) - - diff --git a/evaluate_material/README.md b/evaluate_material/README.md index b474b20..70116a0 100644 --- a/evaluate_material/README.md +++ b/evaluate_material/README.md @@ -2,7 +2,7 @@ ## Compilation -``` +```bash mkdir build && cd build cmake .. make @@ -10,7 +10,7 @@ make ## Usage -``` +```bash ./evaluate_easi -m mesh.h5 -e easi.yaml -o material ``` diff --git a/evaluate_material/src/main.cpp b/evaluate_material/src/main.cpp index b9ab80e..be79504 100644 --- a/evaluate_material/src/main.cpp +++ b/evaluate_material/src/main.cpp @@ -1,14 +1,13 @@ -#include -#include - -#include - #include "mesh.h" #include "parameterDB.h" #include "utils/args.h" #include "utils/logger.h" #include "writer.h" +#include +#include +#include + int main(int argc, char** argv) { MPI_Init(&argc, &argv); diff --git a/evaluate_material/src/mesh.cpp b/evaluate_material/src/mesh.cpp index 49d2f17..e504e40 100644 --- a/evaluate_material/src/mesh.cpp +++ b/evaluate_material/src/mesh.cpp @@ -21,7 +21,8 @@ Mesh::Mesh(const std::string& fileName) { for (unsigned int i = 0; i < nElements; i++) { const auto elementVertex = puml.originalCells()[i]; - std::array tmp({elementVertex[0], elementVertex[1], elementVertex[2], elementVertex[3]}); + std::array tmp( + {elementVertex[0], elementVertex[1], elementVertex[2], elementVertex[3]}); elementVertices.emplace_back(tmp); } diff --git a/evaluate_material/src/mesh.h b/evaluate_material/src/mesh.h index 8193355..2cb5365 100644 --- a/evaluate_material/src/mesh.h +++ b/evaluate_material/src/mesh.h @@ -7,11 +7,15 @@ class Mesh { public: - explicit Mesh(std::string const& fileName); + explicit Mesh(const std::string& fileName); ~Mesh() = default; [[nodiscard]] std::vector> getElementBarycenters() const; - [[nodiscard]] const std::vector>& getVertexCoordinates() const { return vertexCoordinates; }; - [[nodiscard]] const std::vector>& getElementVertices() const { return elementVertices; }; + [[nodiscard]] const std::vector>& getVertexCoordinates() const { + return vertexCoordinates; + }; + [[nodiscard]] const std::vector>& getElementVertices() const { + return elementVertices; + }; [[nodiscard]] const std::vector& getElementGroups() const { return elementGroups; }; private: diff --git a/evaluate_material/src/parameterDB.cpp b/evaluate_material/src/parameterDB.cpp index 5d6ae2b..d7e20ea 100644 --- a/evaluate_material/src/parameterDB.cpp +++ b/evaluate_material/src/parameterDB.cpp @@ -25,14 +25,15 @@ easi::Query ParameterDB::generateQuery(const std::vector return query; } -std::pair, std::vector>> ParameterDB::evaluate(easi::Query& query) const { +std::pair, std::vector>> + ParameterDB::evaluate(easi::Query& query) const { auto supplied = model->suppliedParameters(); auto parameters = std::vector(supplied.begin(), supplied.end()); auto adapter = easi::ArraysAdapter<>{}; auto material = std::vector>(parameters.size()); auto it = material.begin(); - for (auto const& p : parameters) { + for (const auto& p : parameters) { it->resize(query.numPoints()); adapter.addBindingPoint(p, it->data()); ++it; diff --git a/evaluate_material/src/parameterDB.h b/evaluate_material/src/parameterDB.h index 63e60a6..751f826 100644 --- a/evaluate_material/src/parameterDB.h +++ b/evaluate_material/src/parameterDB.h @@ -1,10 +1,9 @@ #ifndef PARAMETER_DB_H_ #define PARAMETER_DB_H_ -#include - #include #include +#include class ParameterDB { public: @@ -12,7 +11,8 @@ class ParameterDB { ~ParameterDB(); [[nodiscard]] easi::Query generateQuery(const std::vector>& points, const std::vector& groups) const; - std::pair, std::vector>> evaluate(easi::Query& query) const; + std::pair, std::vector>> + evaluate(easi::Query& query) const; private: easi::Component* model; diff --git a/evaluate_material/src/writer.cpp b/evaluate_material/src/writer.cpp index acdc704..efd841a 100644 --- a/evaluate_material/src/writer.cpp +++ b/evaluate_material/src/writer.cpp @@ -1,20 +1,26 @@ #include "writer.h" -void Writer::writeData( - hid_t h5file, const std::string& name, const std::vector& sizes, Hdf5DataType type, const void* data) { +void Writer::writeData(hid_t h5file, + const std::string& name, + const std::vector& sizes, + Hdf5DataType type, + const void* data) { hid_t h5space = H5Screate_simple(sizes.size(), sizes.data(), nullptr); checkH5Err(h5space); hid_t h5data = 0; switch (type) { case Hdf5DataType::hdf5UInt64: - h5data = H5Dcreate(h5file, name.c_str(), H5T_STD_U64LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + h5data = H5Dcreate( + h5file, name.c_str(), H5T_STD_U64LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); break; case Hdf5DataType::hdf5Int32: - h5data = H5Dcreate(h5file, name.c_str(), H5T_STD_I32LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + h5data = H5Dcreate( + h5file, name.c_str(), H5T_STD_I32LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); break; case Hdf5DataType::hdf5Double: - h5data = H5Dcreate(h5file, name.c_str(), H5T_IEEE_F64LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + h5data = H5Dcreate( + h5file, name.c_str(), H5T_IEEE_F64LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); break; } checkH5Err(h5data); @@ -57,21 +63,24 @@ void Writer::writeXdmf(const std::vector& parameterNames) { << std::endl // This should be UInt but for some reason this does not work with // binary data - << R"( )" - << hdfFileName << R"(:/connect)" << std::endl + << R"( )" << hdfFileName << R"(:/connect)" << std::endl << R"( )" << std::endl - << R"( )" << std::endl - << R"( )" << hdfFileName << R"(:/geometry)" << std::endl + << R"( )" + << std::endl + << R"( )" << hdfFileName + << R"(:/geometry)" << std::endl << R"( )" << std::endl << R"( )" << std::endl - << R"( )" - << hdfFileName << R"(:/group)" << std::endl + << R"( )" << hdfFileName << R"(:/group)" << std::endl << R"( )" << std::endl; for (const auto& param : parameterNames) { xdmf << R"( )" << std::endl - << R"( )" << hdfFileName << R"(:/)" << param << R"()" << std::endl + << R"( )" << hdfFileName << R"(:/)" + << param << R"()" << std::endl << R"( )" << std::endl; } xdmf << " " << std::endl << " " << std::endl << "" << std::endl; diff --git a/evaluate_material/src/writer.h b/evaluate_material/src/writer.h index 4e0a5d2..428c11f 100644 --- a/evaluate_material/src/writer.h +++ b/evaluate_material/src/writer.h @@ -1,13 +1,12 @@ #ifndef WRITER_H_ #define WRITER_H_ +#include "mesh.h" + #include #include -#include - #include - -#include "mesh.h" +#include enum class Hdf5DataType { hdf5UInt64, hdf5Int32, hdf5Double }; @@ -22,9 +21,10 @@ static void checkH5ErrImpl(TT status, const char* file, int line) { class Writer { public: explicit Writer(const std::string& fileNamePrefix, const Mesh& mesh) - : mesh(mesh), xdmfFileName(fileNamePrefix + ".xdmf"), hdfFileName(fileNamePrefix + ".h5"){}; + : mesh(mesh), xdmfFileName(fileNamePrefix + ".xdmf"), hdfFileName(fileNamePrefix + ".h5") {}; - void write(const std::vector& parameterNames, const std::vector>& materialValues) { + void write(const std::vector& parameterNames, + const std::vector>& materialValues) { assert(parameterNames.size() == materialValues.size()); writeXdmf(parameterNames); writeHdf5(parameterNames, materialValues); @@ -38,7 +38,10 @@ class Writer { void writeHdf5(const std::vector& parameterNames, const std::vector>& materialValues); void writeXdmf(const std::vector& parameterNames); - void - writeData(hid_t h5file, const std::string& name, const std::vector& sizes, Hdf5DataType type, const void* data); + void writeData(hid_t h5file, + const std::string& name, + const std::vector& sizes, + Hdf5DataType type, + const void* data); }; #endif diff --git a/extract_timesteps/SConstruct b/extract_timesteps/SConstruct index ee47ac3..adcbdbe 100644 --- a/extract_timesteps/SConstruct +++ b/extract_timesteps/SConstruct @@ -1,4 +1,3 @@ -#! /usr/bin/python ## # @file # This file is part of SeisSol. @@ -8,21 +7,21 @@ # @section LICENSE # Copyright (c) 2017, SeisSol Group # All rights reserved. -# +# # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: -# +# # 1. Redistributions of source code must retain the above copyright notice, # this list of conditions and the following disclaimer. -# +# # 2. Redistributions in binary form must reproduce the above copyright notice, # this list of conditions and the following disclaimer in the documentation # and/or other materials provided with the distribution. -# +# # 3. Neither the name of the copyright holder nor the names of its # contributors may be used to endorse or promote products derived from this # software without specific prior written permission. -# +# # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE diff --git a/extract_timesteps/src/SConscript b/extract_timesteps/src/SConscript index 797de3a..5419d09 100644 --- a/extract_timesteps/src/SConscript +++ b/extract_timesteps/src/SConscript @@ -1,4 +1,3 @@ -#!/usr/bin/python ## # @file # This file is part of SeisSol. @@ -8,21 +7,21 @@ # @section LICENSE # Copyright (c) 2017, SeisSol Group # All rights reserved. -# +# # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: -# +# # 1. Redistributions of source code must retain the above copyright notice, # this list of conditions and the following disclaimer. -# +# # 2. Redistributions in binary form must reproduce the above copyright notice, # this list of conditions and the following disclaimer in the documentation # and/or other materials provided with the distribution. -# +# # 3. Neither the name of the copyright holder nor the names of its # contributors may be used to endorse or promote products derived from this # software without specific prior written permission. -# +# # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE diff --git a/extract_timesteps/src/SeismicVelocity.cpp b/extract_timesteps/src/SeismicVelocity.cpp index e1558b7..c1f1b30 100644 --- a/extract_timesteps/src/SeismicVelocity.cpp +++ b/extract_timesteps/src/SeismicVelocity.cpp @@ -1,7 +1,6 @@ #include "SeismicVelocity.h" -double landers61(int, double, double, double z) -{ +double landers61(int, double, double, double z) { if (z > -0.1) { return 1925.9102873142; } @@ -26,113 +25,111 @@ double landers61(int, double, double, double z) } /* -BedrockVelModel(1,:) = (/ -6d3, 2550d0,18589500000d0,26571000000d0/) => 5000 m/s -BedrockVelModel(2,:) = (/ -8d3, 2850d0,39016500000d0,42379500000d0/) => 6500 m/s -BedrockVelModel(3,:) = (/ -12d3, 3050d0,50027625000d0,53695250000d0/) => 7100 m/s -BedrockVelModel(4,:) = (/-6d3,2720d0,33320000000d0,31280000000d0/) => 6000 m/s -BedrockVelModel(5,:) = (/-12d3,2860d0,41298400000d0,41984800000d0/) => 6600 m/s -BedrockVelModel(6,:) = (/-23d3,3050d0,46390500000d0,60969500000d0/) => 7100 m/s -BedrockVelModel(7,:) = (/ -5d10, 3330d0,65942325000d0,81235350000d0/) => 8000 m/s */ -double sumatra1223_high(int group, double, double, double z) -{ +BedrockVelModel(1,:) = (/ -6d3, 2550d0,18589500000d0,26571000000d0/) => 5000 +m/s BedrockVelModel(2,:) = (/ -8d3, 2850d0,39016500000d0,42379500000d0/) => +6500 m/s BedrockVelModel(3,:) = (/ -12d3, 3050d0,50027625000d0,53695250000d0/) +=> 7100 m/s BedrockVelModel(4,:) = (/-6d3,2720d0,33320000000d0,31280000000d0/) +=> 6000 m/s BedrockVelModel(5,:) = (/-12d3,2860d0,41298400000d0,41984800000d0/) +=> 6600 m/s BedrockVelModel(6,:) = (/-23d3,3050d0,46390500000d0,60969500000d0/) +=> 7100 m/s BedrockVelModel(7,:) = (/ -5d10, +3330d0,65942325000d0,81235350000d0/) => 8000 m/s */ +double sumatra1223_high(int group, double, double, double z) { double cp = -1.0; switch (group) { - case 4: // LVZ - case 2: // layer 1 - cp = 5000.0; - break; - case 3: // layer 2 - cp = 6500.0; - break; - case 7: // layer 3 + case 4: // LVZ + case 2: // layer 1 + cp = 5000.0; + break; + case 3: // layer 2 + cp = 6500.0; + break; + case 7: // layer 3 + cp = 7100.0; + break; + case 6: // layer 4 + cp = 8000.0; + break; + case 1: // continental + case 5: // box + if (z > -6000.0) { + cp = 6000.0; + } else if (z >= -12000.0) { + cp = 6600.0; + } else if (z > -23000.0) { cp = 7100.0; - break; - case 6: // layer 4 + } else { cp = 8000.0; - break; - case 1: // continental - case 5: // box - if (z > -6000.0) { - cp = 6000.0; - } else if (z >= -12000.0) { - cp = 6600.0; - } else if (z > -23000.0) { - cp = 7100.0; - } else { - cp = 8000.0; - } - break; - default: - break; + } + break; + default: + break; } return cp; } -double sumatra1223_low(int group, double, double, double z) -{ +double sumatra1223_low(int group, double, double, double z) { double cp = -1.0; switch (group) { - case 2: // LVZ - case 7: // layer 1 - cp = 5000.0; - break; - case 4: // layer 2 - cp = 6500.0; - break; - case 3: // layer 3 + case 2: // LVZ + case 7: // layer 1 + cp = 5000.0; + break; + case 4: // layer 2 + cp = 6500.0; + break; + case 3: // layer 3 + cp = 7100.0; + break; + case 6: // layer 4 + cp = 8000.0; + break; + case 1: // continental + case 5: // box + if (z > -6000.0) { + cp = 6000.0; + } else if (z >= -12000.0) { + cp = 6600.0; + } else if (z > -23000.0) { cp = 7100.0; - break; - case 6: // layer 4 + } else { cp = 8000.0; - break; - case 1: // continental - case 5: // box - if (z > -6000.0) { - cp = 6000.0; - } else if (z >= -12000.0) { - cp = 6600.0; - } else if (z > -23000.0) { - cp = 7100.0; - } else { - cp = 8000.0; - } - break; - default: - break; + } + break; + default: + break; } return cp; } -double sumatra1224(int group, double, double, double z) -{ +double sumatra1224(int group, double, double, double z) { double cp = -1.0; switch (group) { - case 1: - cp = 6500.0; - break; - case 2: - cp = 5000.0; - break; - case 3: + case 1: + cp = 6500.0; + break; + case 2: + cp = 5000.0; + break; + case 3: + cp = 7100.0; + break; + case 4: + case 5: + cp = 8000.0; + break; + case 6: + if (z > -6000.0) { + cp = 6000.0; + } else if (z >= -12000.0) { + cp = 6600.0; + } else if (z > -23000.0) { cp = 7100.0; - break; - case 4: - case 5: + } else { cp = 8000.0; - break; - case 6: - if (z > -6000.0) { - cp = 6000.0; - } else if (z >= -12000.0) { - cp = 6600.0; - } else if (z > -23000.0) { - cp = 7100.0; - } else { - cp = 8000.0; - } - break; - default: - break; + } + break; + default: + break; } return cp; } diff --git a/extract_timesteps/src/Writer.cpp b/extract_timesteps/src/Writer.cpp index be79797..b127756 100644 --- a/extract_timesteps/src/Writer.cpp +++ b/extract_timesteps/src/Writer.cpp @@ -1,21 +1,20 @@ #include "Writer.h" + +#include #include #include -#include -static void check_err(int const stat, int const line, char const* file) -{ +static void check_err(const int stat, const int line, const char* file) { if (stat != NC_NOERR) { - fprintf(stderr,"line %d of %s: %s\n", line, file, nc_strerror(stat)); + fprintf(stderr, "line %d of %s: %s\n", line, file, nc_strerror(stat)); fflush(stderr); exit(1); } } -void writeTimesteps(std::vector const& timesteps, std::string const& fileName) -{ - int stat; /* return status */ - int ncid; /* netCDF id */ +void writeTimesteps(const std::vector& timesteps, const std::string& fileName) { + int stat; /* return status */ + int ncid; /* netCDF id */ /* dimension ids */ int element_dim; @@ -27,38 +26,37 @@ void writeTimesteps(std::vector const& timesteps, std::string const& fil int timestep_id; /* rank (number of dimensions) for each variable */ -# define RANK_timestep 1 +#define RANK_timestep 1 /* variable shapes */ int timestep_dims[RANK_timestep]; /* enter define mode */ stat = nc_create(fileName.c_str(), NC_CLOBBER, &ncid); - check_err(stat,__LINE__,__FILE__); + check_err(stat, __LINE__, __FILE__); /* define dimensions */ stat = nc_def_dim(ncid, "element", element_len, &element_dim); - check_err(stat,__LINE__,__FILE__); + check_err(stat, __LINE__, __FILE__); /* define variables */ timestep_dims[0] = element_dim; stat = nc_def_var(ncid, "timestep", NC_DOUBLE, RANK_timestep, timestep_dims, ×tep_id); - check_err(stat,__LINE__,__FILE__); + check_err(stat, __LINE__, __FILE__); /* assign per-variable attributes */ stat = nc_put_att_text(ncid, timestep_id, "units", 1, "s"); - check_err(stat,__LINE__,__FILE__); - + check_err(stat, __LINE__, __FILE__); /* leave define mode */ - stat = nc_enddef (ncid); - check_err(stat,__LINE__,__FILE__); + stat = nc_enddef(ncid); + check_err(stat, __LINE__, __FILE__); /* assign variable data */ stat = nc_put_var_double(ncid, timestep_id, ×teps.data()[0]); - check_err(stat,__LINE__,__FILE__); + check_err(stat, __LINE__, __FILE__); stat = nc_close(ncid); - check_err(stat,__LINE__,__FILE__); + check_err(stat, __LINE__, __FILE__); } diff --git a/extract_timesteps/src/Writer.h b/extract_timesteps/src/Writer.h index 3a98abf..329419c 100644 --- a/extract_timesteps/src/Writer.h +++ b/extract_timesteps/src/Writer.h @@ -1,9 +1,9 @@ #ifndef WRITER_H_ #define WRITER_H_ -#include #include +#include -void writeTimesteps(std::vector const& timesteps, std::string const& fileName); +void writeTimesteps(const std::vector& timesteps, const std::string& fileName); #endif diff --git a/extract_timesteps/src/main.cpp b/extract_timesteps/src/main.cpp index 33e5d6d..152edcd 100644 --- a/extract_timesteps/src/main.cpp +++ b/extract_timesteps/src/main.cpp @@ -1,19 +1,17 @@ -#include -#include -#include -#include - #include "SeismicVelocity.h" #include "Writer.h" #include "common/PartitionReader.h" -double calculateTimestep(double cfl, double insphereRadius, double pVelocity, unsigned order) -{ - return cfl * 2.0 * insphereRadius / (pVelocity) / (2*(order-1)+1); +#include +#include +#include +#include + +double calculateTimestep(double cfl, double insphereRadius, double pVelocity, unsigned order) { + return cfl * 2.0 * insphereRadius / (pVelocity) / (2 * (order - 1) + 1); } -int main(int argc, char** argv) -{ +int main(int argc, char** argv) { utils::Args args; args.addOption("mesh", 'm', "Netcdf mesh file"); args.addOption("velocity-model", 'v', "Velocity model"); @@ -21,17 +19,17 @@ int main(int argc, char** argv) args.addOption("order", 'O', "Convergence order"); args.addOption("output", 'o', "Output file name"); - if (args.parse(argc, argv) != utils::Args::Success) { + if (args.parse(argc, argv) != utils::Args::Success) { return -1; } - + std::string meshFile = args.getArgument("mesh"); std::string velocityModel = args.getArgument("velocity-model"); double cfl = args.getArgument("CFL", 0.5); unsigned order = args.getArgument("order"); std::string outputFile = args.getArgument("output"); - - double (*pWaveVelocityFun)(int,double,double,double); + + double (*pWaveVelocityFun)(int, double, double, double); if (velocityModel == "landers61") { pWaveVelocityFun = &landers61; } else if (velocityModel == "sumatra1223_low") { @@ -44,40 +42,47 @@ int main(int argc, char** argv) std::cerr << "Error: Unknown velocity model." << std::endl; exit(-1); } - + PartitionReader reader(meshFile); std::vector timesteps; unsigned numberOfElements = 0; for (unsigned p = 0; p < reader.partitions; ++p) { - std::cout << "Reading partition " << p+1 << "..." << std::endl; + std::cout << "Reading partition " << p + 1 << "..." << std::endl; reader.readPartition(p); numberOfElements += reader.elementSize[p]; for (unsigned elem = 0; elem < reader.elementSize[p]; ++elem) { - glm::dvec3 barycentre(0.,0.,0.); + glm::dvec3 barycentre(0., 0., 0.); glm::dvec3 x[4]; for (unsigned vtx = 0; vtx < 4; ++vtx) { for (unsigned d = 0; d < 3; ++d) { - x[vtx][d] = reader.vertexCoordinates[3 * reader.elementVertices[4*elem + vtx] + d]; + x[vtx][d] = reader.vertexCoordinates[3 * reader.elementVertices[4 * elem + vtx] + d]; } barycentre += 0.25 * x[vtx]; } - double alpha = determinant(glm::dmat4(glm::dvec4(x[0], 1.0), glm::dvec4(x[1], 1.0), glm::dvec4(x[2], 1.0), glm::dvec4(x[3], 1.0))); - double Nabc = length(cross(x[1]-x[0], x[2]-x[0])); - double Nabd = length(cross(x[1]-x[0], x[3]-x[0])); - double Nacd = length(cross(x[2]-x[0], x[3]-x[0])); - double Nbcd = length(cross(x[2]-x[1], x[3]-x[1])); + double alpha = determinant(glm::dmat4(glm::dvec4(x[0], 1.0), + glm::dvec4(x[1], 1.0), + glm::dvec4(x[2], 1.0), + glm::dvec4(x[3], 1.0))); + double Nabc = length(cross(x[1] - x[0], x[2] - x[0])); + double Nabd = length(cross(x[1] - x[0], x[3] - x[0])); + double Nacd = length(cross(x[2] - x[0], x[3] - x[0])); + double Nbcd = length(cross(x[2] - x[1], x[3] - x[1])); double insphere = std::fabs(alpha) / (Nabc + Nabd + Nacd + Nbcd); - - double timestep = calculateTimestep(cfl, insphere, pWaveVelocityFun(reader.elementGroup[elem], barycentre.x, barycentre.y, barycentre.z), order); - + + double timestep = calculateTimestep( + cfl, + insphere, + pWaveVelocityFun(reader.elementGroup[elem], barycentre.x, barycentre.y, barycentre.z), + order); + timesteps.push_back(timestep); } } - + writeTimesteps(timesteps, outputFile); - + return 0; } diff --git a/mirrorMesh/mirrorMesh.py b/mirrorMesh/mirrorMesh.py old mode 100644 new mode 100755 index ff58240..334b9cb --- a/mirrorMesh/mirrorMesh.py +++ b/mirrorMesh/mirrorMesh.py @@ -83,9 +83,9 @@ def write_xdmfh5(fn, xyzb, tetrab, groupb, boundaryb): xyzb[0:nNodes,:] = xyz for inew in range(nNodes, Nnodesb): i = invnodesLU[inew] - xyzb[inew,:] = xyzb[i,:] + xyzb[inew,:] = xyzb[i,:] xyzb[inew,iN] = xc-(xyzb[i,iN]-xc) - + #create new connect tetrab = np.zeros((nElements*2,4), dtype=int) @@ -130,4 +130,3 @@ def write_xdmfh5(fn, xyzb, tetrab, groupb, boundaryb): boundaryb = boundaryFaceb[:,0]+256*boundaryFaceb[:,1]+256*256*boundaryFaceb[:,2]+256*256*256*boundaryFaceb[:,3] write_xdmfh5(fn, xyzb, tetrab, groupb, boundaryb) - diff --git a/netcdf2xdmf/SConstruct b/netcdf2xdmf/SConstruct index 0243be7..35ed76f 100644 --- a/netcdf2xdmf/SConstruct +++ b/netcdf2xdmf/SConstruct @@ -1,4 +1,3 @@ -#!/usr/bin/python ## # @file # This file is part of SeisSol. @@ -65,7 +64,7 @@ vars.SetCompiler(env) libs.find(env, 'netcdf') hasAsagi = libs.find(env, 'asagi', required=False, parallel=False) if hasAsagi: - env.Append(CPPDEFINES=['USE_ASAGI']) + env.Append(CPPDEFINES=['USE_ASAGI']) env.Tool('cmake') diff --git a/netcdf2xdmf/src/SConscript b/netcdf2xdmf/src/SConscript index 0a1a958..dce56aa 100644 --- a/netcdf2xdmf/src/SConscript +++ b/netcdf2xdmf/src/SConscript @@ -1,4 +1,3 @@ -#!/usr/bin/python ## # @file # This file is part of SeisSol. @@ -44,4 +43,4 @@ sourceFiles = ['main.cpp'] for f in sourceFiles: env.sourceFiles.append(env.Object(f)) -Export('env') \ No newline at end of file +Export('env') diff --git a/netcdf2xdmf/src/main.cpp b/netcdf2xdmf/src/main.cpp index b0acd52..c73b520 100644 --- a/netcdf2xdmf/src/main.cpp +++ b/netcdf2xdmf/src/main.cpp @@ -2,7 +2,8 @@ * @file * This file is part of SeisSol. * - * @author Sebastian Rettenberger (sebastian.rettenberger AT tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger) + * @author Sebastian Rettenberger (sebastian.rettenberger AT tum.de, + * http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger) * * @section LICENSE * Copyright (c) 2016, SeisSol Group @@ -37,369 +38,388 @@ * @section DESCRIPTION */ -#include -#include -#include -#include +#include "easi/YAMLParser.h" +#include "utils/args.h" +#include "utils/logger.h" +#include "utils/path.h" +#include "utils/stringutils.h" #include +#include #include +#include #include +#include +#include +#include #include -#include +int main(int argc, char* argv[]) { + utils::Args args; + args.addOption("boundary", 'b', "Convert only boundary cells", utils::Args::No, false); + args.addOption( + "materialFile", 'm', "Add material parameters to xdmf", utils::Args::Required, false); + args.addAdditionalOption("input", "The netCDF mesh file"); + args.addAdditionalOption("output", "The generated XDMF file", false); + + switch (args.parse(argc, argv)) { + case utils::Args::Error: + return 1; + case utils::Args::Help: + return 2; + } + + bool convertBoundaries = args.isSet("boundary"); + std::string material = args.getArgument("materialFile", ""); + + // Get infile/outfile + std::string infile = args.getAdditionalArgument("input"); + std::string outfilePrefix; + if (args.isSetAdditional("output")) { + outfilePrefix = args.getAdditionalArgument("output"); + } else { + outfilePrefix = infile; + size_t pos = outfilePrefix.find_last_of('.'); + outfilePrefix.resize(pos); + } + + // Remove file ending + if (utils::StringUtils::endsWith(outfilePrefix, ".xdmf")) { + outfilePrefix.resize(outfilePrefix.size() - 5); + } else { + if (convertBoundaries) { + if (!utils::StringUtils::endsWith(outfilePrefix, "_bnd")) + outfilePrefix += "_bnd"; + } + } + + std::string outfilePrefixShort = utils::Path(outfilePrefix).basename(); + + // Open the xml file + std::ofstream xdmfFile((outfilePrefix + ".xdmf").c_str()); + + xdmfFile << "" << std::endl + << "" << std::endl + << "" << std::endl + << " " << std::endl; + xdmfFile << " " << std::endl; + + // + // m_xdmfFile << " " + // << + // backends::Backend::dataItemLocation(backendPrefix.c_str(), "partition") + // << "" << std::endl; + + // Open the netCDF file + int ncFile; + nc_open(infile.c_str(), 0, &ncFile); + + int ncDimPart; + nc_inq_dimid(ncFile, "partitions", &ncDimPart); + size_t partitions; + nc_inq_dimlen(ncFile, ncDimPart, &partitions); + + int ncDimVrtx; + nc_inq_dimid(ncFile, "vertices", &ncDimVrtx); + size_t maxVertices; + nc_inq_dimlen(ncFile, ncDimVrtx, &maxVertices); + + int ncDimElem; + nc_inq_dimid(ncFile, "elements", &ncDimElem); + size_t maxElements; + nc_inq_dimlen(ncFile, ncDimElem, &maxElements); + + // Create netcdf variables + int ncVarElemSize; + nc_inq_varid(ncFile, "element_size", &ncVarElemSize); + + int ncVarElemVertices; + nc_inq_varid(ncFile, "element_vertices", &ncVarElemVertices); + + int ncVarElemBoundaries; + nc_inq_varid(ncFile, "element_boundaries", &ncVarElemBoundaries); + + int ncVarElemGroup; + bool hasGroup = false; + int ncResult = nc_inq_varid(ncFile, "element_group", &ncVarElemGroup); + if (ncResult != NC_ENOTVAR) + hasGroup = true; + + int ncVarVrtxSize; + nc_inq_varid(ncFile, "vertex_size", &ncVarVrtxSize); + + int ncVarVrtxCoords; + nc_inq_varid(ncFile, "vertex_coordinates", &ncVarVrtxCoords); + + unsigned int* vertexSize = new unsigned int[partitions]; + nc_get_var_uint(ncFile, ncVarVrtxSize, vertexSize); + unsigned int* elementSize = new unsigned int[partitions]; + nc_get_var_uint(ncFile, ncVarElemSize, elementSize); + + // Write cells + std::string topoName; + unsigned int topoSize; + + unsigned int nvertices = std::valarray(vertexSize, partitions).sum(); + unsigned int nelements = 0; + if (convertBoundaries) { + topoName = "Triangle"; + topoSize = 3; + + unsigned int* boundaries = new unsigned int[maxElements * 4]; + + size_t start[3] = {0, 0, 0}; + size_t size[3] = {1, 0, 4}; + for (size_t i = 0; i < partitions; i++) { + start[0] = i; + size[1] = elementSize[i]; + nc_get_vara_uint(ncFile, ncVarElemBoundaries, start, size, boundaries); + + for (unsigned int j = 0; j < elementSize[i]; j++) { + for (unsigned int k = 0; k < 4; k++) { + if (boundaries[j * 4 + k] != 0) + nelements++; + } + } + } -#include "utils/args.h" -#include "utils/logger.h" -#include "utils/path.h" -#include "utils/stringutils.h" + delete[] boundaries; + } else { + topoName = "Tetrahedron"; + topoSize = 4; + + nelements = std::valarray(elementSize, partitions).sum(); + } + + xdmfFile << " " + << std::endl + // This should be UInt but for some reason this does not work with + // binary data + << " " << outfilePrefixShort + << "_connect.bin" << std::endl + << " " << std::endl; + + int fd = open((outfilePrefix + "_connect.bin").c_str(), + O_CREAT | O_WRONLY | O_TRUNC, + S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH); + + int* elements = new int[maxElements * 4]; + unsigned int vrtxIdStart = 0; + if (convertBoundaries) { + int* connect = new int[nelements * 3]; + unsigned int nextElement = 0; + + unsigned int* boundaries = new unsigned int[maxElements * 4]; + + for (size_t i = 0; i < partitions; i++) { + size_t start[3] = {i, 0, 0}; + size_t size[3] = {1, elementSize[i], 4}; + nc_get_vara_int(ncFile, ncVarElemVertices, start, size, elements); + nc_get_vara_uint(ncFile, ncVarElemBoundaries, start, size, boundaries); + + for (unsigned int j = 0; j < elementSize[i]; j++) { + if (boundaries[j * 4] != 0) { + connect[nextElement * 3] = vrtxIdStart + elements[j * 4]; + connect[nextElement * 3 + 1] = vrtxIdStart + elements[j * 4 + 1]; + connect[nextElement * 3 + 2] = vrtxIdStart + elements[j * 4 + 2]; + nextElement++; + } + if (boundaries[j * 4 + 1] != 0) { + connect[nextElement * 3] = vrtxIdStart + elements[j * 4]; + connect[nextElement * 3 + 1] = vrtxIdStart + elements[j * 4 + 1]; + connect[nextElement * 3 + 2] = vrtxIdStart + elements[j * 4 + 3]; + nextElement++; + } + if (boundaries[j * 4 + 2] != 0) { + connect[nextElement * 3] = vrtxIdStart + elements[j * 4]; + connect[nextElement * 3 + 1] = vrtxIdStart + elements[j * 4 + 2]; + connect[nextElement * 3 + 2] = vrtxIdStart + elements[j * 4 + 3]; + nextElement++; + } + if (boundaries[j * 4 + 3] != 0) { + connect[nextElement * 3] = vrtxIdStart + elements[j * 4 + 1]; + connect[nextElement * 3 + 1] = vrtxIdStart + elements[j * 4 + 2]; + connect[nextElement * 3 + 2] = vrtxIdStart + elements[j * 4 + 3]; + nextElement++; + } + } -#include "easi/YAMLParser.h" + vrtxIdStart += vertexSize[i]; + } + + delete[] boundaries; + + write(fd, connect, nelements * 3 * sizeof(int)); + } else { + for (size_t i = 0; i < partitions; i++) { + size_t start[3] = {i, 0, 0}; + size_t size[3] = {1, elementSize[i], 4}; + nc_get_vara_int(ncFile, ncVarElemVertices, start, size, elements); + + for (unsigned int j = 0; j < elementSize[i] * 4; j++) + elements[j] = vrtxIdStart + elements[j]; + + vrtxIdStart += vertexSize[i]; -int main(int argc, char* argv[]) -{ - utils::Args args; - args.addOption("boundary", 'b', "Convert only boundary cells", utils::Args::No, false); - args.addOption("materialFile", 'm', "Add material parameters to xdmf", utils::Args::Required, false); - args.addAdditionalOption("input", "The netCDF mesh file"); - args.addAdditionalOption("output", "The generated XDMF file", false); - - switch (args.parse(argc, argv)) { - case utils::Args::Error: - return 1; - case utils::Args::Help: - return 2; - } - - bool convertBoundaries = args.isSet("boundary"); - std::string material = args.getArgument("materialFile", ""); - - // Get infile/outfile - std::string infile = args.getAdditionalArgument("input"); - std::string outfilePrefix; - if (args.isSetAdditional("output")) { - outfilePrefix = args.getAdditionalArgument("output"); - } else { - outfilePrefix = infile; - size_t pos = outfilePrefix.find_last_of('.'); - outfilePrefix.resize(pos); - } - - // Remove file ending - if (utils::StringUtils::endsWith(outfilePrefix, ".xdmf")) { - outfilePrefix.resize(outfilePrefix.size()-5); - } else { - if (convertBoundaries) { - if (!utils::StringUtils::endsWith(outfilePrefix, "_bnd")) - outfilePrefix += "_bnd"; - } - } - - std::string outfilePrefixShort = utils::Path(outfilePrefix).basename(); - - // Open the xml file - std::ofstream xdmfFile((outfilePrefix+".xdmf").c_str()); - - xdmfFile << "" << std::endl - << "" << std::endl - << "" << std::endl - << " " << std::endl; - xdmfFile << " " << std::endl; - -// -// m_xdmfFile << " " -// << backends::Backend::dataItemLocation(backendPrefix.c_str(), "partition") -// << "" << std::endl; - - // Open the netCDF file - int ncFile; - nc_open(infile.c_str(), 0, &ncFile); - - int ncDimPart; - nc_inq_dimid(ncFile, "partitions", &ncDimPart); - size_t partitions; - nc_inq_dimlen(ncFile, ncDimPart, &partitions); - - int ncDimVrtx; - nc_inq_dimid(ncFile, "vertices", &ncDimVrtx); - size_t maxVertices; - nc_inq_dimlen(ncFile, ncDimVrtx, &maxVertices); - - int ncDimElem; - nc_inq_dimid(ncFile, "elements", &ncDimElem); - size_t maxElements; - nc_inq_dimlen(ncFile, ncDimElem, &maxElements); - - // Create netcdf variables - int ncVarElemSize; - nc_inq_varid(ncFile, "element_size", &ncVarElemSize); - - int ncVarElemVertices; - nc_inq_varid(ncFile, "element_vertices", &ncVarElemVertices); - - int ncVarElemBoundaries; - nc_inq_varid(ncFile, "element_boundaries", &ncVarElemBoundaries); - - int ncVarElemGroup; - bool hasGroup = false; - int ncResult = nc_inq_varid(ncFile, "element_group", &ncVarElemGroup); - if (ncResult != NC_ENOTVAR) - hasGroup = true; - - int ncVarVrtxSize; - nc_inq_varid(ncFile, "vertex_size", &ncVarVrtxSize); - - int ncVarVrtxCoords; - nc_inq_varid(ncFile, "vertex_coordinates", &ncVarVrtxCoords); - - unsigned int* vertexSize = new unsigned int[partitions]; - nc_get_var_uint(ncFile, ncVarVrtxSize, vertexSize); - unsigned int* elementSize = new unsigned int[partitions]; - nc_get_var_uint(ncFile, ncVarElemSize, elementSize); - - // Write cells - std::string topoName; - unsigned int topoSize; - - unsigned int nvertices = std::valarray(vertexSize, partitions).sum(); - unsigned int nelements = 0; - if (convertBoundaries) { - topoName = "Triangle"; - topoSize = 3; - - unsigned int* boundaries = new unsigned int[maxElements*4]; - - size_t start[3] = {0, 0, 0}; - size_t size[3] = {1, 0, 4}; - for (size_t i = 0; i < partitions; i++) { - start[0] = i; - size[1] = elementSize[i]; - nc_get_vara_uint(ncFile, ncVarElemBoundaries, start, size, boundaries); - - for (unsigned int j = 0; j < elementSize[i]; j++) { - for (unsigned int k = 0; k < 4; k++) { - if (boundaries[j*4+k] != 0) - nelements++; - } - } - } - - delete [] boundaries; - } else { - topoName = "Tetrahedron"; - topoSize = 4; - - nelements = std::valarray(elementSize, partitions).sum(); - } - - xdmfFile << " " << std::endl - // This should be UInt but for some reason this does not work with binary data - << " " - << outfilePrefixShort << "_connect.bin" << std::endl - << " " << std::endl; - - int fd = open((outfilePrefix+"_connect.bin").c_str(), O_CREAT | O_WRONLY | O_TRUNC, - S_IRUSR |S_IWUSR |S_IRGRP |S_IWGRP | S_IROTH | S_IWOTH); - - int* elements = new int[maxElements*4]; - unsigned int vrtxIdStart = 0; - if (convertBoundaries) { - int* connect = new int[nelements*3]; - unsigned int nextElement = 0; - - unsigned int* boundaries = new unsigned int[maxElements*4]; - - for (size_t i = 0; i < partitions; i++) { - size_t start[3] = {i, 0, 0}; - size_t size[3] = {1, elementSize[i], 4}; - nc_get_vara_int(ncFile, ncVarElemVertices, start, size, elements); - nc_get_vara_uint(ncFile, ncVarElemBoundaries, start, size, boundaries); - - for (unsigned int j = 0; j < elementSize[i]; j++) { - if (boundaries[j*4] != 0) { - connect[nextElement*3] = vrtxIdStart+elements[j*4]; - connect[nextElement*3+1] = vrtxIdStart+elements[j*4+1]; - connect[nextElement*3+2] = vrtxIdStart+elements[j*4+2]; - nextElement++; - } - if (boundaries[j*4+1] != 0) { - connect[nextElement*3] = vrtxIdStart+elements[j*4]; - connect[nextElement*3+1] = vrtxIdStart+elements[j*4+1]; - connect[nextElement*3+2] = vrtxIdStart+elements[j*4+3]; - nextElement++; - } - if (boundaries[j*4+2] != 0) { - connect[nextElement*3] = vrtxIdStart+elements[j*4]; - connect[nextElement*3+1] = vrtxIdStart+elements[j*4+2]; - connect[nextElement*3+2] = vrtxIdStart+elements[j*4+3]; - nextElement++; - } - if (boundaries[j*4+3] != 0) { - connect[nextElement*3] = vrtxIdStart+elements[j*4+1]; - connect[nextElement*3+1] = vrtxIdStart+elements[j*4+2]; - connect[nextElement*3+2] = vrtxIdStart+elements[j*4+3]; - nextElement++; - } - - } - - vrtxIdStart += vertexSize[i]; - } - - delete [] boundaries; - - write(fd, connect, nelements*3*sizeof(int)); - } else { - for (size_t i = 0; i < partitions; i++) { - size_t start[3] = {i, 0, 0}; - size_t size[3] = {1, elementSize[i], 4}; - nc_get_vara_int(ncFile, ncVarElemVertices, start, size, elements); - - for (unsigned int j = 0; j < elementSize[i]*4; j++) - elements[j] = vrtxIdStart + elements[j]; - - vrtxIdStart += vertexSize[i]; - - write(fd, elements, elementSize[i]*4*sizeof(int)); - } - } - delete [] elements; - - close(fd); - - // Write vertices - xdmfFile << " " << std::endl - << " " - << outfilePrefixShort << "_geometry.bin" << std::endl - << " " << std::endl; - - fd = open((outfilePrefix+"_geometry.bin").c_str(), O_CREAT | O_WRONLY | O_TRUNC, - S_IRUSR |S_IWUSR |S_IRGRP |S_IWGRP | S_IROTH | S_IWOTH); - - float* vertices = new float[maxVertices*3]; - for (size_t i = 0; i < partitions; i++) { - size_t start[3] = {i, 0, 0}; - size_t size[3] = {1, vertexSize[i], 3}; - nc_get_vara_float(ncFile, ncVarVrtxCoords, start, size, vertices); - - write(fd, vertices, vertexSize[i]*3*sizeof(float)); - } - - delete [] vertices; - - close(fd); - - // Write Partition/Boundary type - std::string attributeName; - if (convertBoundaries) - attributeName = "boundary"; - else - attributeName = "partition"; - - xdmfFile << " " << std::endl - << " " - << outfilePrefixShort << "_" << attributeName << ".bin" << std::endl - << " " << std::endl; - - fd = open((outfilePrefix+"_"+attributeName+".bin").c_str(), O_CREAT | O_WRONLY | O_TRUNC, - S_IRUSR |S_IWUSR |S_IRGRP |S_IWGRP | S_IROTH | S_IWOTH); - - if (convertBoundaries) { - unsigned int* boundaries = new unsigned int[maxElements*4]; - - int* boundariesOut = new int[nelements]; - unsigned int nextElement = 0; - - for (size_t i = 0; i < partitions; i++) { - size_t start[3] = {i, 0, 0}; - size_t size[3] = {1, elementSize[i], 4}; - nc_get_vara_uint(ncFile, ncVarElemBoundaries, start, size, boundaries); - - for (unsigned int j = 0; j < elementSize[i]; j++) { - for (unsigned int k = 0; k < 4; k++) { - if (boundaries[j*4+k] != 0) - boundariesOut[nextElement++] = boundaries[j*4+k]; - } - } - } - - write(fd, boundariesOut, nelements*sizeof(int)); - - delete [] boundaries; - delete [] boundariesOut; - } else { - int* partition = new int[maxElements]; - - for (size_t i = 0; i < partitions; i++) { - std::fill_n(partition, elementSize[i], i); - - write(fd, partition, elementSize[i]*sizeof(int)); - } - - // Write groups if available - if (hasGroup) { - close(fd); + write(fd, elements, elementSize[i] * 4 * sizeof(int)); + } + } + delete[] elements; + + close(fd); + + // Write vertices + xdmfFile << " " << std::endl + << " " << outfilePrefixShort << "_geometry.bin" << std::endl + << " " << std::endl; + + fd = open((outfilePrefix + "_geometry.bin").c_str(), + O_CREAT | O_WRONLY | O_TRUNC, + S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH); + + float* vertices = new float[maxVertices * 3]; + for (size_t i = 0; i < partitions; i++) { + size_t start[3] = {i, 0, 0}; + size_t size[3] = {1, vertexSize[i], 3}; + nc_get_vara_float(ncFile, ncVarVrtxCoords, start, size, vertices); + + write(fd, vertices, vertexSize[i] * 3 * sizeof(float)); + } + + delete[] vertices; + + close(fd); + + // Write Partition/Boundary type + std::string attributeName; + if (convertBoundaries) + attributeName = "boundary"; + else + attributeName = "partition"; + + xdmfFile << " " << std::endl + << " " << outfilePrefixShort << "_" << attributeName << ".bin" + << std::endl + << " " << std::endl; + + fd = open((outfilePrefix + "_" + attributeName + ".bin").c_str(), + O_CREAT | O_WRONLY | O_TRUNC, + S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH); + + if (convertBoundaries) { + unsigned int* boundaries = new unsigned int[maxElements * 4]; + + int* boundariesOut = new int[nelements]; + unsigned int nextElement = 0; + + for (size_t i = 0; i < partitions; i++) { + size_t start[3] = {i, 0, 0}; + size_t size[3] = {1, elementSize[i], 4}; + nc_get_vara_uint(ncFile, ncVarElemBoundaries, start, size, boundaries); + + for (unsigned int j = 0; j < elementSize[i]; j++) { + for (unsigned int k = 0; k < 4; k++) { + if (boundaries[j * 4 + k] != 0) + boundariesOut[nextElement++] = boundaries[j * 4 + k]; + } + } + } - xdmfFile << " " << std::endl - << " " - << outfilePrefixShort << "_group.bin" << std::endl - << " " << std::endl; + write(fd, boundariesOut, nelements * sizeof(int)); - fd = open((outfilePrefix+"_group.bin").c_str(), O_CREAT | O_WRONLY | O_TRUNC, - S_IRUSR |S_IWUSR |S_IRGRP |S_IWGRP | S_IROTH | S_IWOTH); - - int* groups = new int[maxElements]; - size_t start[3] = {0, 0}; - size_t size[3] = {1, 0}; - for (size_t i = 0; i < partitions; i++) { - start[0] = i; - size[1] = elementSize[i]; - nc_get_vara_int(ncFile, ncVarElemGroup, start, size, groups); - - write(fd, groups, elementSize[i]*sizeof(int)); - } - - delete [] groups; - } - - if (material.size() > 0) { + delete[] boundaries; + delete[] boundariesOut; + } else { + int* partition = new int[maxElements]; + + for (size_t i = 0; i < partitions; i++) { + std::fill_n(partition, elementSize[i], i); + + write(fd, partition, elementSize[i] * sizeof(int)); + } + + // Write groups if available + if (hasGroup) { + close(fd); + + xdmfFile << " " << std::endl + << " " << outfilePrefixShort << "_group.bin" << std::endl + << " " << std::endl; + + fd = open((outfilePrefix + "_group.bin").c_str(), + O_CREAT | O_WRONLY | O_TRUNC, + S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH); + + int* groups = new int[maxElements]; + size_t start[3] = {0, 0}; + size_t size[3] = {1, 0}; + for (size_t i = 0; i < partitions; i++) { + start[0] = i; + size[1] = elementSize[i]; + nc_get_vara_int(ncFile, ncVarElemGroup, start, size, groups); + + write(fd, groups, elementSize[i] * sizeof(int)); + } + + delete[] groups; + } + + if (material.size() > 0) { easi::YAMLParser parser(3); easi::Component* model = parser.parse(material); - - std::vector paramNames{"lambda", "mu", "rho"}; + + std::vector paramNames{"lambda", "mu", "rho"}; int fds[3]; for (unsigned d = 0; d < 3; ++d) { xdmfFile << " " << std::endl - << " " - << outfilePrefixShort << "_" << paramNames[d] << ".bin" << std::endl - << " " << std::endl; - - fds[d] = open((outfilePrefix+"_" + paramNames[d] + ".bin").c_str(), O_CREAT | O_WRONLY | O_TRUNC, - S_IRUSR |S_IWUSR |S_IRGRP |S_IWGRP | S_IROTH | S_IWOTH); + << " " << outfilePrefixShort << "_" << paramNames[d] + << ".bin" << std::endl + << " " << std::endl; + + fds[d] = open((outfilePrefix + "_" + paramNames[d] + ".bin").c_str(), + O_CREAT | O_WRONLY | O_TRUNC, + S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH); } - + easi::ArraysAdapter adapter; double* parameters[3]; for (unsigned d = 0; d < 3; ++d) { parameters[d] = new double[maxElements]; adapter.addBindingPoint(paramNames[d], parameters[d]); } - - int* elements = new int[maxElements*4]; - float* vertices = new float[maxVertices*3]; + + int* elements = new int[maxElements * 4]; + float* vertices = new float[maxVertices * 3]; int* groups; if (hasGroup) { groups = new int[maxElements]; } for (size_t i = 0; i < partitions; ++i) { - size_t startV[3] = {i, 0, 0}; size_t sizeV[3] = {1, elementSize[i], 4}; + size_t startV[3] = {i, 0, 0}; + size_t sizeV[3] = {1, elementSize[i], 4}; nc_get_vara_int(ncFile, ncVarElemVertices, startV, sizeV, elements); - size_t startE[3] = {i, 0, 0}; size_t sizeE[3] = {1, vertexSize[i], 3}; + size_t startE[3] = {i, 0, 0}; + size_t sizeE[3] = {1, vertexSize[i], 3}; nc_get_vara_float(ncFile, ncVarVrtxCoords, startE, sizeE, vertices); if (hasGroup) { - size_t startG[2] = {i, 0}; size_t sizeG[2] = {1, elementSize[i]}; + size_t startG[2] = {i, 0}; + size_t sizeG[2] = {1, elementSize[i]}; nc_get_vara_int(ncFile, ncVarElemGroup, startG, sizeG, groups); } @@ -407,45 +427,43 @@ int main(int argc, char* argv[]) if (hasGroup) { for (unsigned int j = 0; j < elementSize[i]; ++j) { query.group(j) = groups[j]; - } + } } for (unsigned int j = 0; j < elementSize[i]; ++j) { for (unsigned d = 0; d < 3; ++d) { - query.x(j,d) = 0.25f * vertices[ 3*elements[4*j] + d]; + query.x(j, d) = 0.25f * vertices[3 * elements[4 * j] + d]; } for (unsigned v = 1; v < 4; ++v) { for (unsigned d = 0; d < 3; ++d) { - query.x(j,d) += 0.25f * vertices[ 3*elements[4*j+v] + d]; + query.x(j, d) += 0.25f * vertices[3 * elements[4 * j + v] + d]; } } } - + model->evaluate(query, adapter); - + for (unsigned d = 0; d < 3; ++d) { - write(fds[d], parameters[d], elementSize[i]*sizeof(double)); + write(fds[d], parameters[d], elementSize[i] * sizeof(double)); } } - delete [] elements; - delete [] vertices; + delete[] elements; + delete[] vertices; if (hasGroup) { - delete [] groups; + delete[] groups; } for (unsigned d = 0; d < 3; ++d) { delete[] parameters[d]; close(fds[d]); } - + delete model; } - } + } - close(fd); + close(fd); - xdmfFile << " " << std::endl - << " " << std::endl - << "" << std::endl; + xdmfFile << " " << std::endl << " " << std::endl << "" << std::endl; - logInfo() << "Finished successfully"; - return 0; + logInfo() << "Finished successfully"; + return 0; } diff --git a/place_receivers/CMakeLists.txt b/place_receivers/CMakeLists.txt index 9237cbb..082d577 100644 --- a/place_receivers/CMakeLists.txt +++ b/place_receivers/CMakeLists.txt @@ -16,7 +16,7 @@ if(NOT CMAKE_BUILD_TYPE) endif() -add_executable(place_receivers +add_executable(place_receivers ${CMAKE_CURRENT_SOURCE_DIR}/src/Geometry.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/KDTree.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cpp diff --git a/place_receivers/README.md b/place_receivers/README.md index 0e36fc8..595f1b3 100644 --- a/place_receivers/README.md +++ b/place_receivers/README.md @@ -1,15 +1,25 @@ Place receivers =============== -Given receiver locations on a 2D plane, this tool finds the corresponding elevations on the free-surface of the given mesh. + +Given receiver locations on a 2D plane, this tool +finds the corresponding elevations on the free-surface of the given mesh. Method ------ -The program iterates over all faces with free surface boundary condition. For each face, the tool checks if a receiver lies inside the projection of the face on the x-y-plane. Then, the receiver's z coordinate is set such that the receiver lies on the face. Afterwards it is lowered by a configurable amount along the z axis. (The direction is determined by the face normal, i.e. "lowered" means in the opposite direction of the normal.) -Compilation +The program iterates over all faces with free surface boundary +condition. For each face, the tool checks if a receiver lies +inside the projection of the face on the x-y-plane. Then, +the receiver z coordinate is set such that the receiver +lies on the face. Afterwards it is lowered by a configurable +amount along the z axis. (The direction is determined by +the face normal, i.e. "lowered" means in the opposite +direction of the normal.) + +Compilation ----- -``` +```bash mkdir build && cd build cmake .. make @@ -18,29 +28,45 @@ make Usage ----- +```bash +./place_receivers -d DEPTH -r RECEIVERS -m MESH -o OUTPUT ``` -place_receivers -d DEPTH -r RECEIVERS -m MESH -o OUTPUT -``` -- DEPTH should be a positive number that specifies the amount of meters that the receiver shall be placed below the free surface. -- RECEIVERS is the file name of a double column, space separated file where the first column denotes the x coordinate and the second column the y coordinate. For example: +`DEPTH` should be a positive number that specifies +the amount of meters that the receiver shall be +placed below the free surface. + +`RECEIVERS` is the file name of a double column, +space separated file where the first column denotes +the x coordinate and the second column the y coordinate. +For example: -
+```bash
 8.251411e+05 5.900044e+05
 7.654161e+05 3.274531e+05
 7.794138e+05 5.097186e+05
-
+``` -One may specify the depth per receiver with an optional third column, where the depth is the maximum of DEPTH and the third column. For example with DEPTH=0.1: +One may specify the depth per receiver with an +optional third column, where the depth is the +maximum of DEPTH and the third column. For +example with DEPTH=0.1: -
+```bash
 8.251411e+05 5.900044e+05 0.0   # -> depth 0.1
 7.654161e+05 3.274531e+05 100.0 # -> depth 100
-
+``` -- MESH is the file name of the mesh. -- OUTPUT is the file name of the output file. +`MESH` is the file name of the mesh. + +`OUTPUT` is the file name of the output file. Limitations ----------- -The tools does only work if faces with free surface boundary condition are present and if their normals look in +z or -z direction. (This is usually the case if free surfaces due to topography and e.g. enu, ned, seu axis orientations are used.) That is, the tool fails if all respective faces have normals with z=0. + +The tools does only work if faces with free surface +boundary condition are present and if their normals +look in +z or -z direction. (This is usually the case +if free surfaces due to topography and e.g. enu, ned, +seu axis orientations are used.) That is, the tool +fails if all respective faces have normals with z=0. diff --git a/place_receivers/src/Geometry.cpp b/place_receivers/src/Geometry.cpp index ece173e..c5f20de 100644 --- a/place_receivers/src/Geometry.cpp +++ b/place_receivers/src/Geometry.cpp @@ -2,23 +2,24 @@ * @file * This file is part of SeisSol. * - * @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) - * @author Thomas Ulrich + * @author Carsten Uphoff (c.uphoff AT tum.de, + * http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) + * @author Thomas Ulrich * * @section LICENSE * Copyright (c) 2016, SeisSol Group * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, * this list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from this * software without specific prior written permission. @@ -40,13 +41,13 @@ #include "Geometry.h" -#include #include +#include #ifdef USE_NETCDF -int const FACE2NODES[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}}; +const int FACE2NODES[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}}; #else -int const FACE2NODES[4][3] = {{0, 2, 1}, {0, 1, 3}, {1, 2, 3}, {0, 3, 2}}; +const int FACE2NODES[4][3] = {{0, 2, 1}, {0, 1, 3}, {1, 2, 3}, {0, 3, 2}}; #endif struct Support { @@ -60,10 +61,8 @@ struct Support { limits[1][0] = std::numeric_limits::max(); limits[1][1] = std::numeric_limits::min(); } - - double operator()(int splitdim, int side) const { - return limits[splitdim][side]; - } + + double operator()(int splitdim, int side) const { return limits[splitdim][side]; } }; struct Action { @@ -73,41 +72,44 @@ struct Action { double faceNormal[3]; double faceDist; double depth; - + void determineNormals() { for (unsigned side = 0; side < 3; ++side) { double* v0 = vertices[side].coords; - double* v1 = vertices[(side+1)%3].coords; - double* vtest = vertices[(side+2)%3].coords; + double* v1 = vertices[(side + 1) % 3].coords; + double* vtest = vertices[(side + 2) % 3].coords; // normal = (y,-x) - normals[side][0] = v1[1]-v0[1]; - normals[side][1] = -v1[0]+v0[0]; + normals[side][0] = v1[1] - v0[1]; + normals[side][1] = -v1[0] + v0[0]; dist[side] = normals[side][0] * v0[0] + normals[side][1] * v0[1]; - + if (normals[side][0] * vtest[0] + normals[side][1] * vtest[1] > dist[side]) { normals[side][0] *= -1.; normals[side][1] *= -1.; dist[side] *= -1; } } - + double a[3]; double b[3]; for (unsigned d = 0; d < 3; ++d) { - a[d] = vertices[1].coords[d]-vertices[0].coords[d]; - b[d] = vertices[2].coords[d]-vertices[0].coords[d]; + a[d] = vertices[1].coords[d] - vertices[0].coords[d]; + b[d] = vertices[2].coords[d] - vertices[0].coords[d]; } - faceNormal[0] = a[1]*b[2] - a[2]*b[1]; - faceNormal[1] = a[2]*b[0] - a[0]*b[2]; - faceNormal[2] = a[0]*b[1] - a[1]*b[0]; - faceDist = faceNormal[0] * vertices[0].coords[0] + faceNormal[1] * vertices[0].coords[1] + faceNormal[2] * vertices[0].coords[2]; + faceNormal[0] = a[1] * b[2] - a[2] * b[1]; + faceNormal[1] = a[2] * b[0] - a[0] * b[2]; + faceNormal[2] = a[0] * b[1] - a[1] * b[0]; + faceDist = faceNormal[0] * vertices[0].coords[0] + faceNormal[1] * vertices[0].coords[1] + + faceNormal[2] * vertices[0].coords[2]; } - + void operator()(Receiver& receiver) { bool inside = true; for (unsigned side = 0; side < 3; ++side) { - inside = inside && (normals[side][0] * receiver.point.x + normals[side][1] * receiver.point.y <= dist[side]); + inside = + inside && + (normals[side][0] * receiver.point.x + normals[side][1] * receiver.point.y <= dist[side]); } if (inside) { receiver.found = true; @@ -115,7 +117,9 @@ struct Action { if (!std::isnan(receiver.point.z)) { local_depth = std::max(depth, std::fabs(receiver.point.z)); } - receiver.point.z = (faceDist - faceNormal[0] * receiver.point.x - faceNormal[1] * receiver.point.y) / faceNormal[2]; + receiver.point.z = + (faceDist - faceNormal[0] * receiver.point.x - faceNormal[1] * receiver.point.y) / + faceNormal[2]; if (faceNormal[2] >= 0) { receiver.point.z -= local_depth; } else { @@ -125,16 +129,18 @@ struct Action { } }; -void setElevation(int partition, double depth, Mesh const& mesh, KDTree& tree) { +void setElevation(int partition, double depth, const Mesh& mesh, KDTree& tree) { for (unsigned element = 0; element < mesh.elementSize[partition]; ++element) { for (unsigned face = 0; face < 4; ++face) { // Consider only free surface boundaries - if (mesh.elementBoundaries[4*element + face] == 1) { + if (mesh.elementBoundaries[4 * element + face] == 1) { Action act; act.depth = depth; Support sup; for (unsigned node = 0; node < 3; ++node) { - double* vbegin = &mesh.vertexCoordinates[3*mesh.elementVertices[ 4*element + FACE2NODES[face][node] ]]; + double* vbegin = + &mesh.vertexCoordinates[3 * + mesh.elementVertices[4 * element + FACE2NODES[face][node]]]; act.vertices[node].coords[0] = *(vbegin); act.vertices[node].coords[1] = *(vbegin + 1); act.vertices[node].coords[2] = *(vbegin + 2); diff --git a/place_receivers/src/Geometry.h b/place_receivers/src/Geometry.h index a83adaa..e3688a1 100644 --- a/place_receivers/src/Geometry.h +++ b/place_receivers/src/Geometry.h @@ -2,22 +2,23 @@ * @file * This file is part of SeisSol. * - * @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) + * @author Carsten Uphoff (c.uphoff AT tum.de, + * http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) * * @section LICENSE * Copyright (c) 2016, SeisSol Group * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, * this list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from this * software without specific prior written permission. @@ -40,10 +41,9 @@ #ifndef GEOMETRY_H_ #define GEOMETRY_H_ -#include "Reader.h" #include "KDTree.h" +#include "Reader.h" -void setElevation(int partition, double depth, Mesh const& mesh, KDTree& tree); - +void setElevation(int partition, double depth, const Mesh& mesh, KDTree& tree); #endif diff --git a/place_receivers/src/KDTree.cpp b/place_receivers/src/KDTree.cpp index 13f9e91..7beb762 100644 --- a/place_receivers/src/KDTree.cpp +++ b/place_receivers/src/KDTree.cpp @@ -1,96 +1,91 @@ #include "KDTree.h" #include -#include #include +#include + +KDTree::KDTree(const std::vector& points, int maxLeafSize) : p(2), maxLeafN(maxLeafSize) { + const auto n = points.size(); -KDTree::KDTree(std::vector const& points, int maxLeafSize) - : p(2), maxLeafN(maxLeafSize) -{ - const auto n = points.size(); - - // Copy data and change row-major to column-major storage. - data.resize(n); - for (std::size_t i = 0; i < n; ++i) { - data[i].point.x = points[i].x; - data[i].point.y = points[i].y; - data[i].point.z = points[i].z; - data[i].found = false; - } - - idx.resize(n); - for (int i = 0; i < n; ++i) { - idx[i] = i; - } + // Copy data and change row-major to column-major storage. + data.resize(n); + for (std::size_t i = 0; i < n; ++i) { + data[i].point.x = points[i].x; + data[i].point.y = points[i].y; + data[i].point.z = points[i].z; + data[i].found = false; + } - int maxHeight = 1 + ceil(log2(n / static_cast(maxLeafN))); - int maxNodes = (1 << maxHeight) - 1; - nodes.resize(maxNodes); - nodes[0].start = 0; - nodes[0].n = n; - - buildTree(0, 0); + idx.resize(n); + for (int i = 0; i < n; ++i) { + idx[i] = i; + } + + int maxHeight = 1 + ceil(log2(n / static_cast(maxLeafN))); + int maxNodes = (1 << maxHeight) - 1; + nodes.resize(maxNodes); + nodes[0].start = 0; + nodes[0].n = n; + + buildTree(0, 0); } -void KDTree::swap(int i, int j) -{ - if (i != j) { - std::swap(idx[i], idx[j]); +void KDTree::swap(int i, int j) { + if (i != j) { + std::swap(idx[i], idx[j]); std::swap(data[i], data[j]); - } + } } -int KDTree::partition(int left, int right, int pivotIdx, int splitdim) -{ - double pivot = data[pivotIdx].point.coords[splitdim]; - int st = left; - for (int i = left; i < right; ++i) { - if (data[i].point.coords[splitdim] < pivot) { - swap(st, i); - ++st; - } - } - swap(st, right); - return st; +int KDTree::partition(int left, int right, int pivotIdx, int splitdim) { + double pivot = data[pivotIdx].point.coords[splitdim]; + int st = left; + for (int i = left; i < right; ++i) { + if (data[i].point.coords[splitdim] < pivot) { + swap(st, i); + ++st; + } + } + swap(st, right); + return st; } -void KDTree::buildTree(int k, int splitdim) -{ - Node& node = nodes[k]; - node.splitdim = splitdim; - if (node.n > maxLeafN) { - int half = (node.n % 2 != 0) ? (node.n + 1)/2 : node.n/2; - int l = node.start; - int r = l + node.n - 1; - int median_idx = l + half; - if (l != r) { - int pivotIdx; - while (true) { - pivotIdx = r; - pivotIdx = partition(l, r, pivotIdx, splitdim); - if (median_idx == pivotIdx) { - break; - } else if (median_idx < pivotIdx) { - r = pivotIdx - 1; - } else { - l = pivotIdx + 1; - } - } - } - node.pivot = data[median_idx].point.coords[splitdim]; +void KDTree::buildTree(int k, int splitdim) { + Node& node = nodes[k]; + node.splitdim = splitdim; + if (node.n > maxLeafN) { + int half = (node.n % 2 != 0) ? (node.n + 1) / 2 : node.n / 2; + int l = node.start; + int r = l + node.n - 1; + int median_idx = l + half; + if (l != r) { + int pivotIdx; + while (true) { + pivotIdx = r; + pivotIdx = partition(l, r, pivotIdx, splitdim); + if (median_idx == pivotIdx) { + break; + } else if (median_idx < pivotIdx) { + r = pivotIdx - 1; + } else { + l = pivotIdx + 1; + } + } + } + node.pivot = data[median_idx].point.coords[splitdim]; + + Node& left = nodes[leftChild(k)]; + left.start = node.start; + left.n = half; - Node& left = nodes[leftChild(k)]; - left.start = node.start; - left.n = half; + Node& right = nodes[rightChild(k)]; + right.start = median_idx; + right.n = node.n - half; - Node& right = nodes[rightChild(k)]; - right.start = median_idx; - right.n = node.n - half; - - int nextSplitdim = (splitdim + 1) % p; - buildTree(leftChild(k), nextSplitdim); - buildTree(rightChild(k), nextSplitdim); - } else { - node.isLeaf = true; - } + int nextSplitdim = (splitdim + 1) % p; + buildTree(leftChild(k), nextSplitdim); + buildTree(rightChild(k), nextSplitdim); + } else { + node.isLeaf = true; + } } diff --git a/place_receivers/src/KDTree.h b/place_receivers/src/KDTree.h index 180aec6..b549ec8 100644 --- a/place_receivers/src/KDTree.h +++ b/place_receivers/src/KDTree.h @@ -2,22 +2,23 @@ * @file * This file is part of SeisSol. * - * @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) + * @author Carsten Uphoff (c.uphoff AT tum.de, + * http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) * * @section LICENSE * Copyright (c) 2016, SeisSol Group * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, * this list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from this * software without specific prior written permission. @@ -52,67 +53,65 @@ union Point { }; struct Receiver { - Point point; - bool found{false}; + Point point; + bool found{false}; }; class KDTree { -public: - KDTree(std::vector const& points, int maxLeafSize); - - template - void search(Support const& support, Action& action) - { - searchTree(0, support, action); - } - - inline std::vector points() const { return data; } - inline int index(int r) const { return idx[r]; } - inline int numPoints() const { return nodes[0].n;} + public: + KDTree(const std::vector& points, int maxLeafSize); + + template + void search(const Support& support, Action& action) { + searchTree(0, support, action); + } -private: - int leftChild(int k) const { return 2*k + 1; } - int rightChild(int k) const { return 2*k + 2; } + inline std::vector points() const { return data; } + inline int index(int r) const { return idx[r]; } + inline int numPoints() const { return nodes[0].n; } - struct Node { + private: + int leftChild(int k) const { return 2 * k + 1; } + int rightChild(int k) const { return 2 * k + 2; } + + struct Node { Node() : isLeaf(false) {} - double pivot; - int start; - int n; - int splitdim; - bool isLeaf; - }; - std::vector nodes; - - void swap(int i, int j); - int partition(int left, int right, int pivotIdx, int splitdim); - void buildTree(int k, int splitdim); - - template - void searchTree(int k, Support const& support, Action& action); + double pivot; + int start; + int n; + int splitdim; + bool isLeaf; + }; + std::vector nodes; + + void swap(int i, int j); + int partition(int left, int right, int pivotIdx, int splitdim); + void buildTree(int k, int splitdim); + + template + void searchTree(int k, const Support& support, Action& action); - std::vector data; - std::vector idx; - int p; - int maxLeafN; + std::vector data; + std::vector idx; + int p; + int maxLeafN; }; -template -void KDTree::searchTree(int k, Support const& support, Action& action) -{ - const auto& node = nodes[k]; - if (node.isLeaf) { - for (int i = node.start; i < node.start + node.n; ++i) { - action(data[i]); - } - } else { - if (support(node.splitdim, 0) <= node.pivot) { - searchTree(leftChild(k), support, action); - } - if (support(node.splitdim, 1) >= node.pivot) { - searchTree(rightChild(k), support, action); - } - } +template +void KDTree::searchTree(int k, const Support& support, Action& action) { + const auto& node = nodes[k]; + if (node.isLeaf) { + for (int i = node.start; i < node.start + node.n; ++i) { + action(data[i]); + } + } else { + if (support(node.splitdim, 0) <= node.pivot) { + searchTree(leftChild(k), support, action); + } + if (support(node.splitdim, 1) >= node.pivot) { + searchTree(rightChild(k), support, action); + } + } } #endif diff --git a/place_receivers/src/Reader.cpp b/place_receivers/src/Reader.cpp index 42a997b..4fa9d2a 100644 --- a/place_receivers/src/Reader.cpp +++ b/place_receivers/src/Reader.cpp @@ -2,23 +2,24 @@ * @file * This file is part of SeisSol. * - * @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) - * @author Thomas Ulrich + * @author Carsten Uphoff (c.uphoff AT tum.de, + * http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) + * @author Thomas Ulrich * * @section LICENSE * Copyright (c) 2016, SeisSol Group * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, * this list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from this * software without specific prior written permission. @@ -39,19 +40,21 @@ **/ #include "Reader.h" + #include "PUML/PUML.h" -#include + #include -#include -#include +#include #include +#include +#include #include #ifdef USE_NETCDF #include #endif -std::vector readReceiverFile(std::string const& fileName) { +std::vector readReceiverFile(const std::string& fileName) { std::vector locations; std::ifstream in(fileName.c_str()); std::string line; @@ -69,156 +72,154 @@ std::vector readReceiverFile(std::string const& fileName) { exit(-1); } } - + return locations; } -void writeReceiverFile(KDTree const& tree, std::string const& fileName) { +void writeReceiverFile(const KDTree& tree, const std::string& fileName) { std::ofstream out(fileName.c_str()); out << std::scientific << std::setprecision(16); - + int failureCounter = 0; const auto& receivers = tree.points(); std::vector sortedPoints(receivers.size()); for (unsigned p = 0; p < tree.numPoints(); ++p) { sortedPoints[tree.index(p)] = receivers[p]; } - for (auto const& receiver : sortedPoints) { + for (const auto& receiver : sortedPoints) { if (!std::isnan(receiver.point.z) && receiver.found) { out << receiver.point.x << " " << receiver.point.y << " " << receiver.point.z << std::endl; } else { - std::cerr << "Warning: Did not find elevation for receiver at (" << receiver.point.x << ", " << receiver.point.y << ")." << std::endl; + std::cerr << "Warning: Did not find elevation for receiver at (" << receiver.point.x << ", " + << receiver.point.y << ")." << std::endl; ++failureCounter; } } if (failureCounter > 0) { - std::cerr << failureCounter << " points were not written due to missing elevation." << std::endl; + std::cerr << failureCounter << " points were not written due to missing elevation." + << std::endl; } } #ifdef USE_NETCDF -void check_err(const int stat, const int line, const char *file) { +void check_err(const int stat, const int line, const char* file) { if (stat != NC_NOERR) { - fprintf(stderr,"line %d of %s: %s\n", line, file, nc_strerror(stat)); + fprintf(stderr, "line %d of %s: %s\n", line, file, nc_strerror(stat)); fflush(stderr); exit(-1); } } #endif -Mesh::Mesh(std::string const& fileName){ +Mesh::Mesh(const std::string& fileName) { #ifndef USE_NETCDF PUML::TETPUML puml; puml.open((fileName + ":/connect").c_str(), (fileName + ":/geometry").c_str()); puml.addData((fileName + ":/boundary").c_str(), PUML::CELL); - size_t nElements = puml.numOriginalCells(); size_t nVertex = puml.numOriginalVertices(); - elementVertices = new int[nElements*4]; - vertexCoordinates = new double[nVertex*3]; - elementBoundaries = new int[nElements*4]; + elementVertices = new int[nElements * 4]; + vertexCoordinates = new double[nVertex * 3]; + elementBoundaries = new int[nElements * 4]; partitions = 1; elementSize = new int[partitions]; vertexSize = new int[partitions]; elementSize[0] = nElements; - vertexSize[0]= nVertex; - + vertexSize[0] = nVertex; - double* vertexCoordinate = new double [3]; + double* vertexCoordinate = new double[3]; for (unsigned int i = 0; i < nVertex; i++) { - vertexCoordinate = (double*) puml.originalVertices()[i]; + vertexCoordinate = (double*)puml.originalVertices()[i]; for (unsigned int j = 0; j < 3; j++) { - vertexCoordinates[3*i+j] = vertexCoordinate[j]; + vertexCoordinates[3 * i + j] = vertexCoordinate[j]; } } - unsigned long * elementVertice = new unsigned long [4]; + unsigned long* elementVertice = new unsigned long[4]; for (unsigned int i = 0; i < nElements; i++) { - elementVertice = (unsigned long *) puml.originalCells()[i]; + elementVertice = (unsigned long*)puml.originalCells()[i]; for (unsigned int j = 0; j < 4; j++) { - elementVertices[4*i+j] = (int) elementVertice[j]; + elementVertices[4 * i + j] = (int)elementVertice[j]; } } const int* boundaryCond = puml.cellData(0); for (unsigned int i = 0; i < nElements; i++) { for (unsigned int face = 0; face < 4; face++) { - elementBoundaries[4*i + face] = (boundaryCond[i] >> (face*8)) & 0xFF; + elementBoundaries[4 * i + face] = (boundaryCond[i] >> (face * 8)) & 0xFF; } } - + #else int stat; - size_t maxVertices; - size_t maxElements; + size_t maxVertices; + size_t maxElements; stat = nc_open(fileName.c_str(), 0, &ncid); - check_err(stat,__LINE__,__FILE__); - - int ncDimPart; - stat = nc_inq_dimid(ncid, "partitions", &ncDimPart); - check_err(stat,__LINE__,__FILE__); - stat = nc_inq_dimlen(ncid, ncDimPart, &partitions); - check_err(stat,__LINE__,__FILE__); - - int ncDimVrtx; - stat = nc_inq_dimid(ncid, "vertices", &ncDimVrtx); - check_err(stat,__LINE__,__FILE__); - stat = nc_inq_dimlen(ncid, ncDimVrtx, &maxVertices); - check_err(stat,__LINE__,__FILE__); - - int ncDimElem; - stat = nc_inq_dimid(ncid, "elements", &ncDimElem); - check_err(stat,__LINE__,__FILE__); - stat = nc_inq_dimlen(ncid, ncDimElem, &maxElements); - check_err(stat,__LINE__,__FILE__); + check_err(stat, __LINE__, __FILE__); + + int ncDimPart; + stat = nc_inq_dimid(ncid, "partitions", &ncDimPart); + check_err(stat, __LINE__, __FILE__); + stat = nc_inq_dimlen(ncid, ncDimPart, &partitions); + check_err(stat, __LINE__, __FILE__); + + int ncDimVrtx; + stat = nc_inq_dimid(ncid, "vertices", &ncDimVrtx); + check_err(stat, __LINE__, __FILE__); + stat = nc_inq_dimlen(ncid, ncDimVrtx, &maxVertices); + check_err(stat, __LINE__, __FILE__); + + int ncDimElem; + stat = nc_inq_dimid(ncid, "elements", &ncDimElem); + check_err(stat, __LINE__, __FILE__); + stat = nc_inq_dimlen(ncid, ncDimElem, &maxElements); + check_err(stat, __LINE__, __FILE__); int ncVarElemSize; - stat = nc_inq_varid(ncid, "element_size", &ncVarElemSize); - check_err(stat,__LINE__,__FILE__); + stat = nc_inq_varid(ncid, "element_size", &ncVarElemSize); + check_err(stat, __LINE__, __FILE__); - stat = nc_inq_varid(ncid, "element_vertices", &ncVarElemVertices); - check_err(stat,__LINE__,__FILE__); + stat = nc_inq_varid(ncid, "element_vertices", &ncVarElemVertices); + check_err(stat, __LINE__, __FILE__); - stat = nc_inq_varid(ncid, "element_boundaries", &ncVarElemBoundaries); - check_err(stat,__LINE__,__FILE__); + stat = nc_inq_varid(ncid, "element_boundaries", &ncVarElemBoundaries); + check_err(stat, __LINE__, __FILE__); int ncVarVrtxSize; - stat = nc_inq_varid(ncid, "vertex_size", &ncVarVrtxSize); - check_err(stat,__LINE__,__FILE__); + stat = nc_inq_varid(ncid, "vertex_size", &ncVarVrtxSize); + check_err(stat, __LINE__, __FILE__); - stat = nc_inq_varid(ncid, "vertex_coordinates", &ncVarVrtxCoords); - check_err(stat,__LINE__,__FILE__); + stat = nc_inq_varid(ncid, "vertex_coordinates", &ncVarVrtxCoords); + check_err(stat, __LINE__, __FILE__); - elementSize = new int[partitions]; - stat = nc_get_var_int(ncid, ncVarElemSize, elementSize); - check_err(stat,__LINE__,__FILE__); + elementSize = new int[partitions]; + stat = nc_get_var_int(ncid, ncVarElemSize, elementSize); + check_err(stat, __LINE__, __FILE__); - vertexSize = new int[partitions]; - stat = nc_get_var_int(ncid, ncVarVrtxSize, vertexSize); - check_err(stat,__LINE__,__FILE__); - - vertexCoordinates = new double[3*maxVertices]; - elementBoundaries = new int[4*maxElements]; - elementVertices = new int[4*maxElements]; + vertexSize = new int[partitions]; + stat = nc_get_var_int(ncid, ncVarVrtxSize, vertexSize); + check_err(stat, __LINE__, __FILE__); + + vertexCoordinates = new double[3 * maxVertices]; + elementBoundaries = new int[4 * maxElements]; + elementVertices = new int[4 * maxElements]; #endif } - - Mesh::~Mesh() { #ifdef USE_NETCDF int stat = nc_close(ncid); - check_err(stat,__LINE__,__FILE__); + check_err(stat, __LINE__, __FILE__); #endif - + delete[] elementSize; delete[] vertexSize; delete[] elementBoundaries; @@ -231,18 +232,17 @@ void Mesh::readPartition(int partition) { int stat; size_t elementsStart[3] = {partition, 0, 0}; size_t elementsSize[3] = {1, elementSize[partition], 4}; - + stat = nc_get_vara_int(ncid, ncVarElemVertices, elementsStart, elementsSize, elementVertices); - check_err(stat,__LINE__,__FILE__); - + check_err(stat, __LINE__, __FILE__); + stat = nc_get_vara_int(ncid, ncVarElemBoundaries, elementsStart, elementsSize, elementBoundaries); - check_err(stat,__LINE__,__FILE__); - - + check_err(stat, __LINE__, __FILE__); + size_t verticesStart[3] = {partition, 0, 0}; size_t verticesSize[3] = {1, vertexSize[partition], 3}; - + stat = nc_get_vara_double(ncid, ncVarVrtxCoords, verticesStart, verticesSize, vertexCoordinates); - check_err(stat,__LINE__,__FILE__); + check_err(stat, __LINE__, __FILE__); #endif } diff --git a/place_receivers/src/Reader.h b/place_receivers/src/Reader.h index 41fc7a6..17e7b3d 100644 --- a/place_receivers/src/Reader.h +++ b/place_receivers/src/Reader.h @@ -2,22 +2,23 @@ * @file * This file is part of SeisSol. * - * @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) + * @author Carsten Uphoff (c.uphoff AT tum.de, + * http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) * * @section LICENSE * Copyright (c) 2016, SeisSol Group * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, * this list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from this * software without specific prior written permission. @@ -40,33 +41,33 @@ #ifndef READER_H_ #define READER_H_ +#include "KDTree.h" + #include #include -#include "KDTree.h" - class Mesh { -public: - explicit Mesh(std::string const& fileName); + public: + explicit Mesh(const std::string& fileName); ~Mesh(); - + void readPartition(int partition); - size_t partitions; + size_t partitions; int* elementSize; int* vertexSize; int* elementBoundaries; int* elementVertices; double* vertexCoordinates; - -private: + + private: int ncid; - int ncVarElemVertices; - int ncVarElemBoundaries; - int ncVarVrtxCoords; + int ncVarElemVertices; + int ncVarElemBoundaries; + int ncVarVrtxCoords; }; -std::vector readReceiverFile(std::string const& fileName); -void writeReceiverFile(KDTree const& tree, std::string const& fileName); +std::vector readReceiverFile(const std::string& fileName); +void writeReceiverFile(const KDTree& tree, const std::string& fileName); #endif diff --git a/place_receivers/src/main.cpp b/place_receivers/src/main.cpp index b5af6a4..f30a67a 100644 --- a/place_receivers/src/main.cpp +++ b/place_receivers/src/main.cpp @@ -2,23 +2,24 @@ * @file * This file is part of SeisSol. * - * @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) + * @author Carsten Uphoff (c.uphoff AT tum.de, + * http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) * @author Thomas Ulrich * * @section LICENSE * Copyright (c) 2016, SeisSol Group * All rights reserved. - * + * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, * this list of conditions and the following disclaimer. - * + * * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * 3. Neither the name of the copyright holder nor the names of its * contributors may be used to endorse or promote products derived from this * software without specific prior written permission. @@ -38,25 +39,24 @@ * @section DESCRIPTION **/ +#include "Geometry.h" +#include "Reader.h" + #include #include #include -#include "Reader.h" -#include "Geometry.h" - -int main(int argc, char** argv) -{ +int main(int argc, char** argv) { utils::Args args; args.addOption("depth", 'd', "Distance of the receivers to the free surface"); args.addOption("receivers", 'r', "Receiver locations"); args.addOption("mesh", 'm', "mesh file"); args.addOption("output", 'o', "Receiver output file"); - if (args.parse(argc, argv) != utils::Args::Success) { + if (args.parse(argc, argv) != utils::Args::Success) { return -1; } - + double depth = args.getArgument("depth"); std::string receiverFile = args.getArgument("receivers"); std::string meshFile = args.getArgument("mesh"); @@ -70,19 +70,20 @@ int main(int argc, char** argv) } KDTree tree(receivers, 1); - + Mesh mesh(meshFile); - + for (unsigned p = 0; p < mesh.partitions; ++p) { if (p % 100 == 0) { - std::cout << "Processing partition " << p << " to " << std::min(mesh.partitions, static_cast(p+99)) << std::endl; + std::cout << "Processing partition " << p << " to " + << std::min(mesh.partitions, static_cast(p + 99)) << std::endl; } mesh.readPartition(p); setElevation(p, depth, mesh, tree); } - writeReceiverFile(tree, receiverOutputFile); - + writeReceiverFile(tree, receiverOutputFile); + return 0; } diff --git a/remove-groups/README.md b/remove-groups/README.md index 26bc789..307a713 100644 --- a/remove-groups/README.md +++ b/remove-groups/README.md @@ -2,6 +2,8 @@ Remove groups =============== A script that removes certain groups of the mesh. -At the interface between removed cells and other cells, a free surface boundary condition is inserted. +At the interface between removed cells and other +cells, a free surface boundary condition is inserted. -The main use-case is to remove a water layer from a fully-coupled mesh without the need for remeshing. \ No newline at end of file +The main use-case is to remove a water layer from +a fully-coupled mesh without the need for remeshing. diff --git a/remove-groups/src/Filter.h b/remove-groups/src/Filter.h index 7fdd6de..60b4426 100644 --- a/remove-groups/src/Filter.h +++ b/remove-groups/src/Filter.h @@ -1,11 +1,11 @@ #ifndef REMOVE_WATERLAYER_FILTER_H #define REMOVE_WATERLAYER_FILTER_H -#include #include "Reader.h" -void removeGroups(Mesh &mesh, - const std::set &groupsToRemove) { +#include + +void removeGroups(Mesh& mesh, const std::set& groupsToRemove) { int removedCells = 0; // Map old element/vertex id -> new element/vertex id std::unordered_map elementMap; @@ -18,23 +18,23 @@ void removeGroups(Mesh &mesh, if (groupsToRemove.count(mesh.elementGroups[cell]) > 0) { // Remove cell ++removedCells; - const auto &neighbors = mesh.neighbors[cell]; - // For all neighbors, set face type to free surface in slot of current cell + const auto& neighbors = mesh.neighbors[cell]; + // For all neighbors, set face type to free surface in slot of current + // cell for (auto neighborId : neighbors) { if (neighborId) { - auto &neighborBnds = mesh.elementBoundaries[*neighborId]; - auto &neighborNeighbors = mesh.neighbors[*neighborId]; + auto& neighborBnds = mesh.elementBoundaries[*neighborId]; + auto& neighborNeighbors = mesh.neighbors[*neighborId]; // Find which nth neighbor we are of our neighbor auto faceId = -1; for (auto j = 0; j < 4; ++j) { - auto &nn = neighborNeighbors[j]; + auto& nn = neighborNeighbors[j]; if (nn && *nn == cell) { faceId = j; } } neighborBnds[faceId] = 1; // Change to free surface bc - } } } else { @@ -48,7 +48,6 @@ void removeGroups(Mesh &mesh, vertexMap[oldVertexId] = vertexCounter++; } } - } if (cell % (mesh.elementSize / 10) == 0) { std::cout << "Processed " << cell << "\t cells out of " << mesh.elementSize << std::endl; @@ -64,7 +63,8 @@ void removeGroups(Mesh &mesh, auto newElementNeighbors = std::vector(elementCounter); for (auto cell = 0u; cell < mesh.elementSize; ++cell) { - if (elementMap.count(cell) == 0) continue; + if (elementMap.count(cell) == 0) + continue; auto newCellId = elementMap[cell]; newElementBoundaries[newCellId] = mesh.elementBoundaries[cell]; newElementGroups[newCellId] = mesh.elementGroups[cell]; @@ -96,7 +96,8 @@ void removeGroups(Mesh &mesh, auto newVertices = std::vector(vertexCounter * 3); for (unsigned vertex = 0; vertex < mesh.vertexSize; ++vertex) { - if (vertexMap.find(vertex) == vertexMap.end()) continue; + if (vertexMap.find(vertex) == vertexMap.end()) + continue; auto newVertexId = vertexMap[vertex]; for (int i = 0; i < 3; ++i) { newVertices[newVertexId * 3 + i] = mesh.vertices[vertex * 3 + i]; @@ -104,8 +105,6 @@ void removeGroups(Mesh &mesh, } mesh.vertexSize = vertexCounter; mesh.vertices = std::move(newVertices); - - } -#endif //REMOVE_WATERLAYER_FILTER_H +#endif // REMOVE_WATERLAYER_FILTER_H diff --git a/remove-groups/src/Reader.cpp b/remove-groups/src/Reader.cpp index e4e44c7..a340ee4 100644 --- a/remove-groups/src/Reader.cpp +++ b/remove-groups/src/Reader.cpp @@ -1,13 +1,13 @@ #include "Reader.h" -#include -#include - -#include "PUML/PUML.h" #include "PUML/Neighbor.h" +#include "PUML/PUML.h" #include "mpi.h" -Mesh::Mesh(const std::string &fileName) { +#include +#include + +Mesh::Mesh(const std::string& fileName) { PUML::TETPUML puml; puml.setComm(MPI_COMM_WORLD); puml.open((fileName + ":/connect").c_str(), (fileName + ":/geometry").c_str()); @@ -30,8 +30,8 @@ Mesh::Mesh(const std::string &fileName) { neighbors = std::vector(elementSize); std::array neighborsRaw{}; - for (auto &vertex : verticesPuml) { - const auto *coordinate = vertex.coordinate(); + for (auto& vertex : verticesPuml) { + const auto* coordinate = vertex.coordinate(); for (int i = 0; i < 3; ++i) { vertices[verticesOffset++] = coordinate[i]; } @@ -51,15 +51,12 @@ Mesh::Mesh(const std::string &fileName) { PUML::Neighbor::face(puml, cell, neighborsRaw.data()); neighbors_t curNeighbors{}; for (int face = 0; face < 4; ++face) { - curNeighbors[face] = (neighborsRaw[face] >= 0) ? std::optional(neighborsRaw[face]) : std::nullopt; + curNeighbors[face] = + (neighborsRaw[face] >= 0) ? std::optional(neighborsRaw[face]) : std::nullopt; elementBoundaries[cell][face] = decodeBoundaryCondition(puml.cellData(1)[cell], face); } neighbors[cell] = curNeighbors; } - - } -int decodeBoundaryCondition(int encoded, int faceId) { - return (encoded >> (faceId * 8)) & 0xFF; -} \ No newline at end of file +int decodeBoundaryCondition(int encoded, int faceId) { return (encoded >> (faceId * 8)) & 0xFF; } diff --git a/remove-groups/src/Reader.h b/remove-groups/src/Reader.h index e8fe557..be7ae02 100644 --- a/remove-groups/src/Reader.h +++ b/remove-groups/src/Reader.h @@ -1,27 +1,26 @@ #ifndef READER_H_ #define READER_H_ +#include +#include #include #include -#include -#include class Mesh { -public: - using neighbor_t = std::optional; - using neighbors_t = std::array; - using boundaries_t = std::array; - - explicit Mesh(const std::string &fileName); + public: + using neighbor_t = std::optional; + using neighbors_t = std::array; + using boundaries_t = std::array; - unsigned int elementSize; - unsigned int vertexSize; - std::vector elementBoundaries; - std::vector elementGroups; - std::vector connect; - std::vector vertices; - std::vector neighbors; + explicit Mesh(const std::string& fileName); + unsigned int elementSize; + unsigned int vertexSize; + std::vector elementBoundaries; + std::vector elementGroups; + std::vector connect; + std::vector vertices; + std::vector neighbors; }; int decodeBoundaryCondition(int encoded, int faceId); diff --git a/remove-groups/src/Writer.h b/remove-groups/src/Writer.h index d1aef4f..cb42bd5 100644 --- a/remove-groups/src/Writer.h +++ b/remove-groups/src/Writer.h @@ -1,19 +1,16 @@ #ifndef REMOVE_WATERLAYER_WRITER_H #define REMOVE_WATERLAYER_WRITER_H -#include -#include #include "Reader.h" #include "hdf5.h" -enum class Hdf5DataType { - hdf5UInt64, - hdf5Int32, - hdf5Double -}; +#include +#include + +enum class Hdf5DataType { hdf5UInt64, hdf5Int32, hdf5Double }; -template -static void checkH5ErrImpl(TT status, const char *file, int line) { +template +static void checkH5ErrImpl(TT status, const char* file, int line) { if (status < 0) std::cerr << "An HDF5 error occurred (" << file << ": " << line << ")" << std::endl; } @@ -29,136 +26,141 @@ int encodeBoundary(std::array decoded) { } class Writer { -public: - explicit Writer(Mesh *mesh) : mesh(mesh) {}; + public: + explicit Writer(Mesh* mesh) : mesh(mesh) {}; - void writeHdf5(const std::string &filename) { - // Create file - hid_t h5falist = H5Pcreate(H5P_FILE_ACCESS); - checkH5Err(h5falist); + void writeHdf5(const std::string& filename) { + // Create file + hid_t h5falist = H5Pcreate(H5P_FILE_ACCESS); + checkH5Err(h5falist); #ifdef H5F_LIBVER_V18 - checkH5Err(H5Pset_libver_bounds(h5falist, H5F_LIBVER_V18, H5F_LIBVER_V18)); + checkH5Err(H5Pset_libver_bounds(h5falist, H5F_LIBVER_V18, H5F_LIBVER_V18)); #else - checkH5Err(H5Pset_libver_bounds(h5falist, H5F_LIBVER_LATEST, H5F_LIBVER_LATEST)); + checkH5Err(H5Pset_libver_bounds(h5falist, H5F_LIBVER_LATEST, H5F_LIBVER_LATEST)); #endif - checkH5Err(H5Pset_fapl_mpio(h5falist, MPI_COMM_WORLD, MPI_INFO_NULL)); - hid_t h5file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, h5falist); - checkH5Err(h5file); - checkH5Err(H5Pclose(h5falist)); - - writeData(h5file, "/connect", - std::vector{mesh->elementSize, 4}, - Hdf5DataType::hdf5UInt64, - mesh->connect.data()); - - writeData(h5file, "/geometry", - std::vector{mesh->vertexSize, 3}, - Hdf5DataType::hdf5Double, - mesh->vertices.data()); - - writeData(h5file, "/group", - std::vector({mesh->elementSize}), - Hdf5DataType::hdf5Int32, - mesh->elementGroups.data()); - - auto encodedBoundary = std::vector(mesh->elementSize); - for (auto i = 0u; i < mesh->elementBoundaries.size(); ++i) { - encodedBoundary[i] = encodeBoundary(mesh->elementBoundaries[i]); - } - writeData(h5file, "/boundary", - std::vector({mesh->elementSize}), - Hdf5DataType::hdf5Int32, - encodedBoundary.data()); + checkH5Err(H5Pset_fapl_mpio(h5falist, MPI_COMM_WORLD, MPI_INFO_NULL)); + hid_t h5file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, h5falist); + checkH5Err(h5file); + checkH5Err(H5Pclose(h5falist)); + + writeData(h5file, + "/connect", + std::vector{mesh->elementSize, 4}, + Hdf5DataType::hdf5UInt64, + mesh->connect.data()); + + writeData(h5file, + "/geometry", + std::vector{mesh->vertexSize, 3}, + Hdf5DataType::hdf5Double, + mesh->vertices.data()); + + writeData(h5file, + "/group", + std::vector({mesh->elementSize}), + Hdf5DataType::hdf5Int32, + mesh->elementGroups.data()); + + auto encodedBoundary = std::vector(mesh->elementSize); + for (auto i = 0u; i < mesh->elementBoundaries.size(); ++i) { + encodedBoundary[i] = encodeBoundary(mesh->elementBoundaries[i]); } + writeData(h5file, + "/boundary", + std::vector({mesh->elementSize}), + Hdf5DataType::hdf5Int32, + encodedBoundary.data()); + } - void writeXdmf(const std::string &xdmfFileName, const std::string &h5FileName) { - std::ofstream xdmf(xdmfFileName.c_str()); - - const auto globalSize = std::array{mesh->elementSize, mesh->vertexSize}; - xdmf << "" << std::endl - << "" << std::endl - << "" << std::endl - << " " << std::endl - << " " << std::endl - << " " - << std::endl - // This should be UInt but for some reason this does not work with - // binary data - << " " << h5FileName << ":/connect" << std::endl - << " " << std::endl - << " " << std::endl - << " " << h5FileName - << ":/geometry" << std::endl - << " " << std::endl - << " " << std::endl - << " " << h5FileName << ":/group" << std::endl - << " " << std::endl - << " " << std::endl - << " " << h5FileName << ":/boundary" << std::endl - << " " << std::endl - << " " << std::endl - << " " << std::endl - << "" << std::endl; - } + void writeXdmf(const std::string& xdmfFileName, const std::string& h5FileName) { + std::ofstream xdmf(xdmfFileName.c_str()); + + const auto globalSize = std::array{mesh->elementSize, mesh->vertexSize}; + xdmf << "" << std::endl + << "" << std::endl + << "" << std::endl + << " " << std::endl + << " " << std::endl + << " " + << std::endl + // This should be UInt but for some reason this does not work with + // binary data + << " " << h5FileName << ":/connect" << std::endl + << " " << std::endl + << " " << std::endl + << " " << h5FileName + << ":/geometry" << std::endl + << " " << std::endl + << " " << std::endl + << " " << h5FileName << ":/group" << std::endl + << " " << std::endl + << " " << std::endl + << " " << h5FileName << ":/boundary" << std::endl + << " " << std::endl + << " " << std::endl + << " " << std::endl + << "" << std::endl; + } -private: - Mesh *mesh; - - void writeData(hid_t h5file, - const std::string &name, - std::vector sizes, - Hdf5DataType type, - void *data) { - hid_t h5space = H5Screate_simple(sizes.size(), sizes.data(), nullptr); - checkH5Err(h5space); - - hid_t h5data; - switch (type) { - case Hdf5DataType::hdf5UInt64: - h5data = H5Dcreate(h5file, name.c_str(), H5T_STD_U64LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - break; - case Hdf5DataType::hdf5Int32: - h5data = H5Dcreate(h5file, name.c_str(), H5T_STD_I32LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - break; - case Hdf5DataType::hdf5Double: - h5data = H5Dcreate(h5file, name.c_str(), H5T_IEEE_F64LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - break; - } - checkH5Err(h5data); - - hsize_t start[2] = {0, 0}; - checkH5Err(H5Sselect_hyperslab(h5space, H5S_SELECT_SET, start, nullptr, sizes.data(), nullptr)); - - hid_t h5memspace = H5Screate_simple(sizes.size(), sizes.data(), nullptr); - checkH5Err(h5memspace); - - hid_t h5dxlist = H5Pcreate(H5P_DATASET_XFER); - checkH5Err(h5dxlist); - checkH5Err(H5Pset_dxpl_mpio(h5dxlist, H5FD_MPIO_COLLECTIVE)); - - switch (type) { - case Hdf5DataType::hdf5UInt64: - checkH5Err(H5Dwrite(h5data, H5T_NATIVE_ULONG, h5memspace, h5space, h5dxlist, data)); - break; - case Hdf5DataType::hdf5Int32: - checkH5Err(H5Dwrite(h5data, H5T_NATIVE_INT, h5memspace, h5space, h5dxlist, data)); - break; - case Hdf5DataType::hdf5Double: - checkH5Err(H5Dwrite(h5data, H5T_NATIVE_DOUBLE, h5memspace, h5space, h5dxlist, data)); - break; - } + private: + Mesh* mesh; + + void writeData(hid_t h5file, + const std::string& name, + std::vector sizes, + Hdf5DataType type, + void* data) { + hid_t h5space = H5Screate_simple(sizes.size(), sizes.data(), nullptr); + checkH5Err(h5space); + + hid_t h5data; + switch (type) { + case Hdf5DataType::hdf5UInt64: + h5data = H5Dcreate( + h5file, name.c_str(), H5T_STD_U64LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + break; + case Hdf5DataType::hdf5Int32: + h5data = H5Dcreate( + h5file, name.c_str(), H5T_STD_I32LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + break; + case Hdf5DataType::hdf5Double: + h5data = H5Dcreate( + h5file, name.c_str(), H5T_IEEE_F64LE, h5space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + break; } - + checkH5Err(h5data); + + hsize_t start[2] = {0, 0}; + checkH5Err(H5Sselect_hyperslab(h5space, H5S_SELECT_SET, start, nullptr, sizes.data(), nullptr)); + + hid_t h5memspace = H5Screate_simple(sizes.size(), sizes.data(), nullptr); + checkH5Err(h5memspace); + + hid_t h5dxlist = H5Pcreate(H5P_DATASET_XFER); + checkH5Err(h5dxlist); + checkH5Err(H5Pset_dxpl_mpio(h5dxlist, H5FD_MPIO_COLLECTIVE)); + + switch (type) { + case Hdf5DataType::hdf5UInt64: + checkH5Err(H5Dwrite(h5data, H5T_NATIVE_ULONG, h5memspace, h5space, h5dxlist, data)); + break; + case Hdf5DataType::hdf5Int32: + checkH5Err(H5Dwrite(h5data, H5T_NATIVE_INT, h5memspace, h5space, h5dxlist, data)); + break; + case Hdf5DataType::hdf5Double: + checkH5Err(H5Dwrite(h5data, H5T_NATIVE_DOUBLE, h5memspace, h5space, h5dxlist, data)); + break; + } + } }; - -#endif //REMOVE_WATERLAYER_WRITER_H +#endif // REMOVE_WATERLAYER_WRITER_H diff --git a/remove-groups/src/main.cpp b/remove-groups/src/main.cpp index d8bbeab..1ef5eba 100644 --- a/remove-groups/src/main.cpp +++ b/remove-groups/src/main.cpp @@ -1,19 +1,19 @@ -#include -#include -#include -#include -#include -#include +#include "Filter.h" +#include "Reader.h" +#include "Writer.h" +#include "mpi.h" + #include #include +#include #include +#include +#include +#include +#include +#include -#include "mpi.h" -#include "Reader.h" -#include "Filter.h" -#include "Writer.h" - -int main(int argc, char **argv) { +int main(int argc, char** argv) { MPI_Init(&argc, &argv); utils::Args args; args.addOption("input", 'i', "input mesh"); @@ -24,29 +24,25 @@ int main(int argc, char **argv) { return -1; } - auto groupsString = args.getArgument("remove-groups"); const std::regex groupsRegex("\\,"); std::set groupsToRemove{}; - std::for_each(std::sregex_token_iterator(groupsString.begin(), - groupsString.end(), - groupsRegex, - -1), - std::sregex_token_iterator(), - [&groupsToRemove](auto &arg) { - unsigned int groupInt; - auto[ptr, parseError] = std::from_chars( - arg.str().data(), - arg.str().data() + arg.str().size(), - groupInt); - if (parseError == std::errc::invalid_argument - || parseError == std::errc::result_out_of_range) { - std::cerr << "Group string is not a comma separated list of unsigned integers" << std::endl; - std::abort(); - } - groupsToRemove.insert(groupInt); - }); - + std::for_each( + std::sregex_token_iterator(groupsString.begin(), groupsString.end(), groupsRegex, -1), + std::sregex_token_iterator(), + [&groupsToRemove](auto& arg) { + unsigned int groupInt; + auto [ptr, parseError] = + std::from_chars(arg.str().data(), arg.str().data() + arg.str().size(), groupInt); + if (parseError == std::errc::invalid_argument || + parseError == std::errc::result_out_of_range) { + std::cerr << "Group string is not a comma separated list " + "of unsigned integers" + << std::endl; + std::abort(); + } + groupsToRemove.insert(groupInt); + }); // Input const auto meshFile = args.getArgument("input"); diff --git a/submodules/__init__.py b/submodules/__init__.py index 8b13789..e69de29 100644 --- a/submodules/__init__.py +++ b/submodules/__init__.py @@ -1 +0,0 @@ - diff --git a/test_communication_pattern/SConstruct b/test_communication_pattern/SConstruct index 53d788d..ed8d9d3 100644 --- a/test_communication_pattern/SConstruct +++ b/test_communication_pattern/SConstruct @@ -1,4 +1,3 @@ -#! /usr/bin/python ## # @file # This file is part of SeisSol. @@ -8,21 +7,21 @@ # @section LICENSE # Copyright (c) 2017, SeisSol Group # All rights reserved. -# +# # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: -# +# # 1. Redistributions of source code must retain the above copyright notice, # this list of conditions and the following disclaimer. -# +# # 2. Redistributions in binary form must reproduce the above copyright notice, # this list of conditions and the following disclaimer in the documentation # and/or other materials provided with the distribution. -# +# # 3. Neither the name of the copyright holder nor the names of its # contributors may be used to endorse or promote products derived from this # software without specific prior written permission. -# +# # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE diff --git a/test_communication_pattern/src/Comm.cpp b/test_communication_pattern/src/Comm.cpp index 0494e99..eeed998 100644 --- a/test_communication_pattern/src/Comm.cpp +++ b/test_communication_pattern/src/Comm.cpp @@ -1,46 +1,47 @@ -#include #include "Comm.h" -#include +#include #include +#include #include -#include -#include #include +#include +#include inline double usec(struct timeval start, struct timeval end) { - return ((double)(((end.tv_sec * 1000000 + end.tv_usec) - (start.tv_sec * 1000000 + start.tv_usec)))); + return ( + (double)(((end.tv_sec * 1000000 + end.tv_usec) - (start.tv_sec * 1000000 + start.tv_usec)))); } -void testCommunication(counter_t* edgeCut, size_t messageSize, int numTimesteps) -{ +void testCommunication(counter_t* edgeCut, size_t messageSize, int numTimesteps) { struct timeval start_time, end_time, start_time_total, end_time_total; - double time_min = std::numeric_limits::max(), time_max = std::numeric_limits::min(), time_avg; + double time_min = std::numeric_limits::max(), + time_max = std::numeric_limits::min(), time_avg; int rank, numRanks; MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &numRanks); - + unsigned N = 0; counter_t E = 0; for (unsigned r = 0; r < numRanks; ++r) { N += (edgeCut[r] > 0) ? 1 : 0; E += edgeCut[r]; } - - MPI_Request* requests = new MPI_Request[2*N]; - unsigned char* copy = new unsigned char[E*messageSize]; - unsigned char* ghost = new unsigned char[E*messageSize]; - for (unsigned i = 0; i < E*messageSize; ++i) { + + MPI_Request* requests = new MPI_Request[2 * N]; + unsigned char* copy = new unsigned char[E * messageSize]; + unsigned char* ghost = new unsigned char[E * messageSize]; + for (unsigned i = 0; i < E * messageSize; ++i) { copy[i] = static_cast(lrand48() % 256); ghost[i] = static_cast(lrand48() % 256); } - + if (rank == 0) { printf("Rank tavg tmin tmax bw\n"); } MPI_Barrier(MPI_COMM_WORLD); - + gettimeofday(&start_time_total, NULL); for (unsigned t = 0; t < numTimesteps; ++t) { unsigned e = 0; @@ -48,13 +49,25 @@ void testCommunication(counter_t* edgeCut, size_t messageSize, int numTimesteps) gettimeofday(&start_time, NULL); for (unsigned r = 0; r < numRanks; ++r) { if (edgeCut[r] > 0) { - MPI_Isend(copy + e*messageSize, edgeCut[r]*messageSize, MPI_BYTE, r, 0, MPI_COMM_WORLD, requests + 2*n); - MPI_Irecv(ghost + e*messageSize, edgeCut[r]*messageSize, MPI_BYTE, r, 0, MPI_COMM_WORLD, requests + 2*n+1); + MPI_Isend(copy + e * messageSize, + edgeCut[r] * messageSize, + MPI_BYTE, + r, + 0, + MPI_COMM_WORLD, + requests + 2 * n); + MPI_Irecv(ghost + e * messageSize, + edgeCut[r] * messageSize, + MPI_BYTE, + r, + 0, + MPI_COMM_WORLD, + requests + 2 * n + 1); e += edgeCut[r]; ++n; } } - MPI_Waitall(2*N, requests, MPI_STATUSES_IGNORE); + MPI_Waitall(2 * N, requests, MPI_STATUSES_IGNORE); gettimeofday(&end_time, NULL); double time = usec(start_time, end_time); time_min = std::min(time_min, time); @@ -62,8 +75,13 @@ void testCommunication(counter_t* edgeCut, size_t messageSize, int numTimesteps) } gettimeofday(&end_time_total, NULL); time_avg = usec(start_time_total, end_time_total) / numTimesteps; - - printf("%4d %10.2lf %10.2lf %10.2lf %10.2lf\n", rank, time_avg, time_min, time_max, E*messageSize / 1.048576 / time_avg); - - delete[] requests; + + printf("%4d %10.2lf %10.2lf %10.2lf %10.2lf\n", + rank, + time_avg, + time_min, + time_max, + E * messageSize / 1.048576 / time_avg); + + delete[] requests; } diff --git a/test_communication_pattern/src/SConscript b/test_communication_pattern/src/SConscript index 6fec9db..b151869 100644 --- a/test_communication_pattern/src/SConscript +++ b/test_communication_pattern/src/SConscript @@ -1,4 +1,3 @@ -#!/usr/bin/python ## # @file # This file is part of SeisSol. @@ -8,21 +7,21 @@ # @section LICENSE # Copyright (c) 2017, SeisSol Group # All rights reserved. -# +# # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: -# +# # 1. Redistributions of source code must retain the above copyright notice, # this list of conditions and the following disclaimer. -# +# # 2. Redistributions in binary form must reproduce the above copyright notice, # this list of conditions and the following disclaimer in the documentation # and/or other materials provided with the distribution. -# +# # 3. Neither the name of the copyright holder nor the names of its # contributors may be used to endorse or promote products derived from this # software without specific prior written permission. -# +# # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE diff --git a/test_communication_pattern/src/main.cpp b/test_communication_pattern/src/main.cpp index febca32..652db11 100644 --- a/test_communication_pattern/src/main.cpp +++ b/test_communication_pattern/src/main.cpp @@ -1,31 +1,31 @@ -#include -#include +#include "Comm.h" + #include +#include #include +#include #include #include -#include "Comm.h" - -counter_t* readCommMatrix(std::string const& matrixFile) -{ +counter_t* readCommMatrix(const std::string& matrixFile) { int rank, numRanks; MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &numRanks); - + counter_t* edges = NULL; FILE* file = fopen(matrixFile.c_str(), "rb"); fseek(file, 0, SEEK_END); long int size = ftell(file); - + int N = static_cast(sqrt(size / sizeof(counter_t))); - if (N*N*sizeof(counter_t) != size) { + if (N * N * sizeof(counter_t) != size) { logError() << "The communication matrix is not a square matrix."; return NULL; - } + } if (N != numRanks) { - logError() << "Communication matrix has" << N << "ranks whereas this test was started with" << numRanks << "ranks."; + logError() << "Communication matrix has" << N << "ranks whereas this test was started with" + << numRanks << "ranks."; return NULL; } @@ -33,14 +33,13 @@ counter_t* readCommMatrix(std::string const& matrixFile) edges = new counter_t[N]; fread(edges, sizeof(counter_t), N, file); fclose(file); - + return edges; } -int main(int argc, char** argv) -{ +int main(int argc, char** argv) { MPI_Init(&argc, &argv); - + utils::Args args; args.addOption("order", 'o', "Convergence order.", utils::Args::Required, false); args.addOption("num-quantities", 'q', "Number of quantities.", utils::Args::Required, false); @@ -53,15 +52,15 @@ int main(int argc, char** argv) int order = args.getArgument("order", 6); int nq = args.getArgument("num-quantities", 9); std::string matrixFile = args.getAdditionalArgument("matrix"); - + counter_t* edges = readCommMatrix(matrixFile); - - size_t messageSize = sizeof(double) * nq * order * (order+1) * (order+2) / 6; + + size_t messageSize = sizeof(double) * nq * order * (order + 1) * (order + 2) / 6; testCommunication(edges, messageSize); - - delete[] edges; - + + delete[] edges; + MPI_Finalize(); - + return 0; }