diff --git a/cmake/FindNetCDF.cmake b/cmake/FindNetCDF.cmake new file mode 100644 index 0000000..c361d09 --- /dev/null +++ b/cmake/FindNetCDF.cmake @@ -0,0 +1,137 @@ +#[==[ +Provides the following variables: + + * `NetCDF_FOUND`: Whether NetCDF was found or not. + * `NetCDF_INCLUDE_DIRS`: Include directories necessary to use NetCDF. + * `NetCDF_LIBRARIES`: Libraries necessary to use NetCDF. + * `NetCDF_VERSION`: The version of NetCDF found. + * `NetCDF::NetCDF`: A target to use with `target_link_libraries`. +#]==] + +# Stolen from VTK with some small modifications +# https://github.com/Kitware/VTK/blob/master/CMake/FindNetCDF.cmake +# (BSD 3-Clause License) + +#[[========================================================================= + + Program: Visualization Toolkit + Module: Copyright.txt + +Copyright (c) 1993-2015 Ken Martin, Will Schroeder, Bill Lorensen +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + * 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. + + * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names + of any 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 +ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +=========================================================================]] +# Note: The version from VTK tries to first find a CMake built version of CMake and then +# a pkgConf built one. This doesn't work on all systems with custom built versions of +# CMake + HDF5. +# Thus, first try to find it manually, then fall back on CMake/PkgConfig + +find_path(NetCDF_INCLUDE_DIR + NAMES netcdf.h + HINTS $ENV{NETCDF_DIR}/include $ENV{NETCDF_INCLUDE_DIR} + DOC "netcdf include directories") +mark_as_advanced(NetCDF_INCLUDE_DIR) + +find_library(NetCDF_LIBRARY + NAMES netcdf + HINTS $ENV{NETCDF_DIR}/lib $ENV{NETCDF_LIB_DIR} + DOC "netcdf library") +mark_as_advanced(NetCDF_LIBRARY) + +string(REGEX REPLACE "(.*)libnetcdf\.(so|a)" "\\1" NetCDF_LIBRARY_DIR "${NetCDF_LIBRARY}") +mark_as_advanced(NetCDF_LIBRARY_DIR) + +if (NetCDF_INCLUDE_DIR) + file(STRINGS "${NetCDF_INCLUDE_DIR}/netcdf_meta.h" _netcdf_version_lines + REGEX "#define[ \t]+NC_VERSION_(MAJOR|MINOR|PATCH|NOTE)") + string(REGEX REPLACE ".*NC_VERSION_MAJOR *\([0-9]*\).*" "\\1" _netcdf_version_major "${_netcdf_version_lines}") + string(REGEX REPLACE ".*NC_VERSION_MINOR *\([0-9]*\).*" "\\1" _netcdf_version_minor "${_netcdf_version_lines}") + string(REGEX REPLACE ".*NC_VERSION_PATCH *\([0-9]*\).*" "\\1" _netcdf_version_patch "${_netcdf_version_lines}") + string(REGEX REPLACE ".*NC_VERSION_NOTE *\"\([^\"]*\)\".*" "\\1" _netcdf_version_note "${_netcdf_version_lines}") + set(NetCDF_VERSION "${_netcdf_version_major}.${_netcdf_version_minor}.${_netcdf_version_patch}${_netcdf_version_note}") + unset(_netcdf_version_major) + unset(_netcdf_version_minor) + unset(_netcdf_version_patch) + unset(_netcdf_version_note) + unset(_netcdf_version_lines) +endif () + +if (NetCDF_FOUND) + set(NetCDF_INCLUDE_DIRS "${NetCDF_INCLUDE_DIR}") + set(NetCDF_LIBRARIES "${NetCDF_LIBRARY}") + + if (NOT TARGET NetCDF::NetCDF) + add_library(NetCDF::NetCDF UNKNOWN IMPORTED) + set_target_properties(NetCDF::NetCDF PROPERTIES + IMPORTED_LOCATION "${NetCDF_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${NetCDF_INCLUDE_DIR}") + endif () +else() + # Try to find a CMake-built NetCDF. + find_package(netCDF CONFIG QUIET) + if (netCDF_FOUND) + # Forward the variables in a consistent way. + set(NetCDF_FOUND "${netCDF_FOUND}") + set(NetCDF_INCLUDE_DIRS "${netCDF_INCLUDE_DIR}") + set(NetCDF_LIBRARIES "${netCDF_LIBRARIES}") + set(NetCDF_VERSION "${NetCDFVersion}") + if (NOT TARGET NetCDF::NetCDF) + add_library(NetCDF::NetCDF INTERFACE IMPORTED) + set_target_properties(NetCDF::NetCDF PROPERTIES + INTERFACE_LINK_LIBRARIES "${NetCDF_LIBRARIES}") + endif () + # Skip the rest of the logic in this file. + return () + endif () + + find_package(PkgConfig QUIET) + if (PkgConfig_FOUND) + pkg_check_modules(_NetCDF QUIET netcdf IMPORTED_TARGET) + if (_NetCDF_FOUND) + # Forward the variables in a consistent way. + set(NetCDF_FOUND "${_NetCDF_FOUND}") + set(NetCDF_INCLUDE_DIRS "${_NetCDF_INCLUDE_DIRS}") + set(NetCDF_LIBRARIES "${_NetCDF_LIBRARIES}") + set(NetCDF_VERSION "${_NetCDF_VERSION}") + if (NOT TARGET NetCDF::NetCDF) + add_library(NetCDF::NetCDF INTERFACE IMPORTED) + set_target_properties(NetCDF::NetCDF PROPERTIES + INTERFACE_LINK_LIBRARIES "PkgConfig::_NetCDF") + endif () + # Skip the rest of the logic in this file. + return () + endif () + endif () +endif () + + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(NetCDF + REQUIRED_VARS NetCDF_LIBRARY NetCDF_INCLUDE_DIR + VERSION_VAR NetCDF_VERSION) + diff --git a/cmake/README.md b/cmake/README.md new file mode 100644 index 0000000..a49b1e5 --- /dev/null +++ b/cmake/README.md @@ -0,0 +1 @@ +A selective copy of the CMake files from the SeisSol main repository. diff --git a/cube/cubeGenerator.py b/cube/cubeGenerator.py new file mode 100755 index 0000000..8d010e7 --- /dev/null +++ b/cube/cubeGenerator.py @@ -0,0 +1,68 @@ +#!/usr/bin/python +## +# @file +# This file is part of SeisSol. +# +# @author Sebastian Rettenberger (rettenbs AT in.tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger,_M.Sc.) +# +# @section LICENSE +# Copyright (c) 2013, 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 +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +from lib.args import Args +from lib.gambit import GambitWriter +from lib.mesh import Mesh + +import operator + +def main(): + # Parse command line arguements + args = Args() + + print 'Generating mesh with size', ', '.join(map(str, args.size())) + + # Create Mesh + mesh = Mesh(args.size(), args.boundary()) + + print 'Number of elements:', len(mesh.elements()) + + # Write mesh + if args.netcdf(): + try: + from lib.netcdf import NetcdfWriter + except ImportError, e: + print 'netcdf4-python could not be loaded:', e + return + print 'Mesh will contain', ' * '.join(map(str, args.partitions())), '=', reduce(operator.mul, args.partitions(), 1), 'partitions' + NetcdfWriter(mesh, args.partitions(), args.outputFile()) + else: + GambitWriter(mesh, args.outputFile()) + +if __name__ == '__main__': + main() diff --git a/cube/lib/__init__.py b/cube/lib/__init__.py new file mode 100644 index 0000000..4464e3f --- /dev/null +++ b/cube/lib/__init__.py @@ -0,0 +1,36 @@ +#!/usr/bin/python +## +# @file +# This file is part of SeisSol. +# +# @author Sebastian Rettenberger (rettenbs AT in.tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger,_M.Sc.) +# +# @section LICENSE +# Copyright (c) 2013, 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 +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. diff --git a/cube/lib/args.py b/cube/lib/args.py new file mode 100644 index 0000000..bea797a --- /dev/null +++ b/cube/lib/args.py @@ -0,0 +1,111 @@ +#! /usr/bin/python +## +# @file +# This file is part of SeisSol. +# +# @author Sebastian Rettenberger (rettenbs AT in.tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger,_M.Sc.) +# +# @section LICENSE +# Copyright (c) 2013, 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 +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +import argparse +import operator +import sys + +class Args: + """Parses the command line arguments using + argparse.ArgumentParser""" + + def __init__(self): + parser = argparse.ArgumentParser(prog='cubeGenerator') + parser.add_argument('-s', '--size', type=int, + help='number of cubes in each dimensions') + parser.add_argument('-x', '--size-x', type=int, + help='number of cubes in x dimension') + parser.add_argument('-y', '--size-y', type=int, + help='number of cubes in y dimension') + parser.add_argument('-z', '--size-z', type=int, + help='number of cubes in z dimension') + parser.add_argument('--px', type=int, + help='number of partitions x dimension') + parser.add_argument('--py', type=int, + help='number of partitions in y dimension') + parser.add_argument('--pz', type=int, + help='number of partitions in z dimension') + parser.add_argument('-o', '--output', type=argparse.FileType('w'), + required=True, help='output file for resulting Gambit or netCDF mesh') + parser.add_argument('-n', '--netcdf', action='store_true', + help='Create netCDF files'), + parser.add_argument('-b', '--boundary', type=int, + help='boundary condition (default 106 = periodic)', default=106) + parser.add_argument('-v', '--version', action='version', version='%(prog)s 0.1') + + try: + # Parse cmd line options + self.__options = parser.parse_args(sys.argv[1:]) + if self.__options.size == None \ + and (self.__options.size_x == None \ + or self.__options.size_y == None \ + or self.__options.size_z == None): + raise IOError('Either size or x/y/z is required') + + if self.__options.netcdf: + partitions = [getattr(self.__options, 'p'+i) for i in ['x', 'y', 'z']] + p = sum(p != None for p in partitions) + if p != 3: + raise IOError('Number of partitions is required for netCDF files') + + except IOError, e: + parser.error(str(e)) + + # Set correct size options + self.__options.size = [self.__options.size for _ in range(3)] + + if self.__options.size_x != None: + self.__options.size[0] = self.__options.size_x + if self.__options.size_y != None: + self.__options.size[1] = self.__options.size_y + if self.__options.size_z != None: + self.__options.size[2] = self.__options.size_z + + def size(self): + return self.__options.size + + def outputFile(self): + return self.__options.output + + def boundary(self): + return self.__options.boundary + + def netcdf(self): + return self.__options.netcdf + + def partitions(self): + return (self.__options.px, self.__options.py, self.__options.pz) diff --git a/cube/lib/gambit.py b/cube/lib/gambit.py new file mode 100644 index 0000000..21efbe7 --- /dev/null +++ b/cube/lib/gambit.py @@ -0,0 +1,442 @@ +#!/usr/bin/python +## +# @file +# This file is part of SeisSol. +# +# @author Sebastian Rettenberger (rettenbs AT in.tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger,_M.Sc.) +# +# @section LICENSE +# Copyright (c) 2013, 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 +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +import collections +import math +import re + +FLOAT_REGEX_STRING = '[+\-0-9\.eE]' +ENDSECTION = 'ENDOFSECTION' +ENDSECTION_REGEX = re.compile('^ENDOFSECTION$') + +class GambitFile(collections.Iterable): + """Provides an file like iterator but allows seek/tell and removes comments/empty lines""" + + def __init__(self, file): + self.__file = file + + def __del__(self): + self.__file.close() + + def __iter__(self): + line = self.__file.readline() + while line: + line = GambitFile.__prepareLine(line) + if line == False: + continue + + yield line + + line = self.__file.readline() + + def seek(self, pos): + self.__file.seek(pos) + + def tell(self): + return self.__file.tell() + + def skipSection(self): + for line in self: + # End of section? + if line == ENDSECTION: + break + + @staticmethod + def __prepareLine(line): + # Comment? + if line[0] == '/': + return False + + return line.strip() + +class GambitReader(object): + """Reads a gambit neutral file""" + + class Cells(collections.Iterable): + def __init__(self, len, file): + self.__len = len + self.__file = file + self.__seek = file.tell() + + # Iterator finished? + self.__finished = True + + # Skip cells now, so the reader can continue + file.skipSection() + + def __len__(self): + return self.__len + + def __iter__(self): + file = self.__file + + # Reset iterator + if self.__finished: + self.__finished = False + file.seek(self.__seek) + + # This regex only supports tetrahedra + tetCell = re.compile('^('+FLOAT_REGEX_STRING+'+) +6 +4 +(\d+) +(\d+) +(\d+) +(\d+)$') + + nextCellNumber = 1 + for line in file: + # Check for cell + match = re.match(tetCell, line) + if match: + # Check for correct cell number + if float(match.group(1)) != nextCellNumber: + raise IOError('Cells are not sorted, wrong cell number after '+str(nextCellNumber)) + nextCellNumber += 1 + + yield [int(match.group(x))-1 for x in [2, 3, 4, 5]] + continue + + # End of section? + if line == ENDSECTION: + self.__finished = True + return + + raise IOError('Error in cell section in line: '+line) + + + def __init__(self, file, sparse = False): + """sparse currently does not support coords/groups/boundaries""" + + # Open file if a string is specified + if isinstance(file, basestring): + file = open(file, 'rU') + + # We need to close it later + self.__file = GambitFile(file) + + # Read header + self.__readHeader(self.__file) + + # Create array for groups and boundaries + self.__groups = [] + self.__boundaries = [] + + # Read section by section + for line in self.__file: + if line == 'NODAL COORDINATES '+self.version: + self.__readCoords(self.__file, not sparse) + elif line == 'ELEMENTS/CELLS '+self.version: + self.__readCells(self.__file, sparse) + elif line == 'ELEMENT GROUP '+self.version: + self.__readGroup(self.__file, not sparse) + elif line == 'BOUNDARY CONDITIONS '+self.version: + self.__readBoundary(self.__file, not sparse) + else: + raise IOError('Could not parse line: '+line) + + # Check header sizes + if self.problemSize[0] != len(self.__coords): + raise IOError('Number of coordinates not correct: number in header = ' + +str(self.problemSize[0])+' != actual number = '+str(len(self.__coords))) + + if self.problemSize[1] != len(self.__cells): + raise IOError('Number of cells not correct: number in header = ' + +str(self.problemSize[1])+' != actual number = '+str(len(self.__cells))) + + if self.problemSize[2] != len(self.__groups): + raise IOError('Number of groups not correct: number in header = ' + +str(self.problemSize[2])+' != actual number = '+str(len(self.__groups))) + + if self.problemSize[3] != len(self.__boundaries): + raise IOError('Number of boundaries not correct: number in header = ' + +str(self.problemSize[3])+' != actual number = '+str(len(self.__boundaries))) + + def __del__(self): + del self.__file + + def __readHeader(self, file): + header = collections.deque([ + (re.compile('^CONTROL INFO (\d+\.\d+\.\d+)$'), 'version', lambda m: m.group(1)), + (re.compile('^\*\* GAMBIT NEUTRAL FILE$'), None), + (re.compile('^(.*)$'), 'name', lambda m: m.group(1)), + (re.compile('^PROGRAM: .{20} +(VERSION: +)?\d+\.\d+\.\d+$'), None), + (re.compile('^.*$'), 'date', lambda m: m.group(0)), + (re.compile('^NUMNP {5}NELEM {5}NGRPS {4}NBSETS {5}NDFCD {5}NDFVL$'), None), + (re.compile('^(\d+) +(\d+) +(\d+) +(\d+) +(\d+) +(\d+)$'), + 'problemSize', lambda m: [int(m.group(x+1)) for x in range(6)]), + (ENDSECTION_REGEX, None) + ]) + + for line in file: + # Check if we found the correct header + match = re.match(header[0][0], line) + if not match: + raise IOError('Could not parse header line: '+line) + + # Save matched expression? + if header[0][1]: + self.__setattr__(header[0][1], header[0][2](match)) + + # Header found -> goto next header + header.popleft() + + # All headers -> stop + if not header: + break + + def __readCoords(self, file, save): + if save: + coord = re.compile('^('+FLOAT_REGEX_STRING+'+) +('+FLOAT_REGEX_STRING+'+) +(' + +FLOAT_REGEX_STRING+'+) +('+FLOAT_REGEX_STRING+'+)$') + + self.__coords = [] + + nextCoordNumber = 1 + for line in file: + # Check for a coordinate + match = re.match(coord, line) + if match: + # Check for correct coordinate number + if float(match.group(1)) != nextCoordNumber: + raise IOError('Coordinates are not sorted, wrong coordinate number after '+str(nextCoordNumber)) + nextCoordNumber += 1 + + self.__coords.append([float(match.group(x)) for x in [2, 3, 4]]) + continue + + # End of section? + if line == ENDSECTION: + break + + raise IOError('Error in coordinate section in line: '+line) + else: + self.__coords = xrange(self.problemSize[0]) + file.skipSection() + + def __readCells(self, file, sparse): + if sparse: + self.__cells = GambitReader.Cells(self.problemSize[1], file) + return + + self.__cells = [c for c in GambitReader.Cells(self.problemSize[1], file)] + + def __readGroup(self, file, save): + if save: + header = collections.deque([ + (re.compile('^GROUP: +\d+ ELEMENTS: +(\d+) ?MATERIAL: +(\d+) NFLAGS: +\d+$'), + 'header', lambda m: [int(m.group(i)) for i in [1, 2]]), + (re.compile('^.*$'), 'name', lambda m: m.group(0)), + (re.compile('^.*$'), None) # Ignore flags + ]) + + # New group + group = dict() + + # Read headers + for line in file: + # Check if we found the correct header + match = re.match(header[0][0], line) + if not match: + raise IOError('Could not parse group header line: '+line) + + # Save matched expression? + if header[0][1]: + group[header[0][1]] = header[0][2](match) + + # Header found -> goto next header + header.popleft() + + # All headers -> stop + if not header: + break + + # Split header information + group['size'] = group['header'][0] + group['material'] = group['header'][1] + del group['header'] + + # Add cell list + group['cells'] = [] + + # Read cells + for line in file: + # End of section? + if line == ENDSECTION: + break + + cells = line.split() + if len(cells) > 10: + raise IOError('Too many cells in one line in group section: '+line) + + group['cells'].extend(map(lambda x: int(x)-1, cells)) + + # Check group size + if group['size'] != len(group['cells']): + raise IOError('Size of group '+group['name']+' not correct: size in header = ' + +str(group['size'])+' != actual size = '+str(len(group['cells']))) + + # Add group + self.__groups.append(group) + + else: + # Add dummy group + self.__groups.append(None) + file.skipSection() + + def __readBoundary(self, file, save): + if save: + # Very strict at the moment, we might need to relax this + header = re.compile('^(\d+) +1 +(\d+) +0 +6$') + value = re.compile('^(\d+) +6 +(\d+)') + + # New boundary + boundary = dict() + + # Read header + for line in file: + # Check for correct header + match = re.match(header, line) + if not match: + raise IOError('Could not parse boundary header line: '+line) + + # Get values + boundary['name'] = match.group(1) + boundary['size'] = int(match.group(2)) + + # We only have one header line + break + + boundary['sides'] = [] + + # Read values + for line in file: + # Check for boundary + match = re.match(value, line) + if match: + boundary['sides'].append([int(match.group(1))-1, int(match.group(2))]) + continue + + # End of section? + if line == ENDSECTION: + break; + + raise IOError('Error in cell section in line: '+line) + + # Check size + if boundary['size'] != len(boundary['sides']): + raise IOError('Size of boundary '+boundary['name']+' not correct: size in header = ' + +str(boundary['size'])+' != actual size = '+str(len(boundary['sides']))) + + # Add new boundary type + self.__boundaries.append(boundary) + + else: + # Add dummy boundary + self.__boundaries.append(None) + file.skipSection() + + def coords(self): + return self.__coords + + def elements(self): + return self.__cells + + def groups(self): + return self.__groups + + def boundaries(self): + return self.__boundaries + +class GambitWriter: + """Writes a GAMBIT mesh""" + + def __init__(self, mesh, file): + # Open file if a string is specified + if isinstance(file, basestring): + file = open(file, 'w') + + minElementIdSize = int(math.floor(math.log10(len(mesh.elements())))) + 1 + minVertexIdSize = int(math.floor(math.log10(len(mesh.coords())))) + 1 + + # Header + print >> file, '%20s' % ('CONTROL INFO'), mesh.version + print >> file, '** GAMBIT NEUTRAL FILE' + print >> file, mesh.name + print >> file, 'PROGRAM: %20s VERSION: %s' % ('Gambit', mesh.version) + print >> file, ' ' + mesh.date + print >> file, ' %9s %9s %9s %9s %9s %9s' % ('NUMNP', 'NELEM', 'NGRPS', 'NBSETS', 'NDFCD', 'NDFVL') + print >> file, ' %9d %9d %9d %9d %9d %9d' % (len(mesh.coords()), len(mesh.elements()), + len(mesh.groups()), len(mesh.boundaries()), 3, 3) + print >> file, ENDSECTION + + # Coords + print >> file, '%20s' % ('NODAL COORDINATES'), mesh.version + for i, coord in enumerate(mesh.coords()): + print >> file, '%10d %19.10e %19.10e %19.10e' % (i+1, coord[0], coord[1], coord[2]) + print >> file, ENDSECTION + + # Elements + elementIdSize = str(max(7, minElementIdSize)) + vertexIdSize = str(max(7, minVertexIdSize)) + print >> file, '%20s' % ('ELEMENTS/CELLS'), mesh.version + for i, element in enumerate(mesh.elements()): + print >> file, (' %'+elementIdSize+'d %2d %2d' \ + + ' %'+vertexIdSize+'d %'+vertexIdSize+'d %'+vertexIdSize+'d %'+vertexIdSize+'d') \ + % ((i+1, 6, 4) + tuple([x+1 for x in element])) + print >> file, ENDSECTION + + # Groups + elementIdSize = str(max(7, minElementIdSize+1)) + for i, group in enumerate(mesh.groups()): + print >> file, '%20s' % ('ELEMENT GROUP'), mesh.version + print >> file, 'GROUP: %10d ELEMENTS: %10d MATERIAL: %10d NFLAGS: %10d' % (i+1, group['size'], group['material'], 1) + print >> file, '%32s' % (group['name']) + print >> file, '%8d' % (0) + for i, cell in enumerate(group['cells']): + # Whitespace at beginning of each line + if i % 10 == 0: + print >> file, '', + print >> file, ('%'+elementIdSize+'d') % (cell+1), + # Newline after 10 elements or when finished + if (i+1) % 10 == 0 or i+1 == len(group['cells']): + print >> file + print >> file, ENDSECTION + + # Boundary conditions + for boundary in mesh.boundaries(): + print >> file, '%20s' % ('BOUNDARY CONDITIONS'), mesh.version + print >> file, '%32s%8d%8d%8d%8d' % (boundary['name'], 1, boundary['size'], 0, 6) + for side in boundary['sides']: + print >> file, '%10d %4d %4d' % (side[0]+1, 6, side[1]) + print >> file, ENDSECTION + + file.close() diff --git a/cube/lib/mesh.py b/cube/lib/mesh.py new file mode 100644 index 0000000..21aa35c --- /dev/null +++ b/cube/lib/mesh.py @@ -0,0 +1,341 @@ +#!/usr/bin/python +## +# @file +# This file is part of SeisSol. +# +# @author Sebastian Rettenberger (rettenbs AT in.tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger,_M.Sc.) +# +# @section LICENSE +# Copyright (c) 2013, 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 +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. + +import collections +import datetime +import numpy + +class Coords(collections.Sequence): + def __init__(self, mesh): + self.__mesh = mesh + self.__len1 = [2*i + 1 for i in self.__mesh.size()] # number of coords in each dimension + self.__len = reduce(lambda x, y: x*y, self.__len1) + + def len1(self): + return self.__len1 + + def __len__(self): + return self.__len + + def __getitem__(self, index): + if index >= self.__len: + raise StopIteration + + def scale(x, dim): + return (x - (self.__len1[dim]-1)/2.) * 5. / self.__mesh.size()[dim] + + x = scale(index % self.__len1[0], 0) + y = scale((index / self.__len1[0]) % self.__len1[1], 1) + z = scale(index / (self.__len1[0] * self.__len1[1]), 2) + + return (x, y, z) + +class Elements(collections.Sequence): + def __init__(self, mesh): + self.__mesh = mesh + self.__len = 40 * reduce(lambda x, y: x*y, self.__mesh.size()) + self.__len1 = [2 * i for i in self.__mesh.size()] # number of cubes in one dimension + + self.__checkOrientation() + + def __checkOrientation(self): + """Checks the correct orientation for the different tetrahedra""" + + for i in range(10): + vertices = self[i] + + coords = [None]*4 + for j in range(4): + coords[j] = numpy.array(self.__mesh.coords()[vertices[j]]) + + n = numpy.cross(-coords[0]+coords[1], -coords[0]+coords[2]) + orientation = numpy.dot(n, -coords[0]+coords[3]) + + if orientation <= 0: + raise Exception('Wrong orientation in element '+str(i)) + + def __len__(self): + return self.__len + + def tetUnitCube(self, index, even): + """Returns the coordinates of one of the 5 tetrahedron in an even or odd unit cube""" + + if even: + if index == 0: + return ((0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)) + if index == 1: + return ((1, 0, 0), (0, 1, 0), (1, 1, 1), (1, 1, 0)) + if index == 2: + return ((1, 0, 0), (1, 1, 1), (0, 0, 1), (1, 0, 1)) + if index == 3: + return ((0, 1, 0), (0, 1, 1), (0, 0, 1), (1, 1, 1)) + # index == 4 + return ((1, 0, 0), (0, 1, 0), (0, 0, 1), (1, 1, 1)) + + if index == 0: + return ((0, 0, 0), (0, 1, 0), (0, 1, 1), (1, 1, 0)) + if index == 1: + return ((0, 0, 0), (1, 1, 0), (1, 0, 1), (1, 0, 0)) + if index == 2: + return ((0, 0, 0), (1, 0, 1), (0, 1, 1), (0, 0, 1)) + if index == 3: + return ((1, 1, 0), (1, 0, 1), (1, 1, 1), (0, 1, 1)) + # index == 4 + return ((0, 0, 0), (1, 1, 0), (0, 1, 1), (1, 0, 1)) + + def __getitem__(self, index): + if index >= self.__len: + raise StopIteration + + # Index of the cube we are currently working on + cIndex = index/5 + + cx = cIndex % self.__len1[0] + cy = (cIndex / self.__len1[0]) % self.__len1[1] + cz = cIndex / (self.__len1[0] * self.__len1[1]) + + # The index inside the cube + i = index % 5 + + # Odd cube? + odd = (cx+cy+cz) % 2 + + # Number of coords in one dimension + coordLength = self.__mesh.coords().len1() + + def crd((x, y, z)): + return x+cx + ((y+cy) + (z+cz) * coordLength[1]) * coordLength[0] + + return map(crd, self.tetUnitCube(i, not odd)) + +class Group: + def __init__(self, name, material, mesh): + self._name = name + self._size = len(mesh.elements()) + self._material = material + self._cells = xrange(self._size) + + def __getitem__(self, key): + return getattr(self, '_'+str(key)) + +class Boundary: + class Sides(collections.Sequence): + """Generates the boundary condition for all face elements""" + + def __init__(self, mesh): + self.__len1 = [8 * i for i in # Number of faces on one side + [mesh.size()[0]*mesh.size()[1], # Top / Bottom + mesh.size()[0]*mesh.size()[2], # Left / Right behind + mesh.size()[1]*mesh.size()[2]]] # Right / Left behind + self.__len = 2 * reduce(lambda x, y: x+y, self.__len1) + self.__mesh = mesh + + def face(points): + """Returns the number of the face or None if its not a face""" + + if points == [0, 0, 0, 1]: + return 1 + if points == [0, 0, 1, 0]: + return 2 + if points == [0, 1, 0, 0]: + return 4 + if points == [1, 0, 0, 0]: + return 3 + + return None + + self.__top = [] + self.__left = [] + self.__right = [] + self.__bot = [] + self.__rightb = [] + self.__leftb = [] + + for i in range(5): + element1 = mesh.elements().tetUnitCube(i, 1) + + f = face(map(lambda e: e[0], element1)) + if f: + self.__right.append((i, f)) + f = face(map(lambda e: e[1], element1)) + if f: + self.__left.append((i, f)) + f = face(map(lambda e: e[2], element1)) + if f: + self.__top.append((i, f)) + + for i in range(5): + element0 = mesh.elements().tetUnitCube(i, 0) + + f = face(map(lambda e: e[0], element0)) + if f: + self.__right.append((i, f)) + f = face(map(lambda e: e[1], element0)) + if f: + self.__left.append((i, f)) + f = face(map(lambda e: e[2], element0)) + if f: + self.__top.append((i, f)) + + f = face(map(lambda e: 1-e[0], element0)) + if f: + self.__leftb.append((i, f)) + f = face(map(lambda e: 1-e[1], element0)) + if f: + self.__rightb.append((i, f)) + f = face(map(lambda e: 1-e[2], element0)) + if f: + self.__bot.append((i, f)) + + for i in range(5): + element1 = mesh.elements().tetUnitCube(i, 1) + + f = face(map(lambda e: 1-e[0], element1)) + if f: + self.__leftb.append((i, f)) + f = face(map(lambda e: 1-e[1], element1)) + if f: + self.__rightb.append((i, f)) + f = face(map(lambda e: 1-e[2], element1)) + if f: + self.__bot.append((i, f)) + + assert(len(self.__top) == 4) + assert(len(self.__left) == 4) + assert(len(self.__right) == 4) + assert(len(self.__bot) == 4) + assert(len(self.__rightb) == 4) + assert(len(self.__leftb) == 4) + + def __len__(self): + return self.__len + + def __getitem__(self, index): + if index >= self.__len: + raise StopIteration + + len1All = self.__len1 + self.__len1 + for side in range(6): + if index < len1All[side]: + break + index -= len1All[side] + + def getFace(index, sizeX, sizeY): + x = (index / 2) % sizeX + y = (index / (2 * sizeX)) % sizeY + + face = index % 4 + + if y%2 == 1: + # odd rows ... + face = (face+2)%4 + + return (x, y, face) + + size = [2 * i for i in self.__mesh.size()] # number of cubes in each dimension + + if side == 0: # top (x = x; y = y) + (x, y, face) = getFace(index, size[0], size[1]) + + offset = (x + y*size[0]) * 5 + f = self.__top[face] + elif side == 1: # left (x = x; y = z) + (x, y, face) = getFace(index, size[0], size[2]) + + offset = (x + y*size[1]*size[0]) * 5 + f = self.__left[face] + elif side == 2: # right (x = y; y = z) + (x, y, face) = getFace(index, size[1], size[2]) + + offset = (x + y*size[1]) * size[0] * 5 + f = self.__right[face] + elif side == 3: # bottom (x = x; y = y) + (x, y, face) = getFace(index, size[0], size[1]) + + offset = (x + (y + (size[2]-1)*size[1])*size[0]) * 5 + f = self.__bot[face] + elif side == 4: # right behind (x = x; y = z) + (x, y, face) = getFace(index, size[0], size[2]) + + offset = (x + (y*size[1] + size[1]-1)*size[0]) * 5 + f = self.__rightb[face] + else: # side == 5 # left behind (x = y; y = z) + (x, y, face) = getFace(index, size[1], size[2]) + + offset = ((x + y*size[1]) * size[0] + size[0]-1) * 5 + f = self.__leftb[face] + + return (f[0]+offset, f[1]) + + def __init__(self, name, mesh): + self._name = name + self._sides = Boundary.Sides(mesh) + self._size = len(self._sides) + + def __getitem__(self, key): + return getattr(self, '_'+str(key)) + +class Mesh: + version = '2.0.0' + name = 'cube' + date = datetime.datetime.now().strftime('%b %Y') + + def __init__(self, size, boundaryCond): + # The smallest cube for this mesh has size 2*2*2 + self.__size = [i / 2 for i in size] + + self.name += '_'.join(map(str, size)) + + self.__coords = Coords(self) + self.__elements = Elements(self) + self.__groups = [Group('fluid', 2, self)] + self.__boundaries = [Boundary(boundaryCond, self)] + + def size(self): + return self.__size + + def coords(self): + return self.__coords + + def elements(self): + return self.__elements + + def groups(self): + return self.__groups + + def boundaries(self): + return self.__boundaries diff --git a/cube/lib/netcdf.py b/cube/lib/netcdf.py new file mode 100644 index 0000000..d153fd6 --- /dev/null +++ b/cube/lib/netcdf.py @@ -0,0 +1,509 @@ +#!/usr/bin/python +## +# @file +# This file is part of SeisSol. +# +# @author Sebastian Rettenberger (rettenbs AT in.tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger,_M.Sc.) +# +# @section LICENSE +# Copyright (c) 2013, 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 +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# TODO Currently all partitions must have the same size !!! +# + +import operator +import functools +from netCDF4 import Dataset + +def unique(seq, idfun=None): + # order preserving + if idfun is None: + def idfun(x): return x + seen = {} + result = [] + for item in seq: + marker = idfun(item) + # in old Python versions: + # if seen.has_key(marker) + # but in new ones: + if marker in seen: continue + seen[marker] = 1 + result.append(item) + return result + +class NetcdfWriter: + def __init__(self, mesh, partitions, file): + file.close() # Close the file object because argsparse opens it automatically + + numPartitions = reduce(operator.mul, partitions, 1) + partitionCubes = [(mesh.size()[i] + partitions[i] - 1) / partitions[i] * 2 for i in range(3)] + partitionSize = reduce(operator.mul, partitionCubes, 5) # 5*8 elements per double cube + for i in range(3): + if mesh.size()[i] % partitions[i] != 0: + raise Exception('Different partition sizes currently not supported') + + file = Dataset(file.name, 'w', format='NETCDF4') + + # Dimensions + dimDimension = file.createDimension('dimension', 3) + dimPartitions = file.createDimension('partitions', numPartitions) + + dimElem = file.createDimension('elements', partitionSize) + dimElemSides = file.createDimension('element_sides', 4) + dimElemVertices = file.createDimension('element_vertices_dim', 4) # Python bindings are to stupid so we need another name for the dimension + + vrtxSize = reduce(operator.mul, map(lambda i: i+1, partitionCubes), 1) + dimVrtx = file.createDimension('vertices', vrtxSize) + + boundarySizes = [2*partitionCubes[1]*partitionCubes[2], 2*partitionCubes[0]*partitionCubes[2], 2*partitionCubes[0]*partitionCubes[1]] + boundarySizes = boundarySizes[::-1] + boundarySizes + dimBnd = file.createDimension('boundaries', 6) + dimBndElem = file.createDimension('boundary_elements', max(boundarySizes)) + + # Variables + varElemSize = file.createVariable('element_size', 'i4', ('partitions',)) + varElemVertices = file.createVariable('element_vertices', 'i4', ('partitions', 'elements', 'element_vertices_dim')) + varElemNeighbors = file.createVariable('element_neighbors', 'i4', ('partitions', 'elements', 'element_sides')) + varElemBoundaries = file.createVariable('element_boundaries', 'i4', ('partitions', 'elements', 'element_sides')) + varElemNeighborSides = file.createVariable('element_neighbor_sides', 'i4', ('partitions', 'elements', 'element_sides')) + varElemSideOrientations = file.createVariable('element_side_orientations', 'i4', ('partitions', 'elements', 'element_sides')) + varElemNeighborRanks = file.createVariable('element_neighbor_ranks', 'i4', ('partitions', 'elements', 'element_sides')) + varElemMPIIndices = file.createVariable('element_mpi_indices', 'i4', ('partitions', 'elements', 'element_sides')) + + varVrtxSize = file.createVariable('vertex_size', 'i4', ('partitions',)) + varVrtxCoords = file.createVariable('vertex_coordinates', 'f8', ('partitions', 'vertices', 'dimension')) + + varBndSize = file.createVariable('boundary_size', 'i4', ('partitions',)) + varBndElemSize = file.createVariable('boundary_element_size', 'i4', ('partitions', 'boundaries')) + varBndElemRank = file.createVariable('boundary_element_rank', 'i4', ('partitions', 'boundaries')) + varBndElemLocalIds = file.createVariable('boundary_element_localids', 'i4', ('partitions', 'boundaries', 'boundary_elements')) + + # Elements + varElemSize[:] = numPartitions*[partitionSize] + +# for z in range(partitions[2]): +# for y in range(partitions[1]): +# for x in range(partitions[0]): + vertices = [] + + for zz in range(partitionCubes[2]): + for yy in range(partitionCubes[1]): + for xx in range(partitionCubes[0]): + for i in range(5): + #vertices.append(mesh.elements()[(((z*partitionCubes[2]+zz)*mesh.size()[1]*2 + (y*partitionCubes[1]+yy))*mesh.size()[0]*2 + (x*partitionCubes[0]+xx))*5 + i]) + vertices.append(mesh.elements()[((zz*mesh.size()[1]*2 + yy)*mesh.size()[0]*2 + xx)*5 + i]) + + verticesMap = {} + for i in range(len(vertices)): + for j in range(len(vertices[i])): + if vertices[i][j] in verticesMap: + vertices[i][j] = verticesMap[vertices[i][j]] + else: + n = len(verticesMap) + verticesMap[vertices[i][j]] = n + vertices[i][j] = n + + for i in range(numPartitions): + varElemVertices[i,:,:] = vertices + #varElemVertices[(z*partitions[1] + y)*partitions[0] + x,:,:] = vertices + + print 'element_vertices done' + + neighbors = [] + for z in range(partitionCubes[2]): + for y in range(partitionCubes[1]): + for x in range(partitionCubes[0]): + odd = (x + y + z) % 2 + + if odd: + n = [[partitionSize if x == 0 else -4, + partitionSize if z == 0 else -partitionCubes[1]*partitionCubes[0]*5+3, + 4, + partitionSize if y == partitionCubes[1]-1 else partitionCubes[0]*5], + [4, + partitionSize if z == 0 else -partitionCubes[1]*partitionCubes[0]*5+2, + partitionSize if y == 0 else -partitionCubes[0]*5+1, + partitionSize if x == partitionCubes[0]-1 else 5], + [4, + partitionSize if y == 0 else -partitionCubes[0]*5+3, + partitionSize if x == 0 else -3, + partitionSize if z == partitionCubes[2]-1 else partitionCubes[1]*partitionCubes[0]*5], + [partitionSize if x == partitionCubes[0]-1 else 8, + 4, + partitionSize if y == partitionCubes[1]-1 else partitionCubes[0]*5+2, + partitionSize if z == partitionCubes[2]-1 else partitionCubes[1]*partitionCubes[0]*5+1], + [0, 1, 2, 3]] + else: + n = [[partitionSize if z == 0 else -partitionCubes[1]*partitionCubes[0]*5+2, + partitionSize if y == 0 else -partitionCubes[0]*5, + partitionSize if x == 0 else -4, + 4], + [4, + partitionSize if z == 0 else -partitionCubes[1]*partitionCubes[0]*5+3, + partitionSize if x == partitionCubes[0]-1 else 5, + partitionSize if y == partitionCubes[1]-1 else partitionCubes[0]*5+1], + [4, + partitionSize if x == partitionCubes[0]-1 else 7, + partitionSize if y == 0 else -partitionCubes[0]*5+3, + partitionSize if z == partitionCubes[2]-1 else partitionCubes[1]*partitionCubes[0]*5+1], + [partitionSize if x == 0 else -2, + partitionSize if y == partitionCubes[1]-1 else partitionCubes[0]*5+2, + 4, + partitionSize if z == partitionCubes[2]-1 else partitionCubes[1]*partitionCubes[0]*5], + [0, 1, 2, 3]] + + offset = ((z*partitionCubes[1] + y)*partitionCubes[0] + x)*5 + n = map(lambda a: map(lambda i: i + offset if i < partitionSize else i, a), n) + neighbors.extend(n) + + for i in range(numPartitions): + varElemNeighbors[i,:,:] = neighbors + + print 'element_neighbors done' + + for z in range(partitions[2]): + for y in range(partitions[1]): + for x in range(partitions[0]): + boundaries = [] + + for zz in range(partitionCubes[2]): + for yy in range(partitionCubes[1]): + for xx in range(partitionCubes[0]): + odd = (xx + yy + zz) % 2 + + # TODO support different boundaries + + if odd: + b = [[1 if x == 0 and xx == 0 else 0, + 1 if z == 0 and zz == 0 else 0, + 0, + 1 if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 0], + [0, + 1 if z == 0 and zz == 0 else 0, + 1 if y == 0 and yy == 0 else 0, + 1 if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 0], + [0, + 1 if y == 0 and yy == 0 else 0, + 1 if x == 0 and xx == 0 else 0, + 1 if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 0], + [1 if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 0, + 0, + 1 if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 0, + 1 if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 0], + [0, 0, 0, 0]] + else: + b = [[1 if z == 0 and zz == 0 else 0, + 1 if y == 0 and yy == 0 else 0, + 1 if x == 0 and xx == 0 else 0, + 0], + [0, + 1 if z == 0 and zz == 0 else 0, + 1 if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 0, + 1 if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 0], + [0, + 1 if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 0, + 1 if y == 0 and yy == 0 else 0, + 1 if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 0], + [1 if x == 0 and xx == 0 else 0, + 1 if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 0, + 0, + 1 if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 0], + [0, 0, 0, 0]] + + boundaries.extend(b) + + varElemBoundaries[(z*partitions[1] + y)*partitions[0] + x,:,:] = boundaries + + print 'element_boundaries done' + + for z in range(partitions[2]): + for y in range(partitions[1]): + for x in range(partitions[0]): + sides = [] + + for zz in range(partitionCubes[2]): + for yy in range(partitionCubes[1]): + for xx in range(partitionCubes[0]): + odd = (xx + yy + zz) % 2 + + if odd: + s = [[0 if x == 0 and xx == 0 else 2, + 0 if z == 0 and zz == 0 else 3, + 0, + 0 if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 1], + [1, + 0 if z == 0 and zz == 0 else 3, + 0 if y == 0 and yy == 0 else 3, + 0 if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 2], + [2, + 0 if y == 0 and yy == 0 else 1, + 0 if x == 0 and xx == 0 else 1, + 0],# if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 0], + [0,# if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 0, + 3, + 0 if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 2, + 0 if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 1], + [2, 0, 0, 1]] + else: + s = [[0 if z == 0 and zz == 0 else 3, + 0 if y == 0 and yy == 0 else 3, + 0 if x == 0 and xx == 0 else 3, + 0], + [1, + 0 if z == 0 and zz == 0 else 3, + 0,# if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 0, + 0 if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 2], + [2, + 0 if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 2, + 0 if y == 0 and yy == 0 else 2, + 0 if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 1], + [0 if x == 0 and xx == 0 else 0, + 0 if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 1, + 3, + 0 if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 1], + [3, 0, 0, 2]] + + sides.extend(s) + + varElemNeighborSides[(z*partitions[1] + y)*partitions[0] + x,:,:] = sides + + print 'element_neighbor_sides done' + + for z in range(partitions[2]): + for y in range(partitions[1]): + for x in range(partitions[0]): + orientations = [] + + for zz in range(partitionCubes[2]): + for yy in range(partitionCubes[1]): + for xx in range(partitionCubes[0]): + odd = (xx + yy + zz) % 2 + + if odd: + o = [[0,# if x == 0 and xx == 0 else 0, + 0 if z == 0 and zz == 0 else 1, + 0, + 0],# if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 0], + [0, + 0 if z == 0 and zz == 0 else 1, + 0,# if y == 0 and yy == 0 else 0, + 0 if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 2], + [0, + 0,# if y == 0 and yy == 0 else 0, + 0,# if x == 0 and xx == 0 else 0, + 0 if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 2], + [0,# if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 0, + 0, + 0,# if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 0, + 0],# if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 0], + [0, 0, 0, 0]] + else: + o = [[0 if z == 0 and zz == 0 else 2, + 0,# if y == 0 and yy == 0 else 0, + 0 if x == 0 and xx == 0 else 2, + 0], + [0, + 0,# if z == 0 and zz == 0 else 0, + 0,# if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 0, + 0],# if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 0], + [0, + 0,# if x == partitions[0]-1 and xx == partitionCubes[0]-1 else 0, + 0,# if y == 0 and yy == 0 else 0, + 0 if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 1], + [0,# if x == 0 and xx == 0 else 0, + 0,# if y == partitions[1]-1 and yy == partitionCubes[1]-1 else 0, + 0, + 0 if z == partitions[2]-1 and zz == partitionCubes[2]-1 else 1], + [0, 0, 0, 0]] + + orientations.extend(o) + + varElemSideOrientations[(z*partitions[1] + y)*partitions[0] + x,:,:] = orientations + + print 'element_side_orientations done' + + for z in range(partitions[2]): + for y in range(partitions[1]): + for x in range(partitions[0]): + myrank = (z*partitions[1] + y)*partitions[0] + x + ranks = [] + + for zz in range(partitionCubes[2]): + for yy in range(partitionCubes[1]): + for xx in range(partitionCubes[0]): + odd = (xx + yy + zz) % 2 + + if odd: + r = [[myrank if xx != 0 else myrank if x == 0 else myrank-1, + myrank if zz != 0 else myrank if z == 0 else myrank-partitions[1]*partitions[0], + myrank, + myrank if yy != partitionCubes[1]-1 else myrank if y == partitions[1]-1 else myrank+partitions[0]], + [myrank, + myrank if zz != 0 else myrank if z == 0 else myrank-partitions[1]*partitions[0], + myrank if yy != 0 else myrank if y == 0 else myrank-partitions[0], + myrank if xx != partitionCubes[0]-1 else myrank if x == partitions[0]-1 else myrank+1], + [myrank, + myrank if yy != 0 else myrank if y == 0 else myrank-partitions[0], + myrank if xx != 0 else myrank if x == 0 else myrank-1, + myrank if zz != partitionCubes[2]-1 else myrank if z == partitions[2]-1 else myrank+partitions[1]*partitions[0]], + [myrank if xx != partitionCubes[0]-1 else myrank if x == partitions[0]-1 else myrank+1, + myrank, + myrank if yy != partitionCubes[1]-1 else myrank if y == partitions[1]-1 else myrank+partitions[0], + myrank if zz != partitionCubes[2]-1 else myrank if z == partitions[2]-1 else myrank+partitions[1]*partitions[0]], + [myrank, myrank, myrank, myrank]] + else: + r = [[myrank if zz != 0 else myrank if z == 0 else myrank-partitions[1]*partitions[0], + myrank if yy != 0 else myrank if y == 0 else myrank-partitions[0], + myrank if xx != 0 else myrank if x == 0 else myrank-1, + myrank], + [myrank, + myrank if zz != 0 else myrank if z == 0 else myrank-partitions[1]*partitions[0], + myrank if xx != partitionCubes[0]-1 else myrank if x == partitions[0]-1 else myrank+1, + myrank if yy != partitionCubes[1]-1 else myrank if y == partitions[1]-1 else myrank+partitions[0]], + [myrank, + myrank if xx != partitionCubes[0]-1 else myrank if x == partitions[0]-1 else myrank+1, + myrank if yy != 0 else myrank if y == 0 else myrank-partitions[0], + myrank if zz != partitionCubes[2]-1 else myrank if z == partitions[2]-1 else myrank+partitions[1]*partitions[0]], + [myrank if xx != 0 else myrank if x == 0 else myrank-1, + myrank if yy != partitionCubes[1]-1 else myrank if y == partitions[1]-1 else myrank+partitions[0], + myrank, + myrank if zz != partitionCubes[2]-1 else myrank if z == partitions[2]-1 else myrank+partitions[1]*partitions[0]], + [myrank, myrank, myrank, myrank]] + + ranks.extend(r) + + lastIndices = [-1]*6 + bndLocalIds = [None]*6 + for i in range(6): + bndLocalIds[i] = [] + def mpiIndex(e, r): + if r == myrank: + return 0 + if r == myrank-partitions[0]*partitions[1]: + lastIndices[0] += 1 + if ((lastIndices[0] % (partitionCubes[0]*4)) / (partitionCubes[0]*2) == 0 and lastIndices[0] % 4 == 2) \ + or ((lastIndices[0] % (partitionCubes[0]*4)) / (partitionCubes[0]*2) == 1 and (lastIndices[0] % 4) == 0): + bndLocalIds[0].extend([None, e]) + return lastIndices[0]+1 + elif ((lastIndices[0] % (partitionCubes[0]*4)) / (partitionCubes[0]*2) == 0 and lastIndices[0] % 4 == 3) \ + or ((lastIndices[0] % (partitionCubes[0]*4)) / (partitionCubes[0]*2) == 1 and (lastIndices[0] % 4) == 1): + bndLocalIds[0][len(bndLocalIds[0])-2] = e + return lastIndices[0]-1 + bndLocalIds[0].append(e) + return lastIndices[0] + if r == myrank-partitions[0]: + lastIndices[1] += 1 + bndLocalIds[1].append(e) + return lastIndices[1] + if r == myrank-1: + lastIndices[2] += 1 + bndLocalIds[2].append(e) + return lastIndices[2] + if r == myrank+1: + lastIndices[3] += 1 + bndLocalIds[3].append(e) + return lastIndices[3] + if r == myrank+partitions[0]: + lastIndices[4] += 1 + bndLocalIds[4].append(e) + return lastIndices[4] + if r == myrank+partitions[0]*partitions[1]: + lastIndices[5] += 1 + bndLocalIds[5].append(e) + return lastIndices[5] + raise Exception('Invalid neighbor rank found') + indices = map(lambda (e, r): map(functools.partial(mpiIndex, e), r), enumerate(ranks)) + varElemNeighborRanks[(z*partitions[1] + y)*partitions[0] + x,:,:] = ranks + varElemMPIIndices[(z*partitions[1] + y)*partitions[0] + x,:,:] = indices + + k = 0 + for i in range(6): + if bndLocalIds[i]: + varBndElemLocalIds[(z*partitions[1] + y)*partitions[0] + x,k,0:len(bndLocalIds[i])] = bndLocalIds[i] + k += 1 + + print 'element_neighbor_ranks, element_mpi_indices and boundary_element_localids done' + + varVrtxSize[:] = vrtxSize + + for z in range(partitions[2]): + for y in range(partitions[1]): + for x in range(partitions[0]): + globalVrtxIds = [] + + for zz in range(partitionCubes[2]): + for yy in range(partitionCubes[1]): + for xx in range(partitionCubes[0]): + for i in range(5): + globalVrtxIds += mesh.elements()[(((z*partitionCubes[2]+zz)*mesh.size()[1]*2 + (y*partitionCubes[1]+yy))*mesh.size()[0]*2 + (x*partitionCubes[0]+xx))*5 + i] + + coords = map(lambda v: mesh.coords()[v], unique(globalVrtxIds)) + varVrtxCoords[(z*partitions[1] + y)*partitions[0] + x,0:vrtxSize,:] = coords + + print 'vertex_coordinates done' + + boundaries = [None]*numPartitions + for z in range(partitions[2]): + for y in range(partitions[1]): + for x in range(partitions[0]): + b = [True]*6 + if x == 0: + b[2] = False + if y == 0: + b[1] = False + if z == 0: + b[0] = False + if x == partitions[0]-1: + b[3] = False + if y == partitions[1]-1: + b[4] = False + if z == partitions[2]-1: + b[5] = False + boundaries[(z*partitions[1] + y)*partitions[0] + x] = b + + varBndSize[:] = map(sum, boundaries) + + boundaryRankDiff = [1, partitions[0], partitions[0]*partitions[1]] + boundaryRankDiff = map(lambda i: -i, boundaryRankDiff[::-1]) + boundaryRankDiff + + for i in range(numPartitions): + bndElemSize = [] + bndElemRank = [] + + for j in range(6): + if boundaries[i][j]: + bndElemSize.append(boundarySizes[j]) + bndElemRank.append(boundaryRankDiff[j]+i) + + varBndElemSize[i,0:len(bndElemSize)] = bndElemSize + varBndElemRank[i,0:len(bndElemRank)] = bndElemRank + + file.close() diff --git a/cube_c/CMakeLists.txt b/cube_c/CMakeLists.txt new file mode 100644 index 0000000..3d53f8a --- /dev/null +++ b/cube_c/CMakeLists.txt @@ -0,0 +1,58 @@ +cmake_minimum_required(VERSION 3.10) + +# set the project name +project(cube_c) + +# specify the C++ standard +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED True) + +cmake_policy(SET CMP0074 NEW) + + +add_executable(cubeGenerator src/main.cpp) + +#logging +set(LOG_LEVEL "warning" CACHE STRING "Log level for the code") +set(LOG_LEVEL_OPTIONS "debug" "info" "warning" "error") +set_property(CACHE LOG_LEVEL PROPERTY STRINGS ${LOG_LEVEL_OPTIONS}) +if("${LOG_LEVEL}" STREQUAL "debug") + target_compile_definitions(cubeGenerator PUBLIC LOG_LEVEL=3) +elseif("${LOG_LEVEL}" STREQUAL "info") + target_compile_definitions(cubeGenerator PUBLIC LOG_LEVEL=2) +elseif("${LOG_LEVEL}" STREQUAL "warning") + target_compile_definitions(cubeGenerator PUBLIC LOG_LEVEL=1) +elseif("${LOG_LEVEL}" STREQUAL "error") + target_compile_definitions(cubeGenerator PUBLIC LOG_LEVEL=0) +endif() + +#build and link libraries and executable +set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH}) +find_package(NetCDF REQUIRED) +target_include_directories(cubeGenerator PUBLIC ${NetCDF_INCLUDE_DIRS}) +target_link_libraries(cubeGenerator PUBLIC ${NetCDF_LIBRARY}) + +find_package(HDF5 REQUIRED + COMPONENTS C HL) +target_include_directories(cubeGenerator PUBLIC ${HDF5_INCLUDE_DIRS}) +target_link_libraries(cubeGenerator PUBLIC ${HDF5_C_HL_LIBRARIES} ${HDF5_C_LIBRARIES}) + +find_package(MPI REQUIRED) +target_link_libraries(cubeGenerator PUBLIC MPI::MPI_C) + +find_package(OpenMP REQUIRED) +target_link_libraries(cubeGenerator PUBLIC OpenMP::OpenMP_CXX) + +#add some compiler specific flags +if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") + # using GCC + target_compile_options(cubeGenerator PUBLIC -fopenmp -pedantic $<$,$>:-Wall -Wextra -Wno-unused-parameter -Wno-unknown-pragmas>) + target_link_libraries(cubeGenerator PUBLIC "-fopenmp") +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") + target_compile_options(cubeGenerator PUBLIC -qopenmp -pedantic $<$,$>:-Wall -w3 -diag-disable:remark>) + target_link_libraries(cubeGenerator PUBLIC "-qopenmp") +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp=libomp -Wall -Wextra -pedantic") +endif() + +target_include_directories(cubeGenerator PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/submodules) diff --git a/cube_c/cmake b/cube_c/cmake new file mode 120000 index 0000000..eca6da1 --- /dev/null +++ b/cube_c/cmake @@ -0,0 +1 @@ +/home/david/Dokumente/uni/phd/software/meshing/Meshing/cmake \ No newline at end of file diff --git a/cube_c/src/main.cpp b/cube_c/src/main.cpp new file mode 100644 index 0000000..9c5ba98 --- /dev/null +++ b/cube_c/src/main.cpp @@ -0,0 +1,1290 @@ +/** + * @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) + * + * @section LICENSE + * Copyright (c) 2015, 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 + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * + * @section DESCRIPTION + **/ + +#include "utils/args.h" +#include "utils/logger.h" + +#include + +#include + +#include +#include +#include +#include +#include + +typedef int t_vertex[3]; + +// Index of the vertices of a tetraedra in a cube +// even/odd, index of the tetrahedra, index of vertex, offset of the vertices in x/y/z +static const t_vertex TET_VERTICES[2][5][4] = { + { + {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, + {{1, 0, 0}, {0, 1, 0}, {1, 1, 1}, {1, 1, 0}}, + {{1, 0, 0}, {1, 1, 1}, {0, 0, 1}, {1, 0, 1}}, + {{0, 1, 0}, {0, 1, 1}, {0, 0, 1}, {1, 1, 1}}, + {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {1, 1, 1}} + }, + { + {{0, 0, 0}, {0, 1, 0}, {0, 1, 1}, {1, 1, 0}}, + {{0, 0, 0}, {1, 1, 0}, {1, 0, 1}, {1, 0, 0}}, + {{0, 0, 0}, {1, 0, 1}, {0, 1, 1}, {0, 0, 1}}, + {{1, 1, 0}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}}, + {{0, 0, 0}, {1, 1, 0}, {0, 1, 1}, {1, 0, 1}} + } +}; + +struct Vertex +{ + t_vertex v; + + bool operator<(const Vertex& other) const { + return (v[0] < other.v[0]) + || ((v[0] == other.v[0]) && (v[1] < other.v[1])) + || ((v[0] == other.v[0]) && (v[1] == other.v[1]) && (v[2] < other.v[2])); + } +}; + +static const int TET_SIDE_NEIGHBORS[2][5*4] = { + {3, 3, 3, 0, + 1, 3, 0, 2, + 2, 2, 2, 1, + 0, 1, 3, 1, + 3, 0, 0, 2}, + {2, 3, 0, 1, + 1, 3, 3, 2, + 2, 1, 1, 0, + 0, 3, 2, 1, + 2, 0, 0, 1} +}; + +static const int TET_SIDE_ORIENTATIONS[2][5*4] = { + {2, 0, 2, 0, + 0, 0, 0, 0, + 0, 0, 0, 1, + 0, 0, 0, 1, + 0, 0, 0, 0}, + {0, 1, 0, 0, + 0, 1, 0, 2, + 0, 0, 0, 2, + 0, 0, 0, 0, + 0, 0, 0, 0} +}; + +// Process has done i out of n rounds, +// and we want a bar of width w and resolution r. +// TODO put this i an library +static inline void loadBar(int x, int n, int r = 100, int w = 50) +{ + r = std::min(r, n); + + // Only update r times. + if ( x % (n/r) != 0 ) return; + + // Calculuate the ratio of complete-to-incomplete. + float ratio = x/(float)n; + int c = ratio * w; + + // Show the percentage complete. + std::cout << std::setw(3) << (int)(ratio*100) << "% ["; + + // Show the load bar. + for (int x=0; x +std::pair flip_pair(const std::pair &p) +{ + return std::pair(p.second, p.first); +} + +void checkNcError(int error) +{ + if (error != NC_NOERR) + logError() << "Error while writing netCDF file:" << nc_strerror(error); +} + +int main(int argc, char* argv[]) +{ + // Parse command line arguments + utils::Args args; + args.addOption("boundary", 'b', "boundary condition (default: 1)", utils::Args::Required, false); + args.addOption("minx", 0, "boundary condition at the start of x dimension", utils::Args::Required, false); + args.addOption("maxx", 0, "boundary condition at the end of x dimension", utils::Args::Required, false); + args.addOption("miny", 0, "boundary condition at the start of y dimension", utils::Args::Required, false); + args.addOption("maxy", 0, "boundary condition at the end of y dimension", utils::Args::Required, false); + args.addOption("minz", 0, "boundary condition at the start of z dimension", utils::Args::Required, false); + args.addOption("maxz", 0, "boundary condition at the end of z dimension", utils::Args::Required, false); + args.addOption("netcdf", 'n', "for compatibility with the Python script", utils::Args::No, false); + args.addOption("x", 'x', "number of cubes in x dimension"); + args.addOption("y", 'y', "number of cubes in y dimension"); + args.addOption("z", 'z', "number of cubes in z dimension"); + args.addOption("px", 0, "number of partitions x dimension"); + args.addOption("py", 0, "number of partitions y dimension"); + args.addOption("pz", 0, "number of partitions z dimension"); + args.addOption("output", 'o', "output file for resulting netCDF mesh"); + args.addOption("scale", 's', "size of the domain = [-s/2, s/2]^3 (default: 100)", utils::Args::Required, false); + args.addOption("sx", 0, "scale in x direction (overwrites scale in x direction)", utils::Args::Required, false); + args.addOption("sy", 0, "scale in y direction (overwrites scale in y direction)", utils::Args::Required, false); + args.addOption("sz", 0, "scale in z direction (overwrites scale in z direction)", utils::Args::Required, false); + args.addOption("tx", 0, "mesh translation in x direction (after scaling)", utils::Args::Required, false); + args.addOption("ty", 0, "mesh translation in y direction (after scaling)", utils::Args::Required, false); + args.addOption("tz", 0, "mesh translation in z direction (after scaling)", utils::Args::Required, false); + + unsigned int numCubes[4]; + unsigned int numPartitions[4]; + + // Parse/check command line arguments + if (args.parse(argc, argv) == utils::Args::Success) { + numCubes[0] = args.getArgument("x"); + numCubes[1] = args.getArgument("y"); + numCubes[2] = args.getArgument("z"); + + numPartitions[0] = args.getArgument("px"); + numPartitions[1] = args.getArgument("py"); + numPartitions[2] = args.getArgument("pz"); + + for (int i = 0; i < 3; i++) { + if (numCubes[i] < 2) + logError() << "Number of cubes in" << dim2str(i) << "dimension must be at least 2"; + if (numCubes[i] % numPartitions[i] != 0) + logError() << "Number of cubes in" << dim2str(i) + << "dimension can not be distribute to" << numPartitions[i] << "partitions"; + if ((numCubes[i] / numPartitions[i]) % 2 != 0) + logError() << "Number of cubes per partition in" << dim2str(i) << "dimension must be a multiple of 2"; + } + } else { + return 1; + } + + unsigned int boundary = args.getArgument("boundary", 1u); + if (boundary > 100) + boundary -= 100; + unsigned int boundaryMinx = args.getArgument("minx", boundary); + unsigned int boundaryMaxx = args.getArgument("maxx", boundary); + unsigned int boundaryMiny = args.getArgument("miny", boundary); + unsigned int boundaryMaxy = args.getArgument("maxy", boundary); + unsigned int boundaryMinz = args.getArgument("minz", boundary); + unsigned int boundaryMaxz = args.getArgument("maxz", boundary); + logInfo() << "Boundary condition:" << boundary; + logInfo() << "boundaryMinx condition:" << boundaryMinx; + logInfo() << "boundaryMaxx condition:" << boundaryMaxx; + logInfo() << "boundaryMiny condition:" << boundaryMiny; + logInfo() << "boundaryMaxy condition:" << boundaryMaxy; + logInfo() << "boundaryMinz condition:" << boundaryMinz; + logInfo() << "boundaryMaxz condition:" << boundaryMaxz; + + // Compute additional sizes + numCubes[3] = numCubes[0] * numCubes[1] * numCubes[2]; + numPartitions[3] = numPartitions[0] * numPartitions[1] * numPartitions[2]; + + unsigned int numCubesPerPart[4]; + unsigned long numElemPerPart[4]; + for (int i = 0; i < 4; i++) { + numCubesPerPart[i] = numCubes[i] / numPartitions[i]; + numElemPerPart[i] = numCubesPerPart[i] * 5; + } + + unsigned int numVrtxPerPart[4]; + for (int i = 0; i < 3; i++) + numVrtxPerPart[i] = numCubesPerPart[i] + 1; + numVrtxPerPart[3] = numVrtxPerPart[0] * numVrtxPerPart[1] * numVrtxPerPart[2]; + + unsigned int numBndElements[3]; + numBndElements[0] = 2*numCubesPerPart[1]*numCubesPerPart[2]; + numBndElements[1] = 2*numCubesPerPart[0]*numCubesPerPart[2]; + numBndElements[2] = 2*numCubesPerPart[0]*numCubesPerPart[1]; + + logInfo() << "Total number of cubes:" << numCubes[0] << 'x' << numCubes[1] << 'x' << numCubes[2] << '=' << numCubes[3]; + logInfo() << "Total number of partitions" << numPartitions[0] << 'x' << numPartitions[1] << 'x' << numPartitions[2] + << '=' << numPartitions[3]; + logInfo() << "Total number of cubes per partition:" << numCubesPerPart[0] << 'x' << numCubesPerPart[1] << 'x' + << numCubesPerPart[2] << '=' << numCubesPerPart[3]; + logInfo() << "Total number of elements per partition:" << numElemPerPart[0] << 'x' << numElemPerPart[1] << 'x' + << numElemPerPart[2] << '=' << numElemPerPart[3]; + logInfo() << "Using" << omp_get_max_threads() << "threads"; + + int netcdf_writes = 2 + numPartitions[3]*8; + + // Create the netcdf file + int ncFile; + checkNcError(nc_create(args.getArgument("output").c_str(), NC_NETCDF4, &ncFile)); + + int ncDimDimension; + nc_def_dim(ncFile, "dimension", 3, &ncDimDimension); + + int ncDimPart; + nc_def_dim(ncFile, "partitions", numPartitions[3], &ncDimPart); + + int ncDimElem, ncDimElemSides, ncDimElemVertices; + nc_def_dim(ncFile, "elements", numElemPerPart[3], &ncDimElem); + nc_def_dim(ncFile, "element_sides", 4, &ncDimElemSides); + nc_def_dim(ncFile, "element_vertices_dim", 4, &ncDimElemVertices); + + int ncDimVrtx; + nc_def_dim(ncFile, "vertices", numVrtxPerPart[3], &ncDimVrtx); + + int ncDimBnd, ncDimBndElem; + nc_def_dim(ncFile, "boundaries", 6, &ncDimBnd); + nc_def_dim(ncFile, "boundary_elements", *std::max_element(numBndElements, numBndElements+3), &ncDimBndElem); + + // Create netcdf variables + int ncVarElemSize; + checkNcError(nc_def_var(ncFile, "element_size", NC_INT, 1, &ncDimPart, &ncVarElemSize)); +// checkNcError(nc_var_par_access(ncFile, ncVarElemSize, NC_COLLECTIVE)); + + int ncVarElemVertices; + int dimsElemVertices[] = {ncDimPart, ncDimElem, ncDimElemVertices}; + checkNcError(nc_def_var(ncFile, "element_vertices", NC_INT, 3, dimsElemVertices, &ncVarElemVertices)); +// checkNcError(nc_var_par_access(ncFile, ncVarElemVertices, NC_COLLECTIVE)); + + int ncVarElemNeighbors; + int dimsElemSides[] = {ncDimPart, ncDimElem, ncDimElemSides}; + checkNcError(nc_def_var(ncFile, "element_neighbors", NC_INT, 3, dimsElemSides, &ncVarElemNeighbors)); +// checkNcError(nc_var_par_access(ncFile, ncVarElemNeighbors, NC_COLLECTIVE)); + + int ncVarElemBoundaries; + checkNcError(nc_def_var(ncFile, "element_boundaries", NC_INT, 3, dimsElemSides, &ncVarElemBoundaries)); +// checkNcError(nc_var_par_access(ncFile, ncVarElemBoundaries, NC_COLLECTIVE)); + + int ncVarElemNeighborSides; + checkNcError(nc_def_var(ncFile, "element_neighbor_sides", NC_INT, 3, dimsElemSides, &ncVarElemNeighborSides)); +// checkNcError(nc_var_par_access(ncFile, ncVarElemNeighborSides, NC_COLLECTIVE)); + + int ncVarElemSideOrientations; + checkNcError(nc_def_var(ncFile, "element_side_orientations", NC_INT, 3, dimsElemSides, &ncVarElemSideOrientations)); +// checkNcError(nc_var_par_access(ncFile, ncVarElemSideOrientations, NC_COLLECTIVE)); + + int ncVarElemNeighborRanks; + checkNcError(nc_def_var(ncFile, "element_neighbor_ranks", NC_INT, 3, dimsElemSides, &ncVarElemNeighborRanks)); +// checkNcError(nc_var_par_access(ncFile, ncVarElemNeighborRanks, NC_COLLECTIVE)); + + int ncVarElemMPIIndices; + checkNcError(nc_def_var(ncFile, "element_mpi_indices", NC_INT, 3, dimsElemSides, &ncVarElemMPIIndices)); +// checkNcError(nc_var_par_access(ncFile, ncVarElemMPIIndices, NC_COLLECTIVE)); + + int ncVarElemGroup; + checkNcError(nc_def_var(ncFile, "element_group", NC_INT, 2, dimsElemSides, &ncVarElemGroup)); +// checkNcError(nc_var_par_access(ncFile, ncVarElemGroup, NC_COLLECTIVE)); + + int ncVarVrtxSize; + checkNcError(nc_def_var(ncFile, "vertex_size", NC_INT, 1, &ncDimPart, &ncVarVrtxSize)); +// checkNcError(nc_var_par_access(ncFile, ncVarVrtxSize, NC_COLLECTIVE)); + + int ncVarVrtxCoords; + int dimsVrtxCoords[] = {ncDimPart, ncDimVrtx, ncDimDimension}; + checkNcError(nc_def_var(ncFile, "vertex_coordinates", NC_DOUBLE, 3, dimsVrtxCoords, &ncVarVrtxCoords)); +// checkNcError(nc_var_par_access(ncFile, ncVarVrtxCoords, NC_COLLECTIVE)); + + int ncVarBndSize; + checkNcError(nc_def_var(ncFile, "boundary_size", NC_INT, 1, &ncDimPart, &ncVarBndSize)); +// checkNcError(nc_var_par_access(ncFile, ncVarBndSize, NC_COLLECTIVE)); + + int ncVarBndElemSize; + int dimsBndElemSize[] = {ncDimPart, ncDimBnd}; + checkNcError(nc_def_var(ncFile, "boundary_element_size", NC_INT, 2, dimsBndElemSize, &ncVarBndElemSize)); +// checkNcError(nc_var_par_access(ncFile, ncVarBndElemSize, NC_COLLECTIVE)); + + int ncVarBndElemRank; + checkNcError(nc_def_var(ncFile, "boundary_element_rank", NC_INT, 2, dimsBndElemSize, &ncVarBndElemRank)); +// checkNcError(nc_var_par_access(ncFile, ncVarBndElemRank, NC_COLLECTIVE)); + + int ncVarBndElemLocalIds; + int dimsBndElemLocalIds[] = {ncDimPart, ncDimBnd, ncDimBndElem}; + checkNcError(nc_def_var(ncFile, "boundary_element_localids", NC_INT, 3, dimsBndElemLocalIds, &ncVarBndElemLocalIds)); +// checkNcError(nc_var_par_access(ncFile, ncVarBndElemLocalIds, NC_COLLECTIVE)); + + int writes_done = 0; + loadBar(writes_done, netcdf_writes); + + // Elements + int *elemSize = new int[numPartitions[3]]; + std::fill(elemSize, elemSize+numPartitions[3], numElemPerPart[3]); + checkNcError(nc_put_var_int(ncFile, ncVarElemSize, elemSize)); + delete [] elemSize; + writes_done++; loadBar(writes_done, netcdf_writes); + + std::vector vertices; + vertices.resize(numElemPerPart[3]*4); + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + unsigned int c = ((zz*numCubesPerPart[1] + yy)*numCubesPerPart[0] + xx) * 20; + int odd = (zz+yy+xx) % 2; + + for (unsigned int i = 0; i < 5; i++) { + for (unsigned int j = 0; j < 4; j++) { + Vertex v; + v.v[0] = TET_VERTICES[odd][i][j][0] + xx; + v.v[1] = TET_VERTICES[odd][i][j][1] + yy; + v.v[2] = TET_VERTICES[odd][i][j][2] + zz; + vertices[c] = v; + c++; + } + } + } + } + } + + int *elemVertices = new int[numElemPerPart[3]*4]; + std::map vertexMap; + for (unsigned int i = 0; i < vertices.size(); i++) { + std::map::iterator it = vertexMap.find(vertices[i]); + if (it != vertexMap.end()) { + elemVertices[i] = it->second; + } else { + int n = vertexMap.size(); + vertexMap[vertices[i]] = n; + elemVertices[i] = n; + } + } + + for (unsigned int i = 0; i < numPartitions[3]; i++) { + size_t start[3] = {i, 0, 0}; + size_t count[3] = {1, numElemPerPart[3], 4}; + checkNcError(nc_put_vara_int(ncFile, ncVarElemVertices, start, count, elemVertices)); + writes_done++; loadBar(writes_done, netcdf_writes); + } + delete [] elemVertices; + + int *elemNeighbors = new int[numElemPerPart[3]*4]; + const int TET_NEIGHBORS[2][5*4] = { + {-static_cast(numCubesPerPart[1]*numCubesPerPart[0])*5+2, -static_cast(numCubesPerPart[0])*5, -4, 4, + 4, -static_cast(numCubesPerPart[1]*numCubesPerPart[0])*5+3, 5, static_cast(numCubesPerPart[0])*5+1, + 4, 7, -static_cast(numCubesPerPart[0])*5+3, static_cast(numCubesPerPart[1]*numCubesPerPart[0])*5+1, + -2, static_cast(numCubesPerPart[0])*5+2, 4, static_cast(numCubesPerPart[1]*numCubesPerPart[0])*5, + 0, 1, 2, 3}, + {-4, -static_cast(numCubesPerPart[1]*numCubesPerPart[0])*5+3, 4, static_cast(numCubesPerPart[0])*5, + 4, -static_cast(numCubesPerPart[1]*numCubesPerPart[0])*5+2, -static_cast(numCubesPerPart[0])*5+1, 5, + 4, -static_cast(numCubesPerPart[0])*5+3, -3, static_cast(numCubesPerPart[1]*numCubesPerPart[0])*5, + 8, 4, static_cast(numCubesPerPart[0])*5+2, static_cast(numCubesPerPart[1]*numCubesPerPart[0])*5+1, + 0, 1, 2, 3} + }; + + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + unsigned int c = ((zz*numCubesPerPart[1] + yy)*numCubesPerPart[0] + xx) * 20; + int odd = (zz+yy+xx) % 2; + + memcpy(&elemNeighbors[c], TET_NEIGHBORS[odd], sizeof(int)*20); + int offset = ((zz*numCubesPerPart[1] + yy)*numCubesPerPart[0] + xx)*5; + for (int i = 0; i < 20; i++) + elemNeighbors[c+i] += offset; + + if (xx == 0) { // first cube in a partition in x dimension + if (odd) { + if (boundaryMinx == 6 && numPartitions[0] == 1) { + elemNeighbors[c] += numCubesPerPart[0]*5; + elemNeighbors[c+10] += numCubesPerPart[0]*5; + } else { + elemNeighbors[c] = numElemPerPart[3]; + elemNeighbors[c+10] = numElemPerPart[3]; + } + } else { + if (boundaryMinx == 6 && numPartitions[0] == 1) { + elemNeighbors[c+2] += numCubesPerPart[0]*5; + elemNeighbors[c+12] += numCubesPerPart[0]*5; + } else { + elemNeighbors[c+2] = numElemPerPart[3]; + elemNeighbors[c+12] = numElemPerPart[3]; + } + } + } else if (xx == numCubesPerPart[0]-1) { // last cube in a partition in x dimension + if (odd) { + if (boundaryMaxx == 6 && numPartitions[0] == 1) { + elemNeighbors[c+7] -= numCubesPerPart[0]*5; + elemNeighbors[c+12] -= numCubesPerPart[0]*5; + } else { + elemNeighbors[c+7] = numElemPerPart[3]; + elemNeighbors[c+12] = numElemPerPart[3]; + } + } else { + if (boundaryMaxx == 6 && numPartitions[0] == 1) { + elemNeighbors[c+6] -= numCubesPerPart[0]*5; + elemNeighbors[c+9] -= numCubesPerPart[0]*5; + } else { + elemNeighbors[c+6] = numElemPerPart[3]; + elemNeighbors[c+9] = numElemPerPart[3]; + } + } + } + if (yy == 0) { // first cube in a partition in y dimension + if (odd) { + if (boundaryMiny == 6 && numPartitions[1] == 1) { + elemNeighbors[c+6] += numCubesPerPart[0]*numCubesPerPart[1]*5; + elemNeighbors[c+9] += numCubesPerPart[0]*numCubesPerPart[1]*5; + } else { + elemNeighbors[c+6] = numElemPerPart[3]; + elemNeighbors[c+9] = numElemPerPart[3]; + } + } else { + if (boundaryMiny == 6 && numPartitions[1] == 1) { + elemNeighbors[c+1] += numCubesPerPart[0]*numCubesPerPart[1]*5; + elemNeighbors[c+10] += numCubesPerPart[0]*numCubesPerPart[1]*5; + } else { + elemNeighbors[c+1] = numElemPerPart[3]; + elemNeighbors[c+10] = numElemPerPart[3]; + } + } + } else if (yy == numCubesPerPart[1]-1) { // last cube in a partition in y dimension + if (odd) { + if (boundaryMaxy == 6 && numPartitions[1] == 1) { + elemNeighbors[c+3] -= numCubesPerPart[0]*numCubesPerPart[1]*5; + elemNeighbors[c+14] -= numCubesPerPart[0]*numCubesPerPart[1]*5; + } else { + elemNeighbors[c+3] = numElemPerPart[3]; + elemNeighbors[c+14] = numElemPerPart[3]; + } + } else { + if (boundaryMaxy == 6 && numPartitions[1] == 1) { + elemNeighbors[c+7] -= numCubesPerPart[0]*numCubesPerPart[1]*5; + elemNeighbors[c+13] -= numCubesPerPart[0]*numCubesPerPart[1]*5; + } else { + elemNeighbors[c+7] = numElemPerPart[3]; + elemNeighbors[c+13] = numElemPerPart[3]; + } + } + } + if (zz == 0) { // first cube in a partition in z dimension + if (odd) { + if (boundaryMinz == 6 && numPartitions[2] == 1) { + elemNeighbors[c+1] += numCubesPerPart[0]*numCubesPerPart[1]*numCubesPerPart[2]*5; + elemNeighbors[c+5] += numCubesPerPart[0]*numCubesPerPart[1]*numCubesPerPart[2]*5; + } else { + elemNeighbors[c+1] = numElemPerPart[3]; + elemNeighbors[c+5] = numElemPerPart[3]; + } + } else { + if (boundaryMinz == 6 && numPartitions[2] == 1) { + elemNeighbors[c] += numCubesPerPart[0]*numCubesPerPart[1]*numCubesPerPart[2]*5; + elemNeighbors[c+5] += numCubesPerPart[0]*numCubesPerPart[1]*numCubesPerPart[2]*5; + } else { + elemNeighbors[c] = numElemPerPart[3]; + elemNeighbors[c+5] = numElemPerPart[3]; + } + } + } else if (zz == numCubesPerPart[2]-1) { // last cube in a partition in z dimension + if (odd) { + if (boundaryMaxz == 6 && numPartitions[2] == 1) { + elemNeighbors[c+11] -= numCubesPerPart[0]*numCubesPerPart[1]*numCubesPerPart[2]*5; + elemNeighbors[c+15] -= numCubesPerPart[0]*numCubesPerPart[1]*numCubesPerPart[2]*5; + } else { + elemNeighbors[c+11] = numElemPerPart[3]; + elemNeighbors[c+15] = numElemPerPart[3]; + } + } else { + if (boundaryMaxz == 6 && numPartitions[2] == 1) { + elemNeighbors[c+11] -= numCubesPerPart[0]*numCubesPerPart[1]*numCubesPerPart[2]*5; + elemNeighbors[c+15] -= numCubesPerPart[0]*numCubesPerPart[1]*numCubesPerPart[2]*5; + } else { + elemNeighbors[c+11] = numElemPerPart[3]; + elemNeighbors[c+15] = numElemPerPart[3]; + } + } + } + } + } + } + + for (unsigned int i = 0; i < numPartitions[3]; i++) { + size_t start[3] = {i, 0, 0}; + size_t count[3] = {1, numElemPerPart[3], 4}; + checkNcError(nc_put_vara_int(ncFile, ncVarElemNeighbors, start, count, elemNeighbors)); + writes_done++; loadBar(writes_done, netcdf_writes); + } + delete [] elemNeighbors; + + int *elemBoundaries = new int[numElemPerPart[3]*4]; + for (unsigned int z = 0; z < numPartitions[2]; z++) { + for (unsigned int y = 0; y < numPartitions[1]; y++) { + for (unsigned int x = 0; x < numPartitions[0]; x++) { + memset(elemBoundaries, 0, sizeof(int)*numElemPerPart[3]*4); + + if (x == 0) { // first partition in x dimension + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + int odd = (zz+yy) % 2; + if (odd) { + elemBoundaries[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20] = boundaryMinx; + elemBoundaries[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 10] = boundaryMinx; + } else { + elemBoundaries[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 2] = boundaryMinx; + elemBoundaries[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 12] = boundaryMinx; + } + } + } + } + if (x == numPartitions[0]-1) { // last partition in x dimension + + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + int odd = (zz+yy+1) % 2; + if (odd) { + elemBoundaries[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 7] = boundaryMaxx; + elemBoundaries[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 12] = boundaryMaxx; + } else { + elemBoundaries[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 6] = boundaryMaxx; + elemBoundaries[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 9] = boundaryMaxx; + } + } + } + } + if (y == 0) { // first partition in y dimension + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (zz+xx) % 2; + if (odd) { + elemBoundaries[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 6] = boundaryMiny; + elemBoundaries[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 9] = boundaryMiny; + } else { + elemBoundaries[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 1] = boundaryMiny; + elemBoundaries[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 10] = boundaryMiny; + } + } + } + } + if (y == numPartitions[1]-1) { // last partition in y dimension + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (zz+xx+1) % 2; + if (odd) { + elemBoundaries[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 3] = boundaryMaxy; + elemBoundaries[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 14] = boundaryMaxy; + } else { + elemBoundaries[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 7] = boundaryMaxy; + elemBoundaries[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 13] = boundaryMaxy; + } + } + } + } + if (z == 0) { // first partition in z dimension + #pragma omp parallel for + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (yy+xx) % 2; + if (odd) { + elemBoundaries[(yy*numCubesPerPart[0]+xx) * 20 + 1] = boundaryMinz; + elemBoundaries[(yy*numCubesPerPart[0]+xx) * 20 + 5] = boundaryMinz; + } else { + elemBoundaries[(yy*numCubesPerPart[0]+xx) * 20] = boundaryMinz; + elemBoundaries[(yy*numCubesPerPart[0]+xx) * 20 + 5] = boundaryMinz; + } + } + } + } + if (z == numPartitions[2]-1) { // last partition in z dimension + #pragma omp parallel for + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { +// int odd = (yy+xx+1) % 2; +// if (odd) { + elemBoundaries[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 11] = boundaryMaxz; + elemBoundaries[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 15] = boundaryMaxz; +// } else { +// elemBoundaries[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 11] = boundaryMaxz; +// elemBoundaries[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 15] = boundaryMaxz; +// } + } + } + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, 0, 0}; + size_t count[3] = {1, numElemPerPart[3], 4}; + checkNcError(nc_put_vara_int(ncFile, ncVarElemBoundaries, start, count, elemBoundaries)); + writes_done++; loadBar(writes_done, netcdf_writes); + } + } + } + delete [] elemBoundaries; + + int *elemNeighborSides = new int[numElemPerPart[3]*4]; + int *elemNeighborSidesDef = new int[numElemPerPart[3]*4]; + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (zz+yy+xx) % 2; + unsigned int c = ((zz*numCubesPerPart[1] + yy)*numCubesPerPart[0] + xx) * 20; + memcpy(&elemNeighborSidesDef[c], TET_SIDE_NEIGHBORS[odd], sizeof(int)*20); + } + } + } + + for (unsigned int z = 0; z < numPartitions[2]; z++) { + for (unsigned int y = 0; y < numPartitions[1]; y++) { + for (unsigned int x = 0; x < numPartitions[0]; x++) { + memcpy(elemNeighborSides, elemNeighborSidesDef, sizeof(int)*numElemPerPart[3]*4); + + if (boundaryMinx != 6 && x == 0) { //first partition in x dimension + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + int odd = (zz+yy) % 2; + if (odd) { + elemNeighborSides[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20] = 0; + elemNeighborSides[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 10] = 0; + } else { + elemNeighborSides[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 2] = 0; + elemNeighborSides[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 12] = 0; + } + } + } + } + if (boundaryMaxx != 6 && x == numPartitions[0]-1) { // last partition in x dimension + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + int odd = (zz+yy+1) % 2; + if (odd) { + elemNeighborSides[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 7] = 0; + elemNeighborSides[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 12] = 0; + } else { + elemNeighborSides[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 6] = 0; + elemNeighborSides[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 9] = 0; + } + } + } + } + if (boundaryMiny != 6 && y == 0) { // first partition in y dimension + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (zz+xx) % 2; + if (odd) { + elemNeighborSides[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 6] = 0; + elemNeighborSides[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 9] = 0; + } else { + elemNeighborSides[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 1] = 0; + elemNeighborSides[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 10] = 0; + } + } + } + } + if (boundaryMaxy != 6 && y == numPartitions[1]-1) { // last partition in y dimension + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (zz+xx+1) % 2; + if (odd) { + elemNeighborSides[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 3] = 0; + elemNeighborSides[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 14] = 0; + } else { + elemNeighborSides[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 7] = 0; + elemNeighborSides[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 13] = 0; + } + } + } + } + if (boundaryMinz != 6 && z == 0) { // first partition in z dimension + #pragma omp parallel for + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (yy+xx) % 2; + if (odd) { + elemNeighborSides[(yy*numCubesPerPart[0]+xx) * 20 + 1] = 0; + elemNeighborSides[(yy*numCubesPerPart[0]+xx) * 20 + 5] = 0; + } else { + elemNeighborSides[(yy*numCubesPerPart[0]+xx) * 20] = 0; + elemNeighborSides[(yy*numCubesPerPart[0]+xx) * 20 + 5] = 0; + } + } + } + } + if (boundaryMaxz != 6 && z == numPartitions[2]-1) { // last partition in z dimension + #pragma omp parallel for + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { +// int odd = (yy+xx+1) % 2; +// if (odd) { + elemNeighborSides[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 11] = 0; + elemNeighborSides[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 15] = 0; +// } else { +// elemSideNeighbors[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 11] = 0; +// elemSideNeighbors[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 15] = 0; +// } + } + } + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, 0, 0}; + size_t count[3] = {1, numElemPerPart[3], 4}; + checkNcError(nc_put_vara_int(ncFile, ncVarElemNeighborSides, start, count, elemNeighborSides)); + writes_done++; loadBar(writes_done, netcdf_writes); + } + } + } + delete [] elemNeighborSides; + delete [] elemNeighborSidesDef; + + int *elemSideOrientations = new int[numElemPerPart[3]*4]; + int *elemSideOrientationsDef = new int[numElemPerPart[3]*4]; + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (zz+yy+xx) % 2; + unsigned int c = ((zz*numCubesPerPart[1] + yy)*numCubesPerPart[0] + xx) * 20; + memcpy(&elemSideOrientationsDef[c], TET_SIDE_ORIENTATIONS[odd], sizeof(int)*20); + } + } + } + + for (unsigned int z = 0; z < numPartitions[2]; z++) { + for (unsigned int y = 0; y < numPartitions[1]; y++) { + for (unsigned int x = 0; x < numPartitions[0]; x++) { + memcpy(elemSideOrientations, elemSideOrientationsDef, sizeof(int)*numElemPerPart[3]*4); + + if (boundaryMinx != 6 && x == 0) { // first partition in x dimension + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + int odd = (zz+yy) % 2; + if (odd) { + elemSideOrientations[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20] = 0; + elemSideOrientations[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 10] = 0; + } else { + elemSideOrientations[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 2] = 0; + elemSideOrientations[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 12] = 0; + } + } + } + } + if (boundaryMaxx != 6 && x == numPartitions[0]-1) { // last partition in x dimension + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + int odd = (zz+yy+1) % 2; + if (odd) { + elemSideOrientations[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 7] = 0; + elemSideOrientations[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 12] = 0; + } else { + elemSideOrientations[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 6] = 0; + elemSideOrientations[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 9] = 0; + } + } + } + } + // There are zero anyway +// if (boundaryMiny != 6 && y == 0) { // first partition in y dimension +// #pragma omp parallel for +// for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { +// for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { +// int odd = (zz+xx) % 2; +// if (odd) { +// elemSideOrientations[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 6] = 0; +// elemSideOrientations[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 9] = 0; +// } else { +// elemSideOrientations[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 1] = 0; +// elemSideOrientations[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 10] = 0; +// } +// } +// } +// } +// if (boundaryMaxy != 6 && y == numPartitions[1]-1) { // last partition in y dimension +// #pragma omp parallel for +// for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { +// for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { +// int odd = (zz+xx+1) % 2; +// if (odd) { +// elemSideOrientations[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 3] = 0; +// elemSideOrientations[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 14] = 0; +// } else { +// elemSideOrientations[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 7] = 0; +// elemSideOrientations[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 13] = 0; +// } +// } +// } +// } + if (boundaryMinz != 6 && z == 0) { // first partition in z dimension + #pragma omp parallel for + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (yy+xx) % 2; + if (odd) { + elemSideOrientations[(yy*numCubesPerPart[0]+xx) * 20 + 1] = 0; + elemSideOrientations[(yy*numCubesPerPart[0]+xx) * 20 + 5] = 0; + } else { + elemSideOrientations[(yy*numCubesPerPart[0]+xx) * 20] = 0; + elemSideOrientations[(yy*numCubesPerPart[0]+xx) * 20 + 5] = 0; + } + } + } + } + if (boundaryMaxz != 6 && z == numPartitions[2]-1) { // last partition in z dimension + #pragma omp parallel for + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { +// int odd = (yy+xx+1) % 2; +// if (odd) { + elemSideOrientations[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 11] = 0; + elemSideOrientations[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 15] = 0; +// } else { +// elemSideOrientations[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 11] = 0; +// elemSideOrientations[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 15] = 0; +// } + } + } + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, 0, 0}; + size_t count[3] = {1, numElemPerPart[3], 4}; + checkNcError(nc_put_vara_int(ncFile, ncVarElemSideOrientations, start, count, elemSideOrientations)); + writes_done++; loadBar(writes_done, netcdf_writes); + } + } + } + delete [] elemSideOrientations; + delete [] elemSideOrientationsDef; + + int *elemNeighborRanks = new int[numElemPerPart[3]*4]; + + for (unsigned int z = 0; z < numPartitions[2]; z++) { + for (unsigned int y = 0; y < numPartitions[1]; y++) { + for (unsigned int x = 0; x < numPartitions[0]; x++) { + int myrank = (z*numPartitions[1] + y)*numPartitions[0] + x; + + std::fill(elemNeighborRanks, elemNeighborRanks+numElemPerPart[3]*4, myrank); + + if ((boundaryMinx == 6 && numPartitions[0] > 1) || x != 0) { // first partition in x dimension + int rank = (z*numPartitions[1] + y)*numPartitions[0] + (x-1+numPartitions[0])%numPartitions[0]; + + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + int odd = (zz+yy) % 2; + if (odd) { + elemNeighborRanks[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20] = rank; + elemNeighborRanks[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 10] = rank; + } else { + elemNeighborRanks[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 2] = rank; + elemNeighborRanks[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 12] = rank; + } + } + } + } + if ((boundaryMaxx == 6 && numPartitions[0] > 1) || x != numPartitions[0]-1) { // last partition in x dimension + int rank = (z*numPartitions[1] + y)*numPartitions[0] + (x+1)%numPartitions[0]; + + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + int odd = (zz+yy+1) % 2; + if (odd) { + elemNeighborRanks[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 7] = rank; + elemNeighborRanks[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 12] = rank; + } else { + elemNeighborRanks[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 6] = rank; + elemNeighborRanks[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 9] = rank; + } + } + } + } + if ((boundaryMiny == 6 && numPartitions[1] > 1) || y != 0) { // first partition in y dimension + int rank = (z*numPartitions[1] + (y-1+numPartitions[1])%numPartitions[1])*numPartitions[0] + x; + + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (zz+xx) % 2; + if (odd) { + elemNeighborRanks[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 6] = rank; + elemNeighborRanks[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 9] = rank; + } else { + elemNeighborRanks[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 1] = rank; + elemNeighborRanks[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 10] = rank; + } + } + } + } + if ((boundaryMaxy == 6 && numPartitions[1] > 1) || y != numPartitions[1]-1) { // last partition in y dimension + int rank = (z*numPartitions[1] + (y+1)%numPartitions[1])*numPartitions[0] + x; + + #pragma omp parallel for + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (zz+xx+1) % 2; + if (odd) { + elemNeighborRanks[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 3] = rank; + elemNeighborRanks[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 14] = rank; + } else { + elemNeighborRanks[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 7] = rank; + elemNeighborRanks[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 13] = rank; + } + } + } + } + if ((boundaryMinz == 6 && numPartitions[2] > 1) || z != 0) { // first partition in z dimension + int rank = (((z-1+numPartitions[2])%numPartitions[2])*numPartitions[1] + y)*numPartitions[0] + x; + + #pragma omp parallel for + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (yy+xx) % 2; + if (odd) { + elemNeighborRanks[(yy*numCubesPerPart[0]+xx) * 20 + 1] = rank; + elemNeighborRanks[(yy*numCubesPerPart[0]+xx) * 20 + 5] = rank; + } else { + elemNeighborRanks[(yy*numCubesPerPart[0]+xx) * 20] = rank; + elemNeighborRanks[(yy*numCubesPerPart[0]+xx) * 20 + 5] = rank; + } + } + } + } + if ((boundaryMaxz == 6 && numPartitions[2] > 1) || z != numPartitions[2]-1) { // last partition in z dimension + int rank = (((z+1)%numPartitions[2])*numPartitions[1] + y)*numPartitions[0] + x; + + #pragma omp parallel for + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { +// int odd = (yy+xx+1) % 2; +// if (odd) { + elemNeighborRanks[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 11] = rank; + elemNeighborRanks[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 15] = rank; +// } else { +// elemNeighborRanks[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 11] = rank; +// elemNeighborRanks[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 15] = rank; +// } + } + } + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, 0, 0}; + size_t count[3] = {1, numElemPerPart[3], 4}; + checkNcError(nc_put_vara_int(ncFile, ncVarElemNeighborRanks, start, count, elemNeighborRanks)); + writes_done++; loadBar(writes_done, netcdf_writes); + } + } + } + delete [] elemNeighborRanks; + + int *elemMPIIndices = new int[numElemPerPart[3]*4]; + int *bndLocalIds = new int[*std::max_element(numBndElements, numBndElements+3)]; + + for (unsigned int z = 0; z < numPartitions[2]; z++) { + for (unsigned int y = 0; y < numPartitions[1]; y++) { + for (unsigned int x = 0; x < numPartitions[0]; x++) { + memset(elemMPIIndices, 0, sizeof(int)*numElemPerPart[3]*4); + + unsigned int bndSize = 0; + + if ((boundaryMinz == 6 && numPartitions[2] > 1) || z != 0) { // first partition in z dimension + int nextMPIIndex = 0; + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (yy+xx) % 2; + if (odd) { + bndLocalIds[nextMPIIndex] = (yy*numCubesPerPart[0]+xx) * 5 + 1; + elemMPIIndices[(yy*numCubesPerPart[0]+xx) * 20 + 5] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = (yy*numCubesPerPart[0]+xx) * 5; + elemMPIIndices[(yy*numCubesPerPart[0]+xx) * 20 + 1] = nextMPIIndex++; + } else { + bndLocalIds[nextMPIIndex] = (yy*numCubesPerPart[0]+xx) * 5; + elemMPIIndices[(yy*numCubesPerPart[0]+xx) * 20] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = (yy*numCubesPerPart[0]+xx) * 5 + 1; + elemMPIIndices[(yy*numCubesPerPart[0]+xx) * 20 + 5] = nextMPIIndex++; + } + } + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, bndSize, 0u}; + size_t count[3] = {1, 1, static_cast(nextMPIIndex)}; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemSize, start, &nextMPIIndex)); + int rank = (((z-1+numPartitions[2])%numPartitions[2])*numPartitions[1] + y)*numPartitions[0] + x; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemRank, start, &rank)); + checkNcError(nc_put_vara_int(ncFile, ncVarBndElemLocalIds, start, count, bndLocalIds)); + + bndSize++; + } + if ((boundaryMiny == 6 && numPartitions[1] > 1) || y != 0) { // first partition in y dimension + int nextMPIIndex = 0; + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (zz+xx) % 2; + if (odd) { + bndLocalIds[nextMPIIndex] = (zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 5 + 1; + elemMPIIndices[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 6] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = (zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 5 + 2; + elemMPIIndices[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 9] = nextMPIIndex++; + } else { + bndLocalIds[nextMPIIndex] = (zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 5; + elemMPIIndices[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 1] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = (zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 5 + 2; + elemMPIIndices[(zz*numCubesPerPart[1]*numCubesPerPart[0]+xx) * 20 + 10] = nextMPIIndex++; + } + } + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, bndSize, 0u}; + size_t count[3] = {1, 1, static_cast(nextMPIIndex)}; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemSize, start, &nextMPIIndex)); + int rank = (z*numPartitions[1] + (y-1+numPartitions[1])%numPartitions[1])*numPartitions[0] + x; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemRank, start, &rank)); + checkNcError(nc_put_vara_int(ncFile, ncVarBndElemLocalIds, start, count, bndLocalIds)); + + bndSize++; + } + if ((boundaryMinx == 6 && numPartitions[0] > 1) || x != 0) { // first partition in x dimension + int nextMPIIndex = 0; + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + int odd = (zz+yy) % 2; + if (odd) { + bndLocalIds[nextMPIIndex] = (zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 5; + elemMPIIndices[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = (zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 5 + 2; + elemMPIIndices[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 10] = nextMPIIndex++; + } else { + bndLocalIds[nextMPIIndex] = (zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 5; + elemMPIIndices[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 2] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = (zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 5 + 3; + elemMPIIndices[(zz*numCubesPerPart[1]+yy)*numCubesPerPart[0] * 20 + 12] = nextMPIIndex++; + } + } + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, bndSize, 0u}; + size_t count[3] = {1, 1, static_cast(nextMPIIndex)}; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemSize, start, &nextMPIIndex)); + int rank = (z*numPartitions[1] + y)*numPartitions[0] + (x-1+numPartitions[0])%numPartitions[0]; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemRank, start, &rank)); + checkNcError(nc_put_vara_int(ncFile, ncVarBndElemLocalIds, start, count, bndLocalIds)); + + bndSize++; + } + if ((boundaryMaxx == 6 && numPartitions[0] > 1) || x != numPartitions[0]-1) { // last partition in x dimension + int nextMPIIndex = 0; + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + int odd = (zz+yy+1) % 2; + if (odd) { + bndLocalIds[nextMPIIndex] = ((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 5 + 1; + elemMPIIndices[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 7] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = ((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 5 + 3; + elemMPIIndices[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 12] = nextMPIIndex++; + } else { + bndLocalIds[nextMPIIndex] = ((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 5 + 1; + elemMPIIndices[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 6] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = ((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 5 + 2; + elemMPIIndices[((zz*numCubesPerPart[1]+yy)*numCubesPerPart[0]+numCubesPerPart[0]-1) * 20 + 9] = nextMPIIndex++; + } + } + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, bndSize, 0}; + size_t count[3] = {1, 1, static_cast(nextMPIIndex)}; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemSize, start, &nextMPIIndex)); + int rank = (z*numPartitions[1] + y)*numPartitions[0] + (x+1)%numPartitions[0]; + rank = (rank + numPartitions[3]) % numPartitions[3]; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemRank, start, &rank)); + checkNcError(nc_put_vara_int(ncFile, ncVarBndElemLocalIds, start, count, bndLocalIds)); + + bndSize++; + } + if ((boundaryMaxy == 6 && numPartitions[1] > 1) || y != numPartitions[1]-1) { // last partition in y dimension + int nextMPIIndex = 0; + for (unsigned int zz = 0; zz < numCubesPerPart[2]; zz++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { + int odd = (zz+xx+1) % 2; + if (odd) { + bndLocalIds[nextMPIIndex] = ((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 5; + elemMPIIndices[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 3] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = ((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 5 + 3; + elemMPIIndices[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 14] = nextMPIIndex++; + } else { + bndLocalIds[nextMPIIndex] = ((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 5 + 1; + elemMPIIndices[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 7] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = ((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 5 + 3; + elemMPIIndices[((zz*numCubesPerPart[1]+numCubesPerPart[1]-1)*numCubesPerPart[0]+xx) * 20 + 13] = nextMPIIndex++; + } + } + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, bndSize, 0}; + size_t count[3] = {1, 1, static_cast(nextMPIIndex)}; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemSize, start, &nextMPIIndex)); + int rank = (z*numPartitions[1] + (y+1)%numPartitions[1])*numPartitions[0] + x; + rank = (rank + numPartitions[3]) % numPartitions[3]; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemRank, start, &rank)); + checkNcError(nc_put_vara_int(ncFile, ncVarBndElemLocalIds, start, count, bndLocalIds)); + + bndSize++; + } + if ((boundaryMaxz == 6 && numPartitions[2] > 1) || z != numPartitions[2]-1) { // last partition in z dimension + int nextMPIIndex = 0; + for (unsigned int yy = 0; yy < numCubesPerPart[1]; yy++) { + for (unsigned int xx = 0; xx < numCubesPerPart[0]; xx++) { +// int odd = (yy+xx+1) % 2; +// if (odd) { + bndLocalIds[nextMPIIndex] = (((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 5 + 2; + elemMPIIndices[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 11] = nextMPIIndex++; + bndLocalIds[nextMPIIndex] = (((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 5 + 3; + elemMPIIndices[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 15] = nextMPIIndex++; +// } else { +// bndLocalIds[nextMPIIndex] = (((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 5 + 2; +// elemMPIIndices[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 11] = nextMPIIndex++; +// bndLocalIds[nextMPIIndex] = (((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 5 + 3; +// elemMPIIndices[(((numCubesPerPart[2]-1)*numCubesPerPart[1]+yy)*numCubesPerPart[0]+xx) * 20 + 15] = nextMPIIndex++; +// } + } + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, bndSize, 0u}; + size_t count[3] = {1, 1, static_cast(nextMPIIndex)}; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemSize, start, &nextMPIIndex)); + int rank = (((z+1)%numPartitions[2])*numPartitions[1] + y)*numPartitions[0] + x; + rank = (rank + numPartitions[3]) % numPartitions[3]; + checkNcError(nc_put_var1_int(ncFile, ncVarBndElemRank, start, &rank)); + checkNcError(nc_put_vara_int(ncFile, ncVarBndElemLocalIds, start, count, bndLocalIds)); + + bndSize++; + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, 0, 0}; + size_t count[3] = {1, numElemPerPart[3], 4}; + checkNcError(nc_put_vara_int(ncFile, ncVarElemMPIIndices, start, count, elemMPIIndices)); + + checkNcError(nc_put_var1_uint(ncFile, ncVarBndSize, &start[0], &bndSize)); + writes_done++; loadBar(writes_done, netcdf_writes); + } + } + } + delete [] elemMPIIndices; + delete [] bndLocalIds; + + // Set material zone to 1 + int *elemGroup = new int[numElemPerPart[3]]; + std::fill(elemGroup, elemGroup + numElemPerPart[3], 1); + for (unsigned int x = 0; x < numPartitions[3]; x++) { + size_t start[2] = {x, 0}; + size_t count[2] = {1, numElemPerPart[3]}; + checkNcError(nc_put_vara_int(ncFile, ncVarElemGroup, start, count, elemGroup)); + } + delete[] elemGroup; + + // Vertices + std::map uniqueVertices; + transform(vertexMap.begin(), vertexMap.end(), + inserter(uniqueVertices, uniqueVertices.begin()), + flip_pair); + + int *vrtxSize = new int[numPartitions[3]]; + std::fill(vrtxSize, vrtxSize+numPartitions[3], uniqueVertices.size()); + checkNcError(nc_put_var_int(ncFile, ncVarVrtxSize, vrtxSize)); + delete [] vrtxSize; + writes_done++; loadBar(writes_done, netcdf_writes); + + double *vrtxCoords = new double[uniqueVertices.size()*3]; + + + double scale = args.getArgument("scale", 100.0); + double scaleX = args.getArgument("sx", scale); + double scaleY = args.getArgument("sy", scale); + double scaleZ = args.getArgument("sz", scale); + double halfWidthX = scaleX / 2.0; + double halfWidthY = scaleY / 2.0; + double halfWidthZ = scaleZ / 2.0; + + double tx = args.getArgument("tx", 0.0); + double ty = args.getArgument("ty", 0.0); + double tz = args.getArgument("tz", 0.0); + + for (unsigned int z = 0; z < numPartitions[2]; z++) { + for (unsigned int y = 0; y < numPartitions[1]; y++) { + for (unsigned int x = 0; x < numPartitions[0]; x++) { + + #pragma omp parallel for + for (unsigned int i = 0; i < uniqueVertices.size(); i++) { + vrtxCoords[i*3] = static_cast(uniqueVertices.at(i).v[0] + x*numCubesPerPart[0])/static_cast(numCubes[0])*scaleX - halfWidthX + tx; + vrtxCoords[i*3+1] = static_cast(uniqueVertices.at(i).v[1] + y*numCubesPerPart[1])/static_cast(numCubes[1])*scaleY - halfWidthY + ty; + vrtxCoords[i*3+2] = static_cast(uniqueVertices.at(i).v[2] + z*numCubesPerPart[2])/static_cast(numCubes[2])*scaleZ - halfWidthZ + tz; + } + + size_t start[3] = {(z*numPartitions[1] + y)*numPartitions[0] + x, 0, 0}; + size_t count[3] = {1, uniqueVertices.size(), 3}; + checkNcError(nc_put_vara_double(ncFile, ncVarVrtxCoords, start, count, vrtxCoords)); + writes_done++; loadBar(writes_done, netcdf_writes); + } + } + } + + delete [] vrtxCoords; + + checkNcError(nc_close(ncFile)); + + logInfo() << "Finished"; +} diff --git a/cube_c/src/utils/args.h b/cube_c/src/utils/args.h new file mode 100644 index 0000000..03cee2a --- /dev/null +++ b/cube_c/src/utils/args.h @@ -0,0 +1,443 @@ +/** + * @file + * This file is part of SeisSol. + * + * @author Sebastian Rettenberger (rettenbs AT in.tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger,_M.Sc.) + * + * @section LICENSE + * Copyright (c) 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 + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * + * @section DESCRIPTION + * TODO place this in the submodule folder + */ + +#ifndef UTILS_ARGS_H +#define UTILS_ARGS_H + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace utils +{ + +/** + * Parses command line arguments + * @todo comment functions + */ +class Args +{ +private: + struct optionInfo { + std::string longOption; // We need a copy here to get the const char* correct + /** Name of the value in the help command */ + std::string value; + std::string description; + bool required; + }; + + struct additionalOptionInfo { + std::string name; + std::string description; + bool required; + }; + + /** + * Convert a long option into an "argument" that is shown in the help message + */ + struct valueConvert { + void operator()(char& c) + { + c = toupper(static_cast(c)); + switch (c) { + case '-': + c = '_'; + break; + } + } + }; + + /** Program description (can be empty) */ + const std::string m_description; + /** Automatically add help option */ + const bool m_addHelp; + + /** The command line arguments */ + std::vector m_options; + /** + * Additional information for the command line arguments + * required to generate help information + */ + std::vector m_optionInfo; + + std::vector m_additionalOptionInfo; + + /** Maps from short option to index in m_options */ + std::map m_short2option; + + /** Contains the arguments after parse was called */ + std::map m_arguments; + + /** + * Contains additional arguments after parse was called + * @todo Find a better name + */ + std::map m_additionalArguments; + +public: + enum Argument + { + Required = required_argument, + No = no_argument, + Optional = optional_argument + }; + + enum Result + { + Success = 0, + Error, + /** Help message printed */ + Help + }; + +public: + Args(const std::string &description = "", bool addHelp = true) + : m_description(description), + m_addHelp(addHelp) + { + } + + void addOption(const std::string &longOption, + char shortOption = 0, + const std::string &description = "", + Argument argument = Required, + bool required = true) + { + std::string value; + + if (shortOption) + m_short2option[shortOption] = m_options.size(); + + if (argument != No) { + value = longOption; + for_each(value.begin(), value.end(), valueConvert()); + } + + struct optionInfo i = {longOption, value, description, required}; + m_optionInfo.push_back(i); + + struct option o = {m_optionInfo.back().longOption.c_str(), argument, 0, shortOption}; + m_options.push_back(o); + } + + void addAdditionalOption(const std::string &name, + const std::string &description = "", + bool required = true) + { + if (!m_additionalOptionInfo.empty()) { + if (required && !m_additionalOptionInfo.back().required) + // After one optional argument there can only be more optional arguments + return; + } + + struct additionalOptionInfo i = {name, description, required}; + m_additionalOptionInfo.push_back(i); + } + + /** + * @return True of options are successfully parsed, false otherwise + */ + Result parse(int argc, char* const* argv, bool printHelp = true) + { + if (m_addHelp) + addOption("help", 'h', "Show this help message", No, false); + + std::ostringstream shortOptions; + for (std::vector::const_iterator i = m_options.begin(); + i != m_options.end(); i++) { + if (i->val != 0) { + shortOptions << static_cast(i->val); + switch (i->has_arg) + { + case required_argument: + shortOptions << ':'; + break; + case optional_argument: + shortOptions << "::"; + break; + } + } + } + + // Add null option + struct option o = {0, 0, 0, 0}; + m_options.push_back(o); + + // Update const char* in m_options + for (size_t i = 0; i < m_optionInfo.size(); i++) + m_options[i].name = m_optionInfo[i].longOption.c_str(); + + while (true) { + int optionIndex = 0; + + int c = getopt_long(argc, argv, shortOptions.str().c_str(), + &m_options[0], &optionIndex); + + if (c < 0) + break; + + switch (c) { + case '?': + if (printHelp) + helpMessage(argv[0], std::cerr); + return Error; + case 0: + // Nothing to do + break; + default: + optionIndex = m_short2option.at(c); + } + + if (optarg == 0L) + m_arguments[m_options[optionIndex].name] = ""; + else + m_arguments[m_options[optionIndex].name] = optarg; + } + + if (m_addHelp && isSet("help")) { + if (printHelp) + helpMessage(argv[0]); + return Help; + } + + for (std::vector::const_iterator i = m_optionInfo.begin(); + i != m_optionInfo.end(); i++) { + if (i->required && !isSet(i->longOption)) { + if (printHelp) { + std::cerr << argv[0] << ": option --" << i->longOption << " is required" << std::endl; + helpMessage(argv[0], std::cerr); + } + return Error; + } + } + + // Parse additional options and check if all required options are set + int i; + for (i = 0; i < argc-optind; i++) + m_additionalArguments[m_additionalOptionInfo[i].name] = argv[i+optind]; + if (static_cast(i) < m_additionalOptionInfo.size()) { + if (m_additionalOptionInfo[i].required) { + if (printHelp) { + std::cerr << argv[0] << ": option <" << m_additionalOptionInfo[i].name << "> is required" << std::endl; + helpMessage(argv[0], std::cerr); + } + return Error; + } + } + + return Success; + } + + bool isSet(const std::string &option) const + { + return m_arguments.find(option) != m_arguments.end(); + } + + bool isSetAdditional(const std::string &option) const + { + return m_additionalArguments.find(option) != m_additionalArguments.end(); + } + + template + T getArgument(const std::string &option) + { + std::istringstream ss(m_arguments.at(option)); + + T result; + ss >> result; + + return result; + } + + template + T getArgument(const std::string &option, T defaultArgument) + { + if (!isSet(option)) + return defaultArgument; + + return getArgument(option); + } + + + template + T getAdditionalArgument(const std::string &option) + { + std::istringstream ss(m_additionalArguments.at(option)); + + T result; + ss >> result; + + return result; + } + + template + T getAdditionalArgument(const std::string &option, T defaultArgument) + { + if (!isSetAdditional(option)) + return defaultArgument; + + return getAdditionalArgument(option); + } + + void helpMessage(const char* prog, std::ostream &out = std::cout) + { + // First line with all short options + out << "Usage: " << prog; + for (size_t i = 0; i < m_options.size()-1; i++) { + out << ' '; + + if (!m_optionInfo[i].required) + out << '['; + + if (m_options[i].val != 0) + out << '-' << static_cast(m_options[i].val); + else + out << "--" << m_options[i].name; + + argumentInfo(i, out); + + if (!m_optionInfo[i].required) + out << ']'; + } + for (size_t i = 0; i < m_additionalOptionInfo.size(); i++) { + out << ' '; + + if (!m_additionalOptionInfo[i].required) + out << '['; + + out << '<' << m_additionalOptionInfo[i].name << '>'; + + if (!m_additionalOptionInfo[i].required) + out << ']'; + + } + out << std::endl; + + // General program description + if (!m_description.empty()) + out << std::endl << m_description << std::endl; + + // Arguments + out << std::endl << "arguments:" << std::endl; + for (size_t i = 0; i < m_additionalOptionInfo.size(); i++) { + out << " <" << m_additionalOptionInfo[i].name << '>'; + + // Number of characters used for the option + size_t length = 4 + m_additionalOptionInfo[i].name.size(); + + if (length >= 30) { + out << std::endl; + out << std::setw(30) << ' '; + } else + out << std::setw(30-length) << ' '; + + out << m_additionalOptionInfo[i].description << std::endl; + } + + // Optional arguments + out << std::endl << "optional arguments:" << std::endl; + for (size_t i = 0; i < m_options.size()-1; i++) { + out << " "; + + // Number of characters used for the option + size_t length = 2; + + if (m_options[i].val != 0) { + out << '-' << static_cast(m_options[i].val); + out << ", "; + length += 4; + } + + out << "--" << m_options[i].name; + length += m_optionInfo[i].longOption.size() + 2; + length += argumentInfo(i, out); + + if (length >= 30) { + out << std::endl; + out << std::setw(30) << ' '; + } else + out << std::setw(30-length) << ' '; + + out << m_optionInfo[i].description << std::endl; + } + } + +private: + /** + * Writes the argument information to out + * + * @param i The index of the option for which the argument should be generated + * @return The number if characters written + */ + size_t argumentInfo(size_t i, std::ostream &out) + { + switch (m_options[i].has_arg) { + case required_argument: + out << ' ' << m_optionInfo[i].value; + return m_optionInfo[i].value.size() + 1; + case optional_argument: + out << " [" << m_optionInfo[i].value << ']'; + return m_optionInfo[i].value.size() + 3; + } + + return 0; + } +}; + +template<> inline +std::string utils::Args::getArgument(const std::string &option) +{ + return m_arguments.at(option); +} + +template<> inline +std::string utils::Args::getAdditionalArgument(const std::string &option) +{ + return m_additionalArguments.at(option); +} + +} + +#endif // UTILS_ARGS_H diff --git a/cube_c/submodules b/cube_c/submodules new file mode 120000 index 0000000..8c556f9 --- /dev/null +++ b/cube_c/submodules @@ -0,0 +1 @@ +/home/david/Dokumente/uni/phd/software/meshing/Meshing/submodules \ No newline at end of file diff --git a/gmsh2gambit/CMakeLists.txt b/gmsh2gambit/CMakeLists.txt new file mode 100644 index 0000000..4e786ba --- /dev/null +++ b/gmsh2gambit/CMakeLists.txt @@ -0,0 +1,42 @@ +cmake_minimum_required(VERSION 3.10) + +# set the project name +project(gmsh2gambit) + +# specify the C++ standard (11 is enough for now, according to the SCons file) +set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD_REQUIRED True) + +cmake_policy(SET CMP0074 NEW) + + +add_executable(gmsh2gambit + src/main.cpp + src/GMSH.cpp + src/GAMBITWriter.cpp) + +#logging +set(LOG_LEVEL "warning" CACHE STRING "Log level for the code") +set(LOG_LEVEL_OPTIONS "debug" "info" "warning" "error") +set_property(CACHE LOG_LEVEL PROPERTY STRINGS ${LOG_LEVEL_OPTIONS}) +if("${LOG_LEVEL}" STREQUAL "debug") + target_compile_definitions(gmsh2gambit PUBLIC LOG_LEVEL=3) +elseif("${LOG_LEVEL}" STREQUAL "info") + target_compile_definitions(gmsh2gambit PUBLIC LOG_LEVEL=2) +elseif("${LOG_LEVEL}" STREQUAL "warning") + target_compile_definitions(gmsh2gambit PUBLIC LOG_LEVEL=1) +elseif("${LOG_LEVEL}" STREQUAL "error") + target_compile_definitions(gmsh2gambit PUBLIC LOG_LEVEL=0) +endif() + +#add some compiler specific flags +if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") + # using GCC + target_compile_options(gmsh2gambit PUBLIC -pedantic $<$,$>:-Wall -Wextra -Wno-unused-parameter -Wno-unknown-pragmas>) +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") + target_compile_options(gmsh2gambit PUBLIC -pedantic $<$,$>:-Wall -w3 -diag-disable:remark>) +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -pedantic") +endif() + +target_include_directories(gmsh2gambit PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/submodules) diff --git a/gmsh2gambit/README.md b/gmsh2gambit/README.md new file mode 100644 index 0000000..13e21c6 --- /dev/null +++ b/gmsh2gambit/README.md @@ -0,0 +1,26 @@ +gmsh2gambit +=============== +Transform msh file (gmsh msh 2.0 mesh file) to neu file (neutral mesh file) +From gmsh 4.0, the option '-format neu' of gmsh should be used instead. + +Build +=============== +Using CMake (preferred) +``` +cmake +make -j 3 +``` +(this assumes that UNIX Makefiles are the default build system for your CMake installation) + +Using scons (legacy) +``` +CC=gcc CXX=g++ scons +``` +Or, replace `gcc` and `g++` by another compiler (e.g. `icc` and `icpc` for the legacy Intel compilers). +Note that the specification of `CC` and `CXX` is mandatory for the scons script to work. + +Usage +=============== +gmsh2gambit -i test.msh -o test.neu + + diff --git a/gmsh2gambit/SConstruct b/gmsh2gambit/SConstruct new file mode 100644 index 0000000..77c6654 --- /dev/null +++ b/gmsh2gambit/SConstruct @@ -0,0 +1,96 @@ +#! /usr/bin/python +## +# @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.) +# +# @section LICENSE +# Copyright (c) 2015, 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 +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# + +import os + +# +# set possible variables +# +vars = Variables() + +vars.AddVariables( + PathVariable( 'buildDir', 'where to build the code', 'build', PathVariable.PathIsDirCreate ), + + EnumVariable( 'compileMode', 'mode of the compilation', 'release', + allowed_values=('debug', 'release') + ) +) + +# set environment +env = Environment(variables=vars, CC=os.environ['CC'], CXX=os.environ['CXX']) +env['ENV'] = os.environ + +# generate help text +Help(vars.GenerateHelpText(env)) + +# handle unknown, maybe misspelled variables +unknownVariables = vars.UnknownVariables() + +# exit in the case of unknown variables +if unknownVariables: + raise EnvironmentError("*** The following build variables are unknown: " + str(unknownVariables.keys())) + + +# +# precompiler, compiler and linker flags +# +env.Append(CXXFLAGS=['--std=c++11']) +if env['compileMode'] == 'debug': + env.Append(CXXFLAGS=['-O0','-g']) +elif env['compileMode'] == 'release': + env.Append(CPPDEFINES=['NDEBUG']) + env.Append(CXXFLAGS=['-O2']) + +# add pathname to the list of directories which are search for include +env.Append(CPPPATH=['#/src', '#/submodules']) + +# build directory +env['execDir'] = env['buildDir']+'/bin' +env['buildDir'] = env['buildDir'] + +# get the source files +env.sourceFiles = [] + +Export('env') +SConscript('src/SConscript', variant_dir='#/'+env['buildDir'], src_dir='#/', duplicate=0) +Import('env') + +# build executable +env.Program('#/'+env['execDir']+'/gmsh2gambit', env.sourceFiles) diff --git a/gmsh2gambit/src/GAMBITWriter.cpp b/gmsh2gambit/src/GAMBITWriter.cpp new file mode 100644 index 0000000..2e12860 --- /dev/null +++ b/gmsh2gambit/src/GAMBITWriter.cpp @@ -0,0 +1,43 @@ +/** + * @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.) + * + * @section LICENSE + * Copyright (c) 2015, 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 + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * + * @section DESCRIPTION + **/ + +#include "GAMBITWriter.h" + +template<> unsigned const GambitInfo<2>::Type = 3; +template<> unsigned const GambitInfo<3>::Type = 6; diff --git a/gmsh2gambit/src/GAMBITWriter.h b/gmsh2gambit/src/GAMBITWriter.h new file mode 100644 index 0000000..88643fc --- /dev/null +++ b/gmsh2gambit/src/GAMBITWriter.h @@ -0,0 +1,124 @@ +/** + * @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.) + * + * @section LICENSE + * Copyright (c) 2015, 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 + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * + * @section DESCRIPTION + **/ +#ifndef GAMBITWRITER_H_ +#define GAMBITWRITER_H_ +#include +#include + +#include "GMSH.h" + +template +struct GambitInfo { + static unsigned const Type; +}; + +template +void writeNEU(char const* filename, GMSH const& msh) { + FILE* file; + file = fopen(filename, "w"); + fprintf(file, " CONTROL INFO 2.0.0\n"); + fprintf(file, "** GAMBIT NEUTRAL FILE\n"); + fprintf(file, "Gmsh mesh in GAMBIT neutral file format\n"); + fprintf(file, "PROGRAM: Gambit VERSION: 2.0.0\n"); + time_t rawtime; + struct tm* timeinfo; + time(&rawtime); + timeinfo = localtime(&rawtime); + char datestring[256]; + strftime(datestring, sizeof(datestring), "%c", timeinfo); + fprintf(file, "%s\n", datestring); + fprintf(file, " NUMNP NELEM NGRPS NBSETS NDFCD NDFVL\n"); + fprintf(file, " %9d %9d %9d %9d %9d %9d\n", msh.numVertices, msh.numElements, static_cast(msh.materialGroups.size()), static_cast(msh.boundaryConditions.size()), DIM, DIM); + fprintf(file, "ENDOFSECTION\n"); + + // Vertices + fprintf(file, " NODAL COORDINATES 2.0.0\n"); + for (unsigned vertex = 0; vertex < msh.numVertices; ++vertex) { + fprintf(file, "%10d", vertex+1); + for (unsigned c = 0; c < DIM; ++c) { + fprintf(file, "%20.11e", msh.vertices[vertex][c]); + } + fprintf(file, "\n"); + } + fprintf(file, "ENDOFSECTION\n"); + + // Elements + fprintf(file, " ELEMENTS/CELLS 2.0.0\n"); + for (unsigned element = 0; element < msh.numElements; ++element) { + fprintf(file, "%8d %2d %2d ", element+1, GambitInfo::Type, GMSH::Element::NumNodes); + for (unsigned n = 0; n < GMSH::Element::NumNodes; ++n) { + fprintf(file, "%8d", msh.elements[element].nodes[n]+1); + } + fprintf(file, "\n"); + } + fprintf(file, "ENDOFSECTION\n"); + + // Material groups + for (auto group = msh.materialGroups.cbegin(); group != msh.materialGroups.cend(); ++group) { + fprintf(file, " ELEMENT GROUP 2.0.0\n"); + fprintf(file, "GROUP: %10d ELEMENTS: %10d MATERIAL: %10d NFLAGS: %10d\n", group->first, static_cast(group->second.size()), 2, 1); + fprintf(file, "Material group %d\n", group->first); + fprintf(file, " 0"); + for (unsigned gm = 0; gm < group->second.size(); ++gm) { + if (gm % 10 == 0) { + fprintf(file, "\n"); + } + fprintf(file, "%8d", group->second[gm]+1); + } + fprintf(file, "\n"); + fprintf(file, "ENDOFSECTION\n"); + } + + // Boundary conditions + for (auto boundaryCondition = msh.boundaryConditions.cbegin(); boundaryCondition != msh.boundaryConditions.cend(); ++boundaryCondition) { + fprintf(file, " BOUNDARY CONDITIONS 2.0.0\n"); + // collect members in group + fprintf(file, "%32d%8d%8d%8d%8d\n", boundaryCondition->first, 1, static_cast(boundaryCondition->second.size()), 0, 6); + for (auto boundary = boundaryCondition->second.begin(); boundary < boundaryCondition->second.end(); ++boundary) { + fprintf(file, "%10d %5d %5d\n", boundary->element+1, GambitInfo::Type, boundary->side+1); + } + fprintf(file, "ENDOFSECTION\n"); + } + + fclose(file); +} + + + +#endif diff --git a/gmsh2gambit/src/GMSH.cpp b/gmsh2gambit/src/GMSH.cpp new file mode 100644 index 0000000..138e283 --- /dev/null +++ b/gmsh2gambit/src/GMSH.cpp @@ -0,0 +1,54 @@ +/** + * @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.) + * + * @section LICENSE + * Copyright (c) 2015, 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 + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * + * @section DESCRIPTION + **/ + +#include "GMSH.h" + +#include +#include + +void error(std::string const& errMessage) +{ + std::cerr << "MSH parse error: " << errMessage << std::endl; + exit(-1); +} + +// Uses GAMBIT neu conventions. See GAMBIT NEUTRAL FILE FORMAT Appendix C.2. +template<> unsigned const Simplex<2>::Face2Nodes[][2] = {{0, 1}, {1, 2}, {2, 0}}; +template<> unsigned const Simplex<3>::Face2Nodes[][3] = {{1, 0, 2}, {0, 1, 3}, {1, 2, 3}, {2, 0, 3}}; + diff --git a/gmsh2gambit/src/GMSH.h b/gmsh2gambit/src/GMSH.h new file mode 100644 index 0000000..f1eb8a0 --- /dev/null +++ b/gmsh2gambit/src/GMSH.h @@ -0,0 +1,211 @@ +/** + * @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.) + * + * @section LICENSE + * Copyright (c) 2015, 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 + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * + * @section DESCRIPTION + **/ +#ifndef GMSH_H_ +#define GMSH_H_ + +#include +#include +#include +#include + +#define MAX_NUMBER_OF_TAGS 16 + +template +struct Simplex { + /* According to MSH mesh format documentation: + * 1 = 2-node line. + * 2 = 3-node triangle. + * 4 = 4-node tetrahedron. + * Will not work for other types (e.g. 11 = 10-node second order tetrahedron). + */ + static unsigned const Type = 1 << (N-1); + static unsigned const NumNodes = N+1; + static unsigned const NumFaces = N+1; + /* Each N-simplex includes N+1 (N-1)-simplices where each has DIM nodes + * E.g. Triangle (N=2): 3 lines with 2 nodes each + * Tetrahedron (N=3): 4 triangles with 3 nodes each */ + static unsigned const Face2Nodes[N+1][N]; + + unsigned nodes[NumNodes]; + unsigned group; +}; + +struct Boundary { + unsigned element; + unsigned side; +}; + +template +struct GMSH { + typedef Simplex Element; + typedef Simplex Face; + + double (*vertices)[3]; + Element* elements; + unsigned numVertices; + unsigned numElements; + std::unordered_map< unsigned, std::vector > materialGroups; + std::unordered_map< unsigned, std::vector > boundaryConditions; + + explicit GMSH(char const* filename); + + ~GMSH() { + delete[] vertices; + delete[] elements; + } +}; + +void error(std::string const& errMessage); + +template +GMSH::GMSH(char const* filename) +{ + numVertices = 0; + numElements = 0; + + Face* faces = NULL; + unsigned numFaces = 0; + + std::ifstream in(filename, std::ifstream::in); + + unsigned tags[MAX_NUMBER_OF_TAGS]; + + // Parse MSH + while (in.good()) { + std::string block; + in >> block; + if (block.compare("$Nodes") == 0) { + in >> numVertices; + vertices = new double[numVertices][3]; + for (unsigned vertex = 0; vertex < numVertices; ++vertex) { + unsigned id; + double x, y, z; + in >> id >> x >> y >> z; + if (id-1 != vertex) { + error("Ids are not contiguous."); + } + vertices[vertex][0] = x; + vertices[vertex][1] = y; + vertices[vertex][2] = z; + } + } else if (block.compare("$Elements") == 0) { + unsigned numGeoShapes; + in >> numGeoShapes; + faces = new Face[numGeoShapes]; + elements = new Element[numGeoShapes]; + + for (unsigned geoShape = 0; geoShape < numGeoShapes; ++geoShape) { + unsigned id, type, numTags; + in >> id >> type >> numTags; + if (type != Face::Type && type != Element::Type) { + error("Unsuported type. Must be 1 (line, 2D), 2 (triangle, 2D and 3D), or 4 (tetrahedron, 3D)."); + } + if (numTags > MAX_NUMBER_OF_TAGS) { + error("Too many tags. Increase MAX_NUMBER_OF_TAGS."); + } + for (unsigned tag = 0; tag < numTags; ++tag) { + in >> tags[tag]; + } + if (type == Face::Type) { + Face& face = faces[numFaces++]; + for (unsigned node = 0; node < Face::NumNodes; ++node) { + in >> face.nodes[node]; + --face.nodes[node]; + } + face.group = tags[0]; + } else if (type == Element::Type) { + Element& element = elements[numElements]; + for (unsigned node = 0; node < Element::NumNodes; ++node) { + in >> element.nodes[node]; + --element.nodes[node]; + } + element.group = tags[0]; + materialGroups[element.group].push_back(numElements); + ++numElements; + } + } + } + } + + in.close(); + + // vertex to triangle map + std::vector* vertex2face = new std::vector[numVertices]; + for (unsigned face = 0; face < numFaces; ++face) { + for (unsigned node = 0; node < Face::NumNodes; ++node) { + vertex2face[faces[face].nodes[node]].push_back(face); + } + } + for (unsigned vertex = 0; vertex < numVertices; ++vertex) { + std::sort(vertex2face[vertex].begin(), vertex2face[vertex].end()); + } + + for (unsigned element = 0; element < numElements; ++element) { + for (unsigned elmtFace = 0; elmtFace < Element::NumFaces; ++elmtFace) { + unsigned nodes[Face::NumNodes]; + for (unsigned node = 0; node < Face::NumNodes; ++node) { + nodes[node] = elements[element].nodes[ Element::Face2Nodes[elmtFace][node] ]; + } + std::vector intersect[Face::NumNodes-1]; + std::set_intersection( vertex2face[ nodes[0] ].begin(), vertex2face[ nodes[0] ].end(), + vertex2face[ nodes[1] ].begin(), vertex2face[ nodes[1] ].end(), + std::back_inserter(intersect[0]) ); + for (unsigned node = 2; node < Face::NumNodes; ++node) { + std::set_intersection( intersect[node-2].begin(), intersect[node-2].end(), + vertex2face[ nodes[node] ].begin(), vertex2face[ nodes[node] ].end(), + std::back_inserter(intersect[node-1]) ); + } + if (!intersect[Face::NumNodes-2].empty()) { + if (intersect[Face::NumNodes-2].size() > 1) { + error("A face of an element exists multiple times in the surface mesh."); + } + Face& face = faces[ intersect[Face::NumNodes-2][0] ]; + Boundary bnd; + bnd.element = element; + bnd.side = elmtFace; + boundaryConditions[face.group].push_back(bnd); + } + } + } + + delete[] vertex2face; + delete[] faces; +} + +#endif diff --git a/gmsh2gambit/src/SConscript b/gmsh2gambit/src/SConscript new file mode 100644 index 0000000..8b985b2 --- /dev/null +++ b/gmsh2gambit/src/SConscript @@ -0,0 +1,49 @@ +#! /usr/bin/python +## +# @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.) +# +# @section LICENSE +# Copyright (c) 2015, 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 +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# + +Import('env') + +env.sourceFiles.extend([ + env.Object('main.cpp'), + env.Object('GMSH.cpp'), + env.Object('GAMBITWriter.cpp') +]) + +Export('env') diff --git a/gmsh2gambit/src/main.cpp b/gmsh2gambit/src/main.cpp new file mode 100644 index 0000000..f79aa4b --- /dev/null +++ b/gmsh2gambit/src/main.cpp @@ -0,0 +1,86 @@ +/** + * @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.) + * + * @section LICENSE + * Copyright (c) 2015, 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 + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * + * @section DESCRIPTION + **/ + +#include "GMSH.h" +#include "GAMBITWriter.h" + +#include +#include + +template +void convert(std::string const& in, std::string const& out) { + std::cout << "Reading MSH..." << std::flush; + GMSH msh(in.c_str()); + std::cout << "finished." << std::endl; + std::cout << "Writing NEU..." << std::flush; + writeNEU(out.c_str(), msh); + std::cout << "finished." << std::endl; +} + +int main(int argc, char** argv) +{ + utils::Args args; + args.addOption("dim", 'd', "Dimension (Defaults to 3, put 2 here for SeisSol 2D)", utils::Args::Required, false); + args.addOption("input", 'i', "Input file (.msh)"); + args.addOption("output", 'o', "Output file (.neu)"); + + if (args.parse(argc, argv) == utils::Args::Success) { + unsigned dim = args.getArgument("dim", 3); + std::string in = args.getArgument("input"); + std::string out = args.getArgument("output"); + + std::cout << "Reading MSH..." << std::flush; + switch (dim) { + case 3: + convert<3>(in, out); + break; + case 2: + std::cerr << "Warning: The z-coordinate will be ignored. If your mesh is not parallel to the z=0 plane, you will be in trouble." << std::endl; + convert<2>(in, out); + break; + default: + std::cerr << "Error: Only dim=2 and dim=3 are supported." << std::endl; + break; + } + } else { + return -1; + } + + return 0; +} diff --git a/gmsh2gambit/submodules b/gmsh2gambit/submodules new file mode 120000 index 0000000..8c556f9 --- /dev/null +++ b/gmsh2gambit/submodules @@ -0,0 +1 @@ +/home/david/Dokumente/uni/phd/software/meshing/Meshing/submodules \ No newline at end of file diff --git a/gmsh_example/planar-anydip.geo b/gmsh_example/planar-anydip.geo new file mode 100644 index 0000000..ea51a1a --- /dev/null +++ b/gmsh_example/planar-anydip.geo @@ -0,0 +1,163 @@ +/* +/** + * @file + * This file is part of SeisSol. + * + * @author Thomas Ulrich + * + * @section LICENSE + * Copyright (c) 2014-2015, 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 + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + +19.02.2016 +Create a planar fault of given dimension and dip angle, the nucleation is also included in the geometry +To use obtain the mesh: +gmsh planar-anydip.geo -3 -optimize +To convert the mesh: +gmsh2gambit -i planar-anydip.msh -o planar-anydip.neu + + */ + + +lc = 25e3; +lc_fault = 300; + +Fault_length = 40e3; +Fault_width = 20e3; +Fault_dip = 60*Pi/180.; + +//Nucleation in X,Z local coordinates +X_nucl = 0e3; +Width_nucl = 0.5*Fault_width; +R_nucl = 2e3; +lc_nucl = 100; + +Xmax = 60e3; +Xmin = -Xmax; +Ymin = -Xmax + 0.5 * Fault_width *Cos(Fault_dip); +Ymax = Xmax + 0.5 * Fault_width *Cos(Fault_dip); +Zmin = -Xmax; + + +//Create the Volume +Point(1) = {Xmin, Ymin, 0, lc}; +Point(2) = {Xmin, Ymax, 0, lc}; +Point(3) = {Xmax, Ymax, 0, lc}; +Point(4) = {Xmax, Ymin, 0, lc}; +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; +Line Loop(5) = {1,2,3,4}; +Plane Surface(1) = {5}; +Extrude {0,0, Zmin} { Surface{1}; } + +//Create the fault +Point(100) = {-0.5*Fault_length, Fault_width *Cos(Fault_dip), -Fault_width *Sin(Fault_dip), lc}; +Point(101) = {-0.5*Fault_length, 0, 0e3, lc}; +Point(102) = {0.5*Fault_length, 0, 0e3, lc}; +Point(103) = {0.5*Fault_length, Fault_width *Cos(Fault_dip), -Fault_width *Sin(Fault_dip), lc}; +Line(100) = {100, 101}; +Line(101) = {101, 102}; +Line{101} In Surface{1}; +Line(102) = {102, 103}; +Line(103) = {103, 100}; + +//create nucleation patch +Point(200) = {X_nucl, Width_nucl*Cos(Fault_dip), -Width_nucl *Sin(Fault_dip), lc_nucl}; +Point(201) = {X_nucl, (Width_nucl + R_nucl) * Cos(Fault_dip), -(Width_nucl+R_nucl) *Sin(Fault_dip), lc_nucl}; +Point(202) = {X_nucl + R_nucl, Width_nucl*Cos(Fault_dip), -Width_nucl *Sin(Fault_dip), lc_nucl}; +Point(203) = {X_nucl, (Width_nucl - R_nucl) * Cos(Fault_dip), -(Width_nucl-R_nucl) *Sin(Fault_dip), lc_nucl}; +Point(204) = {X_nucl - R_nucl, Width_nucl*Cos(Fault_dip), -Width_nucl *Sin(Fault_dip), lc_nucl}; +Circle(200) = {201,200,202}; +Circle(201) = {202,200,203}; +Circle(202) = {203,200,204}; +Circle(203) = {204,200,201}; +Line Loop(204) = {200,201,202,203}; +Plane Surface(200) = {204}; + +Line Loop(105) = {100,101,102,103,200,201,202,203}; +Plane Surface(100) = {105}; + +//There is a bug in "Attractor", we need to define a Ruled surface in FaceList +Line Loop(106) = {100,101,102,103}; +Ruled Surface(101) = {106}; +Ruled Surface(201) = {204}; + +Surface{100,200} In Volume{1}; + + +//1.2 Managing coarsening away from the fault +// Attractor field returns the distance to the curve (actually, the +// distance to 100 equidistant points on the curve) +Field[1] = Attractor; +Field[1].FacesList = {101}; + +// Matheval field returns "distance squared + lc/20" +Field[2] = MathEval; +//Field[2].F = Sprintf("0.02*F1 + 0.00001*F1^2 + %g", lc_fault); +//Field[2].F = Sprintf("0.02*F1 +(F1/2e3)^2 + %g", lc_fault); +Field[2].F = Sprintf("0.05*F1 +(F1/2.5e3)^2 + %g", lc_fault); + +//3.4.5 Managing coarsening around the nucleation Patch +Field[3] = Attractor; +Field[3].FacesList = {201}; + +Field[4] = Threshold; +Field[4].IField = 3; +Field[4].LcMin = lc_nucl; +Field[4].LcMax = lc_fault; +Field[4].DistMin = R_nucl; +Field[4].DistMax = 3*R_nucl; + +Field[5] = Restrict; +Field[5].IField = 4; +Field[5].FacesList = {100,200} ; + +//equivalent of propagation size on element +Field[6] = Threshold; +Field[6].IField = 1; +Field[6].LcMin = lc_fault; +Field[6].LcMax = lc; +Field[6].DistMin = 2*lc_fault; +Field[6].DistMax = 2*lc_fault+0.001; + + +Field[7] = Min; +Field[7].FieldsList = {2,5,6}; + + +Background Field = 7; + +Physical Surface(101) = {1}; +Physical Surface(103) = {100,200}; +//This ones are read in the gui +Physical Surface(105) = {14,18,22,26,27}; + +Physical Volume(1) = {1}; diff --git a/gmsh_example/tpv33_half.geo b/gmsh_example/tpv33_half.geo new file mode 100644 index 0000000..402df3c --- /dev/null +++ b/gmsh_example/tpv33_half.geo @@ -0,0 +1,139 @@ +lc= 5000; +lc_fault = 100; + +Point(1) = {-50000, 0, -50000, lc}; +Point(2) = {-50000, 0, 0, lc}; +Point(3) = {50000, 0, 0, lc}; +Point(4) = {50000, 0, -50000, lc}; +Point(100) = {-12000, 0, -10000, lc}; +Point(101) = {4000, 0, -10000, lc}; +Point(102) = {4000, 0, 0, lc}; +Point(103) = {-12000, 0, 0, lc}; +Point(200) = {-6000, 0, -6000, lc_fault}; +Point(201) = {-6000, 0, -5450, lc_fault}; +Point(202) = {-5450, 0, -6000, lc_fault}; +Point(203) = {-6000, 0, -6550, lc_fault}; +Point(204) = {-6550, 0, -6000, lc_fault}; +Point(211) = {-6000, 0, -5200, lc_fault}; +Point(212) = {-5200, 0, -6000, lc_fault}; +Point(213) = {-6000, 0, -6800, lc_fault}; +Point(214) = {-6800, 0, -6000, lc_fault}; +Point(1001) = {-50000, 800, -50000, lc}; +Point(1002) = {-50000, 800, 0, lc}; +Point(1003) = {50000, 800, 0, lc}; +Point(1004) = {50000, 800, -50000, lc}; +Point(1011) = {-50000, 50000, -50000, lc}; +Point(1012) = {-50000, 50000, 0, lc}; +Point(1013) = {50000, 50000, 0, lc}; +Point(1014) = {50000, 50000, -50000, lc}; +Point(1100) = {-12000, 800, -10000, lc}; +Point(1101) = {4000, 800, -10000, lc}; +Point(1102) = {4000, 800, 0, lc}; +Point(1103) = {-12000, 800, 0, lc}; +Line(1) = {1, 2}; +Line(2) = {2, 103}; +Line(3) = {103, 100}; +Line(4) = {100, 101}; +Line(5) = {101, 102}; +Line(6) = {102, 3}; +Line(7) = {3, 4}; +Line(8) = {4, 1}; +Line(9) = {103, 102}; +Circle(200) = {201, 200, 202}; +Circle(201) = {202, 200, 203}; +Circle(202) = {203, 200, 204}; +Circle(203) = {204, 200, 201}; +Circle(210) = {211, 200, 212}; +Circle(211) = {212, 200, 213}; +Circle(212) = {213, 200, 214}; +Circle(213) = {214, 200, 211}; +Line(216) = {2, 1002}; +Line(217) = {1002, 1012}; +Line(218) = {3, 1003}; +Line(219) = {1003, 1013}; +Line(220) = {4, 1004}; +Line(221) = {1004, 1014}; +Line(222) = {1014, 1013}; +Line(223) = {1004, 1003}; +Line(224) = {1003, 1102}; +Line(225) = {1103, 1102}; +Line(226) = {1102, 1101}; +Line(227) = {1101, 1100}; +Line(228) = {1100, 1103}; +Line(229) = {101, 1101}; +Line(230) = {102, 1102}; +Line(231) = {103, 1103}; +Line(232) = {1103, 1002}; +Line(233) = {1001, 1}; +Line(234) = {1004, 1001}; +Line(235) = {1011, 1012}; +Line(236) = {1013, 1012}; +Line(237) = {1011, 1014}; +Line(238) = {1011, 1001}; +Line(239) = {100, 1100}; +Line(240) = {1001, 1002}; +Line Loop(1) = {1, 2, 3, 4, 5, 6, 7, 8}; +Plane Surface(1) = {1}; +Line Loop(2) = {3, 4, 5, -9, 210, 211, 212, 213}; +Plane Surface(2) = {2}; +Line Loop(3) = {210, 211, 212, 213, 200, 201, 202, 203}; +Plane Surface(3) = {3}; +Line Loop(4) = {200, 201, 202, 203}; +Plane Surface(4) = {4}; +Line Loop(242) = {1, 216, -240, 233}; +Plane Surface(242) = {242}; +Line Loop(244) = {238, 240, 217, -235}; +Plane Surface(244) = {244}; +Line Loop(246) = {238, -234, 221, -237}; +Plane Surface(246) = {246}; +Line Loop(248) = {234, 233, -8, 220}; +Plane Surface(248) = {248}; +Line Loop(250) = {237, 222, 236, -235}; +Plane Surface(250) = {250}; +Line Loop(252) = {217, -236, -219, 224, -225, 232}; +Plane Surface(252) = {252}; +Line Loop(254) = {221, 222, -219, -223}; +Plane Surface(254) = {254}; +Line Loop(256) = {220, 223, -218, 7}; +Plane Surface(256) = {256}; +Line Loop(258) = {6, 218, 224, -230}; +Plane Surface(258) = {258}; +Line Loop(260) = {9, 230, -225, -231}; +Plane Surface(260) = {260}; +Line Loop(262) = {216, -232, -231, -2}; +Plane Surface(262) = {262}; +Line Loop(264) = {228, -231, 3, 239}; +Plane Surface(264) = {264}; +Line Loop(266) = {227, -239, 4, 229}; +Plane Surface(266) = {266}; +Line Loop(268) = {229, -226, -230, -5}; +Plane Surface(268) = {268}; +Line Loop(270) = {234, 240, -232, -228, -227, -226, -224, -223}; +Plane Surface(270) = {270}; +Line Loop(272) = {228, 225, 226, 227}; +Plane Surface(272) = {272}; +Line Loop(10000) = {3, 4, 5, -9}; +Ruled Surface(10000) = {10000}; +Surface Loop(274) = {244, 246, 254, 250, 252, 270, 272}; +Volume(274) = {274}; +Surface Loop(276) = {272, 268, 266, 264, 2, 260, 3, 4}; +Volume(276) = {276}; +Surface Loop(278) = {1, 242, 262, 248, 256, 258, 270, 268, 266, 264}; +Volume(278) = {278}; +Field[1] = Attractor; +Field[1].FacesList = {10000}; +Field[1].FieldX = -1; +Field[1].FieldY = -1; +Field[1].FieldZ = -1; +Field[1].NNodesByEdge = 20; +Field[2] = MathEval; +//Field[2].F = Sprintf("0.05*F1 +(F1/2.5e3)^2 + %g", lc_fault); +Field[2].F = Sprintf("0.1*F1 +(F1/5.0e3)^2 + %g", lc_fault); +Background Field = 2; + +Physical Surface(101) = {252, 258, 260, 262}; +Physical Surface(105) = {242, 244, 246, 248, 250, 254, 256}; +Physical Surface(103) = {2, 3, 4}; +Physical Volume(2) = {276,278}; +Physical Volume(3) = {274}; + diff --git a/legacy_scripts/File_ts2STL.m b/legacy_scripts/File_ts2STL.m new file mode 100644 index 0000000..41b7ec7 --- /dev/null +++ b/legacy_scripts/File_ts2STL.m @@ -0,0 +1,142 @@ +function [Vertex,Conectivity]=File_ts2STL(name) +%% +% @file +% This file is part of SeisSol. +% +% @author Christian Pelties (pelties AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/pelties) +% +% @section LICENSE +% Copyright (c) 2010, 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 +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. + +disp(' '),disp(' ') +disp('File_ts2STL') +disp(' '),disp(' ') +disp('converts GoCads ts files in STL files read by ICEM') + +% Open file to read +hdr = strcat(name,'.ts'); +fid = fopen(hdr); +% Read the full file character by character and store the ascii code in one +% single array +Ascii = fread(fid); +% Close file +fclose(fid); + + + +nascii = length(Ascii); +icolumn =1; irow = 1; +buildnumber = 0; assignnumber = 0; number = 0; power = 0; +sc = 0; neg=1; +for iascii = 1:nascii + data = Ascii(iascii); + if data==10 %new line + if sc==0 + irow = irow+1; + icolumn = 1; + else + assignnumber=1; + end + elseif data==32 %space + assignnumber = 1; + elseif data==45 %- + neg = -1; + elseif data==46 %. + sign = -1; + elseif data>=48 && data<=57 %decimal numbers + buildnumber = 1; + elseif data==76 && sc == 0 %L + sc = 1; + icolumn =1; irow = 1; + buildnumber = 0; assignnumber = 0; number = 0; power = 0; + end + + if buildnumber==1 + digit = neg*str2num(char(data)); + if sign==1 + number = number*10 + digit; + elseif sign==-1 + power = power+1; + number = number + digit/(10^power); + end + buildnumber = 0; + end + if assignnumber==1 + if sc==0 + Vertex(icolumn,irow) = number; + else + Conectivity(icolumn,irow) = number; + end + icolumn = icolumn+1; + sign = 1; + power = 0; + number = 0; + assignnumber = 0; + neg = 1; + if data==10 && sc==1 + irow = irow+1; + icolumn = 1; + end + end + +end + +file = strcat(name,'.stl'); +fid=fopen(file,'w'); +fprintf(fid,'solid '); +fprintf(fid,name); +fprintf(fid,'\n'); + +for ielem=1:length(Conectivity) + x=Conectivity(2,ielem); + y=Conectivity(3,ielem); + z=Conectivity(4,ielem); + fprintf(fid,'facet normal 0 0 0\n'); + fprintf(fid,'outer loop\n'); + fprintf(fid,'vertex%20.11e%20.11e%20.11e\n',Vertex(3,x),Vertex(4,x),Vertex(5,x)); + fprintf(fid,'vertex%20.11e%20.11e%20.11e\n',Vertex(3,y),Vertex(4,y),Vertex(5,y)); + fprintf(fid,'vertex%20.11e%20.11e%20.11e\n',Vertex(3,z),Vertex(4,z),Vertex(5,z)); + fprintf(fid,'endloop\n'); + fprintf(fid,'endfacet\n'); +end +fprintf(fid,'endsolid '); +fprintf(fid,name); +fclose(fid) +return + + + + + + + + + + diff --git a/legacy_scripts/Hex_Refinement.f90 b/legacy_scripts/Hex_Refinement.f90 new file mode 100644 index 0000000..05566e1 --- /dev/null +++ b/legacy_scripts/Hex_Refinement.f90 @@ -0,0 +1,483 @@ +!> +!! @file +!! This file is part of SeisSol. +!! +!! @section LICENSE +!! Copyright (c) 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 +!! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +!! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!! POSSIBILITY OF SUCH DAMAGE. +PROGRAM Hex_refinement + !-------------------------------------------------------------------------- + IMPLICIT NONE + !-------------------------------------------------------------------------- + ! Type declaration + INTEGER :: uid1, uid2 + INTEGER :: c, ix, iy, iz, nline, nrest, ng1, ng2 + INTEGER :: clocal(32), ind(8) + INTEGER :: nVert, nElem, nl, nv, ntr + INTEGER :: nx, nx2, ny, ny2, nz, nz2 + REAL :: xmin, xmax, ymin, ymax, zmin, zmax, zb + REAL :: dx, dy, dz + REAL, ALLOCATABLE :: X(:,:) + INTEGER, ALLOCATABLE :: HEX(:,:) + INTEGER, ALLOCATABLE :: member(:) + CHARACTER (LEN=350) :: InFilename + CHARACTER (LEN=350) :: OutFilename + CHARACTER (LEN=70) :: StringLine + !-------------------------------------------------------------------------- + ! + WRITE(*,*) ' | -------------------------------------------------- |' + WRITE(*,*) ' | HEX_REFINEMENT |' + WRITE(*,*) ' | -------------------------------------------------- |' + ! + WRITE(*,*) ' ' + WRITE(*,*) ' Enter the block (brick) dimensions! ' + WRITE(*,*) ' ' + WRITE(*,*) ' x_min x_max ' + READ (*,*) xmin, xmax + WRITE(*,*) ' y_min y_max ' + READ (*,*) ymin, ymax + WRITE(*,*) ' z_min z_max ' + READ (*,*) zmin, zmax + WRITE(*,*) ' Enter depth of material boundary under which mesh is coarsened! ' + WRITE(*,*) ' ' + READ (*,*) zb + WRITE(*,*) ' Enter numbers of intervals in x,y,z-directions for the coarse mesh! ' + WRITE(*,*) ' ' + READ (*,*) nx2, ny2, nz2 + WRITE(*,*) ' Enter numbers of intervals in z-directions for the fine mesh! ' + WRITE(*,*) ' ' + READ (*,*) nz + + nx = nx2*3 + ny = ny2*3 + + dx = (xmax-xmin)/nx + dy = (ymax-ymin)/ny + dz = (zmax-zb) /nz + + !OutFilename = TRIM(InFilename)//'.met' + !InFilename = TRIM(InFilename)//'.neu' + WRITE(*,*) ' Give output file name ' + READ (*,*) OutFilename + + !nVert is number of vertices, i.e. of fine mesh, course mesh and transition layer + nVert = (nx+1)*(ny+1)*(nz+1) + (nx2+1)*(ny2+1)*(nz2) + (nx2*ny2*8+nx2*2+ny2*2) + !nVert is number of elements, i.e. of fine mesh, course mesh and transition layer + nElem = (nx*ny*nz) + (nx2*ny2*(nz2-1)) + (nx2*ny2*13) + + + ALLOCATE( X(nVert,3) ) + ALLOCATE( HEX(nElem,8) ) + + !-------------------------------------------------------------------------- + ! Compute vertices of fine mesh + + c = 0 + + DO iz = 1, nz+1 + DO iy = 1, ny+1 + DO ix = 1, nx+1 + + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + X(c,2) = ymin + (iy-1)*dy + X(c,3) = zmax - (iz-1)*dz + + ENDDO + ENDDO + ENDDO + + !-------------------------------------------------------------------------- + ! Compute vertices of transition mesh + + dx = 3*dx + dy = 3*dy + dz = (zb-zmin)/nz2 + + DO iy = 1, ny2 + DO ix = 1, nx2 + + ! four inner vertices + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + dx/3 + X(c,2) = ymin + (iy-1)*dy + dy/3 + X(c,3) = zb - dz/3 + + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + dx*2/3 + X(c,2) = ymin + (iy-1)*dy + dy /3 + X(c,3) = zb - dz /3 + + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + dx*2/3 + X(c,2) = ymin + (iy-1)*dy + dy*2/3 + X(c,3) = zb - dz /3 + + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + dx/3 + X(c,2) = ymin + (iy-1)*dy + dy*2/3 + X(c,3) = zb - dz /3 + + ! four boundary vertices + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + X(c,2) = ymin + (iy-1)*dy + dy*2/3 + X(c,3) = zb - dz/2 + + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + X(c,2) = ymin + (iy-1)*dy + dy/3 + X(c,3) = zb - dz/2 + + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + dx/3 + X(c,2) = ymin + (iy-1)*dy + X(c,3) = zb - dz*2/3 + + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + dx*2/3 + X(c,2) = ymin + (iy-1)*dy + X(c,3) = zb - dz*2/3 + + ENDDO + + ! two boundary vertices at the end of a row in x-direction + c = c+1 + + X(c,1) = xmax + X(c,2) = ymin + (iy-1)*dy + dy*2/3 + X(c,3) = zb - dz/2 + + c = c+1 + + X(c,1) = xmax + X(c,2) = ymin + (iy-1)*dy + dy/3 + X(c,3) = zb - dz/2 + + ENDDO + + DO ix = 1, nx2 + ! two boundary vertices at the end of a row in y-direction + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + dx/3 + X(c,2) = ymax + X(c,3) = zb - dz*2/3 + + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + dx*2/3 + X(c,2) = ymax + X(c,3) = zb - dz*2/3 + + ENDDO + + !-------------------------------------------------------------------------- + ! Compute vertices of coarse mesh + + DO iz = 1, nz2 + DO iy = 1, ny2+1 + DO ix = 1, nx2+1 + + c = c+1 + + X(c,1) = xmin + (ix-1)*dx + X(c,2) = ymin + (iy-1)*dy + X(c,3) = zb - iz*dz + + ENDDO + ENDDO + ENDDO + + WRITE(*,*) ' | -------------------------------------------------- ' + WRITE(*,*) ' |',c,' vertices computed ! ' + WRITE(*,*) ' | -------------------------------------------------- ' + WRITE(*,*) ' ' + + !-------------------------------------------------------------------------- + ! Build connectivity of fine mesh + + c = 0 + nl = (nx+1)*(ny+1) + + DO iz = 1, nz + DO iy = 1, ny + DO ix = 1, nx + + c = c+1 + + HEX(c,1) = iz*nl + (iy-1)*(nx+1) + ix + HEX(c,2) = iz*nl + (iy-1)*(nx+1) + ix + 1 + HEX(c,3) = iz*nl + iy*(nx+1) + ix + HEX(c,4) = iz*nl + iy*(nx+1) + ix + 1 + HEX(c,5) = (iz-1)*nl + (iy-1)*(nx+1) + ix + HEX(c,6) = (iz-1)*nl + (iy-1)*(nx+1) + ix + 1 + HEX(c,7) = (iz-1)*nl + iy*(nx+1) + ix + HEX(c,8) = (iz-1)*nl + iy*(nx+1) + ix + 1 + + ENDDO + ENDDO + ENDDO + + !-------------------------------------------------------------------------- + ! Build connectivity of transition mesh + + nv = nl * nz + ntr = nv + nl + nx2*ny2*8 + (nx2+ny2)*2 + + DO iy = 1, ny2 + DO ix = 1, nx2 + + ! First determine absolute vertex numbers for the + ! local 32 vectices of a refined hexahedron + + clocal(1) = nv + ((iy-1)*(nx+1) + (ix-1))*3 + 1 + clocal(2) = clocal(1) + 1 + clocal(3) = clocal(2) + 1 + clocal(4) = clocal(3) + 1 + clocal(5) = nv + ((iy-1)*(nx+1) + (ix-1))*3 + nx+1 + 1 + clocal(6) = clocal(5) + 1 + clocal(7) = clocal(6) + 1 + clocal(8) = clocal(7) + 1 + clocal(9) = nv + ((iy-1)*(nx+1) + (ix-1))*3 + (nx+1)*2 + 1 + clocal(10) = clocal(9) + 1 + clocal(11) = clocal(10) + 1 + clocal(12) = clocal(11) + 1 + clocal(13) = nv + ((iy-1)*(nx+1) + (ix-1))*3 + (nx+1)*3 + 1 + clocal(14) = clocal(13) + 1 + clocal(15) = clocal(14) + 1 + clocal(16) = clocal(15) + 1 + clocal(17) = nv + nl + (iy-1)*(8*nx2+2) + (ix-1)*8 + 1 + clocal(18) = clocal(17) + 1 + clocal(19) = clocal(18) + 1 + clocal(20) = clocal(19) + 1 + clocal(21) = clocal(20) + 1 + clocal(22) = clocal(21) + 1 + clocal(23) = clocal(22) + 1 + clocal(24) = clocal(23) + 1 + ! Special treatment of boundary vertices at the end of each row in x-direction + IF(ix.LT.nx2)THEN + clocal(25) = nv + nl + (iy-1)*(8*nx2+2) + ix*8 + 5 + ELSE + clocal(25) = nv + nl + (iy-1)*(8*nx2+2) + ix*8 + 1 + END IF + clocal(26) = clocal(25) + 1 + ! Special treatment of boundary vertices at the end of each row in y-direction + IF(iy.LT.ny2)THEN + clocal(27) = nv + nl + iy*(8*nx2+2) + (ix-1)*8 + 7 + ELSE + clocal(27) = nv + nl + iy*(8*nx2+2) + (ix-1)*2 + 1 + END IF + clocal(28) = clocal(27) + 1 + clocal(29) = ntr + (iy-1)*(nx2+1) + ix + clocal(30) = clocal(29) + 1 + clocal(31) = ntr + iy*(nx2+1) + ix + clocal(32) = clocal(31) + 1 + + ! Now build the connectivity of the 13 sub-hexahedrons with the local indices "ind" + ! Local hexahedron 1 + c = c+1 + ind = (/29,23,22,17,1,2,5,6/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 2 + c = c+1 + ind = (/23,24,17,18,2,3,6,7/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 3 + c = c+1 + ind = (/24,30,18,26,3,4,7,8/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 4 + c = c+1 + ind = (/22,17,21,20,5,6,9,10/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 5 + c = c+1 + ind = (/17,18,20,19,6,7,10,11/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 6 + c = c+1 + ind = (/18,26,19,25,7,8,11,12/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 7 + c = c+1 + ind = (/21,20,31,27,9,10,13,14/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 8 + c = c+1 + ind = (/20,19,27,28,10,11,14,15/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 9 + c = c+1 + ind = (/19,25,28,32,11,12,15,16/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 10 + c = c+1 + ind = (/23,24,27,28,17,18,20,19/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 11 + c = c+1 + ind = (/29,23,31,27,22,17,21,20/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 12 + c = c+1 + ind = (/29,30,31,32,23,24,27,28/) + HEX(c,1:8) = clocal(ind) + ! Local hexahedron 13 + c = c+1 + ind = (/24,30,28,32,18,26,19,25/) + HEX(c,1:8) = clocal(ind) + + ENDDO + ENDDO + + !-------------------------------------------------------------------------- + ! Build connectivity of coarse mesh + + nl = (nx2+1)*(ny2+1) + + DO iz = 1, nz2-1 + DO iy = 1, ny2 + DO ix = 1, nx2 + + c = c+1 + + HEX(c,1) = ntr + iz*nl + (iy-1)*(nx2+1) + ix + HEX(c,2) = ntr + iz*nl + (iy-1)*(nx2+1) + ix + 1 + HEX(c,3) = ntr + iz*nl + iy*(nx2+1) + ix + HEX(c,4) = ntr + iz*nl + iy*(nx2+1) + ix + 1 + HEX(c,5) = ntr + (iz-1)*nl + (iy-1)*(nx2+1) + ix + HEX(c,6) = ntr + (iz-1)*nl + (iy-1)*(nx2+1) + ix + 1 + HEX(c,7) = ntr + (iz-1)*nl + iy*(nx2+1) + ix + HEX(c,8) = ntr + (iz-1)*nl + iy*(nx2+1) + ix + 1 + + ENDDO + ENDDO + ENDDO + + WRITE(*,*) ' | -------------------------------------------------- ' + WRITE(*,*) ' |',c,' hexahedrons built ! ' + WRITE(*,*) ' | -------------------------------------------------- ' + WRITE(*,*) ' ' + + + !-------------------------------------------------------------------------- + ! Writing output + WRITE(*,*) ' | -------------------------------------------------- ' + WRITE(*,*) ' | Writing output to ',TRIM(OutFilename) + WRITE(*,*) ' | -------------------------------------------------- ' + WRITE(*,*) ' ' + uid2 = 11 + + + OPEN( UNIT = uid2, FILE = TRIM(OutFilename), STATUS = 'UNKNOWN' ) + ! Write GAMBIT header + WRITE(uid2,'(" CONTROL INFO 2.2.30")') + WRITE(uid2,'("** GAMBIT NEUTRAL FILE")') + WRITE(uid2,'("default_id3556")') + WRITE(uid2,'("PROGRAM: Gambit VERSION: 2.2.30")') + WRITE(uid2,'(" Feb 2007")') + WRITE(uid2,'(" NUMNP NELEM NGRPS NBSETS NDFCD NDFVL")') + WRITE(uid2,'(I10,I10,I10,I10,I10,I10)') nVert, nElem, 2,0,3,3 + WRITE(uid2,'("ENDOFSECTION")') + WRITE(uid2,'(" NODAL COORDINATES 2.2.30")') + ! Write vertices + DO c = 1, nVert + WRITE(uid2,'(I10, E20.11E3, E20.11E3, E20.11E3)') c, X(c,:) + ENDDO + + WRITE(uid2,'("ENDOFSECTION")') + WRITE(uid2,'(" ELEMENTS/CELLS 2.2.30")') + ! Write connectivity + DO c = 1,nElem + WRITE(uid2,'(I8, I3, I3, I9,I8,I8,I8,I8,I8,I8)') c, 4, 8, HEX(c,1:7) + WRITE(uid2,'(I23)') HEX(c,8) + ENDDO + WRITE(uid2,'("ENDOFSECTION")') + + ! Write refined layer as group 1 + ng1 = nx*ny*nz + WRITE(uid2,'(" ELEMENT GROUP 2.2.30")') + WRITE(uid2,'("GROUP: 1 ELEMENTS:",I11," MATERIAL:",I11," NFLAGS:",I11)') ng1,4,1 + WRITE(uid2,'(" layer")') + WRITE(uid2,'(I8)') 0 + + ALLOCATE( member(ng1) ) + + DO c = 1,ng1 + member(c) = c + ENDDO + nline = FLOOR(ng1/10.) + 1 + nrest = MOD(ng1,10) + DO c = 1, nline-1 + WRITE(uid2,'(I8,I8,I8,I8,I8,I8,I8,I8,I8,I8)') member((c-1)*10+1:c*10) + ENDDO + IF(nrest.GE.1)THEN + WRITE(uid2,'(I8,I8,I8,I8,I8,I8,I8,I8,I8,I8)') member(ng1-(nrest-1):ng1) + END IF + ! Write coarse halfspace as group 2 + ng2 = nx2*ny2*((nz2-1)+13) + WRITE(uid2,'("ENDOFSECTION")') + WRITE(uid2,'(" ELEMENT GROUP 2.2.30")') + WRITE(uid2,'("GROUP: 2 ELEMENTS:",I11," MATERIAL:",I11," NFLAGS:",I11)') ng2,4,1 + WRITE(uid2,'(" halfspace")') + WRITE(uid2,'(I8)') 0 + DEALLOCATE( member ) + ALLOCATE( member(ng2) ) + DO c = 1,ng2 + member(c) = ng1 + c + ENDDO + nline = FLOOR(ng2/10.) + 1 + nrest = MOD(ng2,10) + DO c = 1, nline-1 + WRITE(uid2,'(I8,I8,I8,I8,I8,I8,I8,I8,I8,I8)') member((c-1)*10+1:c*10) + ENDDO + IF(nrest.GE.1)THEN + WRITE(uid2,'(I8,I8,I8,I8,I8,I8,I8,I8,I8,I8)') member(ng2-(nrest-1):ng2) + END IF + WRITE(uid2,'("ENDOFSECTION")') + + CLOSE(uid2) + + WRITE(*,*) ' | -------------------------------------------------- |' + WRITE(*,*) ' | HEX_REFINEMENT finished successfully |' + WRITE(*,*) ' | -------------------------------------------------- |' + WRITE(*,*) ' ' + + +END PROGRAM Hex_refinement diff --git a/legacy_scripts/README_hex2tet.txt b/legacy_scripts/README_hex2tet.txt new file mode 100644 index 0000000..208742e --- /dev/null +++ b/legacy_scripts/README_hex2tet.txt @@ -0,0 +1,77 @@ +%% +% @file +% This file is part of SeisSol. +% +% @author Martin Kaeser (martin.kaeser AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/kaeser) +% +% @section LICENSE +% Copyright (c) 2005, 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 +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %% %% + %% GAMBIT_HEX2TET %% + %% Creates regular tetrahedral mesh %% + %% on the base of an hexahedral mesh %% + %% %% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + Description of usage: + + GAMBIT_HEX2TET works only for domains that are meshed with a regular + hexahedral mesh, i.e. the numbers nx, ny, and nz of hexahedral elements + in each space dimension is constant. + + Steps: + + (1) Create a hexagonal mesh with GAMBIT + + (2) Export the mesh as meshfile.neu + + (3) Run GAMBIT_HEX2TET + + - specify the meshfile with the relativ path name + with respect to the path, where GAMBIT_HEX2TET.m is stored + The suffix .neu is assumed as standard, so just use e.g. + ../directory1/directory2/meshfile + + - specify the numbers nx, ny, and nz of hexahedral elements + in x, y, and z dimension + + - this automatically produces a regular tetrahedral mesh + where all hexahedral elements are subdivided into 5 + tetrahedrons and the new tetrahedral mesh is stored in + GAMBIT format in the same directory as meshfile.neu + with the new name meshfile_tetra.neu + + (4) Import the new meshfile_tetra.neu in GAMBIT for further + processing, e.g. boundary conditions, ... + + diff --git a/legacy_scripts/XYinElement.m b/legacy_scripts/XYinElement.m new file mode 100644 index 0000000..6a6a691 --- /dev/null +++ b/legacy_scripts/XYinElement.m @@ -0,0 +1,61 @@ +function inside = XYinElement(x,y,xS,yS) +%% +% @file +% This file is part of SeisSol. +% +% @author Christian Pelties (pelties AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/pelties) +% +% @section LICENSE +% Copyright (c) 2014, 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 +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. +% +% @section DESCRIPTION +% XYinElement +% XYInElement returns 1 if the point (xS/yS) is inside the +% element defined by XY, else it returns 0 +% adapted from SeisSol2Ds XYinElement function + +epsilon = 0.00001; % tolerance + +refFactor = 0.5; +VOL = volume(x,y); + +xi = refFactor/VOL*( ( x(3)*y(1) - x(1)*y(3) ) + xS*(y(3)-y(1)) + yS*(x(1)-x(3)) ); +eta = refFactor/VOL*( ( x(1)*y(2) - x(2)*y(1) ) + xS*(y(1)-y(2)) + yS*(x(2)-x(1)) ); + + +if (xi <(0.0-epsilon)) || (eta <(0.0-epsilon)) || (eta > (1.0-xi+epsilon)) + inside = 0; +else + inside = 1; +end + +function VOL = volume(x,y) +% Sarrus +VOL = 0.5*( (x(2)-x(1))*(y(3)-y(1))-(x(3)-x(1))*(y(2)-y(1)) ); diff --git a/legacy_scripts/gambit_check_insphere.m b/legacy_scripts/gambit_check_insphere.m new file mode 100644 index 0000000..3c81b46 --- /dev/null +++ b/legacy_scripts/gambit_check_insphere.m @@ -0,0 +1,138 @@ +%% +% @file +% This file is part of SeisSol. +% +% @author Martin Kaeser (martin.kaeser AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/kaeser) +% +% @section LICENSE +% Copyright (c) 2005, 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 +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. + +home; +disp(' ') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' %% %%') +disp(' %% GAMBIT_CHECK_INSPHER %%') +disp(' %% TO ANALYSE QUALITY OF 3D_MESH %%') +disp(' %% %%') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' ') + +clear, close all; + +filename = input(' Filename(suffix ".neu" assumed): ','s'); +disp(' ') +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\t Reading data from: \t\t%s',[filename,'.neu'])); +disp(' ') +fid = fopen([filename ,'.neu']); +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +param = fscanf(fid,'%i',[6,1]); NX = param(1); NT = param(2); +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +X = fscanf(fid,'%g',[4,NX]); X = X'; X(:,1) = []; +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +tetra = fscanf(fid,'%g',[7,NT]); tetra = tetra'; tetra(:,1:3) = []; +junk = fgetl(fid); +fclose(fid); +disp(sprintf('\t Read \t\t\t\t\t%i vertices',NX)); +disp(sprintf('\t Read \t\t\t\t\t%i tetrahedral elements',NT)); + +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\n\t Start analysis ...')); + +t=cputime; +%compute edges +a = (X(tetra(:,2),:)-X(tetra(:,1),:)); +b = (X(tetra(:,3),:)-X(tetra(:,1),:)); +c = (X(tetra(:,4),:)-X(tetra(:,1),:)); +d = (X(tetra(:,3),:)-X(tetra(:,2),:)); +e = (X(tetra(:,4),:)-X(tetra(:,2),:)); +%compute volumes +V = abs(sum(a.*cross(b,c),2)/6); +%compute face areas +s = cross(a,b); +A(:,1) = 0.5*sqrt(s(:,1).^2+s(:,2).^2+s(:,3).^2); +s = cross(a,c); +A(:,2) = 0.5*sqrt(s(:,1).^2+s(:,2).^2+s(:,3).^2); +s = cross(b,c); +A(:,3) = 0.5*sqrt(s(:,1).^2+s(:,2).^2+s(:,3).^2); +s = cross(d,e); +A(:,4) = 0.5*sqrt(s(:,1).^2+s(:,2).^2+s(:,3).^2); +%compute total surfaces area +S = sum(A,2); +%compute in-radius +r_in = 3*V./S; + +tt = cputime-t; + +figure(1); plot(r_in); xlabel('Number of Tetrahedron'); ylabel('Insphere Radius') + +[r_min,i] = min(r_in); +disp(sprintf('\n\t Smallest insphere radius : %g ',r_min)); +disp(sprintf('\n\t Smallest insphere found at tetrahedron : %i ',i)); +disp(sprintf('\n\t with vertices : ')); +disp([X(tetra(i,:),:)]); +x_min = min(X(:,1)); x_max = max(X(:,1)); +y_min = min(X(:,2)); y_max = max(X(:,2)); +z_min = min(X(:,3)); z_max = max(X(:,3)); +%plotting tetrahedron with smallest insphere +i = i(1); +s_vert(1,:) = [1,3,2]; s_vert(2,:) = [1,2,4]; s_vert(3,:) = [1,4,3]; s_vert(4,:) = [2,3,4]; +tmp = [X(tetra(i,:),1),X(tetra(i,:),2),X(tetra(i,:),3)]; +tx_min = min(tmp(:,1)); tx_max = max(tmp(:,1)); +ty_min = min(tmp(:,2)); ty_max = max(tmp(:,2)); +tz_min = min(tmp(:,3)); tz_max = max(tmp(:,3)); +tMAX = max(abs([tx_min,tx_max,ty_min,ty_max,tz_min,tz_max])); +figure(2) +subplot(121), hold on +for j = 1:4 + fill3(tmp(s_vert(j,:),1),tmp(s_vert(j,:),2),tmp(s_vert(j,:),3),'r'); +end +axis equal, grid on, axis([x_min x_max y_min y_max z_min z_max]); +xlabel('Location in x');ylabel('Location in y');zlabel('Location in z'); +subplot(122), hold on +for j = 1:4 + fill3(tmp(s_vert(j,:),1),tmp(s_vert(j,:),2),tmp(s_vert(j,:),3),'r'); +end +axis equal, grid on +axis([tx_min-0.1*tMAX tx_max+0.1*tMAX ty_min-0.1*tMAX ty_max+0.1*tMAX tz_min-0.1*tMAX tz_max+0.1*tMAX]); +xlabel('zoomed in x');ylabel('zoomed in y');zlabel('zoomed in z'); +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\n\t Analysis finished successfully! (%g CPU sec)\n',tt)); +disp('-----------------------------------------------------------------------------------') +disp(' ') diff --git a/legacy_scripts/gambit_check_stencil.m b/legacy_scripts/gambit_check_stencil.m new file mode 100644 index 0000000..2cadf7d --- /dev/null +++ b/legacy_scripts/gambit_check_stencil.m @@ -0,0 +1,121 @@ +%% +% @file +% This file is part of SeisSol. +% +% @author Martin Kaeser (martin.kaeser AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/kaeser) +% +% @section LICENSE +% Copyright (c) 2005, 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 +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. + +home; +disp(' ') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' %% %%') +disp(' %% GAMBIT_CHECK_INSPHER %%') +disp(' %% TO ANALYSE QUALITY OF 3D_MESH %%') +disp(' %% %%') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' ') + +close all; + +filename = input(' Filename(suffix ".neu" assumed): ','s'); + +disp(' ') +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\t Reading data from: \t\t%s',[filename,'.neu'])); +disp(' ') +fid = fopen([filename ,'.neu']); +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +param = fscanf(fid,'%i',[6,1]); NX = param(1); NT = param(2); +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +X = fscanf(fid,'%g',[4,NX]); X = X'; X(:,1) = []; +junk = fgetl(fid); +junk = fgetl(fid); +junk = fgetl(fid); +tetra = fscanf(fid,'%g',[7,NT]); tetra = tetra'; tetra(:,1:3) = []; +junk = fgetl(fid); +fclose(fid); +disp(sprintf('\t Read \t\t\t\t\t%i vertices',NX)); +disp(sprintf('\t Read \t\t\t\t\t%i tetrahedral elements',NT)); + +disp('-----------------------------------------------------------------------------------') +eval(['load ',filename,'_stencil.dat']); +stencil = cube_005_stencil; +n = input( 'Given number of tetrahedron: '); +stencil = stencil(n,:); +figure(1) +s_vert(1,:) = [1,3,2]; s_vert(2,:) = [1,2,4]; s_vert(3,:) = [1,4,3]; s_vert(4,:) = [2,3,4]; +for i = 1:length(stencil) + t = stencil(i); + tmp = [X(tetra(t,:),1),X(tetra(t,:),2),X(tetra(t,:),3)]; + for j = 1:4 + if i==1 + fill3(tmp(s_vert(j,:),1),tmp(s_vert(j,:),2),tmp(s_vert(j,:),3),'b'); hold on, grid on + else + fill3(tmp(s_vert(j,:),1),tmp(s_vert(j,:),2),tmp(s_vert(j,:),3),'r'); hold on + end + end + alpha(0.5), axis equal + pause +end + +x_min = min(X(:,1)); x_max = max(X(:,1)); +y_min = min(X(:,2)); y_max = max(X(:,2)); +z_min = min(X(:,3)); z_max = max(X(:,3)); +%plotting tetrahedron with smallest insphere +i = i(1); +s_vert(1,:) = [1,3,2]; s_vert(2,:) = [1,2,4]; s_vert(3,:) = [1,4,3]; s_vert(4,:) = [2,3,4]; +tmp = [X(tetra(i,:),1),X(tetra(i,:),2),X(tetra(i,:),3)]; +tx_min = min(tmp(:,1)); tx_max = max(tmp(:,1)); +ty_min = min(tmp(:,2)); ty_max = max(tmp(:,2)); +tz_min = min(tmp(:,3)); tz_max = max(tmp(:,3)); +tMAX = max(abs([tx_min,tx_max,ty_min,ty_max,tz_min,tz_max])); +figure(2) +subplot(121), hold on +for j = 1:4 + fill3(tmp(s_vert(j,:),1),tmp(s_vert(j,:),2),tmp(s_vert(j,:),3),'r'); +end +axis equal, grid on, axis([x_min x_max y_min y_max z_min z_max]); +xlabel('Location in x');ylabel('Location in y');zlabel('Location in z'); +subplot(122), hold on +for j = 1:4 + fill3(tmp(s_vert(j,:),1),tmp(s_vert(j,:),2),tmp(s_vert(j,:),3),'r'); +end +axis equal, grid on +axis([tx_min-0.1*tMAX tx_max+0.1*tMAX ty_min-0.1*tMAX ty_max+0.1*tMAX tz_min-0.1*tMAX tz_max+0.1*tMAX]); +xlabel('zoomed in x');ylabel('zoomed in y');zlabel('zoomed in z'); diff --git a/legacy_scripts/gambit_hex2tet.m b/legacy_scripts/gambit_hex2tet.m new file mode 100644 index 0000000..afda339 --- /dev/null +++ b/legacy_scripts/gambit_hex2tet.m @@ -0,0 +1,162 @@ +%% +% @file +% This file is part of SeisSol. +% +% @author Martin Kaeser (martin.kaeser AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/kaeser) +% +% @section LICENSE +% Copyright (c) 2005, 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 +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. + +home; +disp(' ') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' %% %%') +disp(' %% GAMBIT_HEX2TET %%') +disp(' %% Creates regular tetrahedral mesh %%') +disp(' %% on the base of an hexahedral mesh %%') +disp(' %% %%') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' ') + +clear, close all; + +filename = input(' Filename(suffix ".neu" assumed) : ','s'); +disp(' ') +nx = input(' Give number of elements in x-dimension : '); +ny = input(' Give number of elements in y-dimension : '); +nz = input(' Give number of elements in z-dimension : '); + +disp(' ') +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\t Reading data from: \t\t%s',[filename,'.neu'])); + +fid_in = fopen([filename ,'.neu']); +fid_out = fopen([filename ,'_tetra.neu'],'w'); +junk = fgetl(fid_in); fprintf(fid_out,'%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +param = fscanf(fid_in,'%i',[6,1]); NX = param(1); NT = param(2); +param(2) = param(2)*5; fprintf(fid_out,'\n%10.0f%10.0f%10.0f%10.0f%10.0f%10.0f',param); +junk = fgetl(fid_in); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +X = fscanf(fid_in,'%g',[4,NX]); + fprintf(fid_out,'\n%10.0f%20.10e%20.10e%20.10e',X); +X = X'; X(:,1) = []; +junk = fgetl(fid_in); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +hex = fscanf(fid_in,'%g',[11,NT]); hex = hex'; hex(:,1:3) = []; +junk = fgetl(fid_in); +junk = fgetl(fid_in); +junk = fgetl(fid_in); +fclose(fid_in); + +disp(sprintf('\t Read \t\t\t\t\t%i vertices',NX)); +disp(sprintf('\t Read \t\t\t\t\t%i hexahedral elements',NT)); + +xmin = min(X(:,1)); xmax = max(X(:,1)); +ymin = min(X(:,2)); ymax = max(X(:,2)); +zmin = min(X(:,3)); zmax = max(X(:,3)); + +dx = (xmax-xmin)/nx; +dy = (ymax-ymin)/ny; +dz = (zmax-zmin)/nz; + +disp(' ') +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\t Creating tetraheral elemenst ... \n')); +t = cputime; + +tet = zeros(NT*5,7); +count = 0; +for i = 1:NT + + center = sum([X(hex(i,:),1), X(hex(i,:),2), X(hex(i,:),3)])/8; + + xpos = ceil( (center(1)-xmin)/dx ); + ypos = ceil( (center(2)-ymin)/dy ); + zpos = ceil( (center(3)-zmin)/dz ); + + if( (mod(xpos,2)==1 & mod(ypos,2)==1 & mod(zpos,2)==1) | ... + (mod(xpos,2)==0 & mod(ypos,2)==0 & mod(zpos,2)==1) | ... + (mod(xpos,2)==1 & mod(ypos,2)==0 & mod(zpos,2)==0) | ... + (mod(xpos,2)==0 & mod(ypos,2)==1 & mod(zpos,2)==0) ) + + count = count+1; + tet(count,:) = [count,6,4,hex(i,1),hex(i,2),hex(i,3),hex(i,5)]; + count = count+1; + tet(count,:) = [count,6,4,hex(i,6),hex(i,2),hex(i,5),hex(i,8)]; + count = count+1; + tet(count,:) = [count,6,4,hex(i,4),hex(i,3),hex(i,2),hex(i,8)]; + count = count+1; + tet(count,:) = [count,6,4,hex(i,7),hex(i,5),hex(i,3),hex(i,8)]; + count = count+1; + tet(count,:) = [count,6,4,hex(i,3),hex(i,5),hex(i,2),hex(i,8)]; + + else + + count = count+1; + tet(count,:) = [count,6,4,hex(i,8),hex(i,7),hex(i,4),hex(i,6)]; + count = count+1; + tet(count,:) = [count,6,4,hex(i,3),hex(i,1),hex(i,4),hex(i,7)]; + count = count+1; + tet(count,:) = [count,6,4,hex(i,2),hex(i,1),hex(i,6),hex(i,4)]; + count = count+1; + tet(count,:) = [count,6,4,hex(i,5),hex(i,6),hex(i,1),hex(i,7)]; + count = count+1; + tet(count,:) = [count,6,4,hex(i,4),hex(i,6),hex(i,7),hex(i,1)]; + + end + +end + +disp(sprintf('\t Writing output-file: \t\t\t%s',[filename,'_tetra.neu'])); + +fprintf(fid_out,'\n%8.0f%3.0f%3.0f%9.0f%8.0f%8.0f%8.0f',tet'); +fprintf(fid_out,'\n%s','ENDOFSECTION'); +fprintf(fid_out,'\n%s',' ELEMENT GROUP 2.0.0'); +fprintf(fid_out,'\n%s','GROUP: 1 ELEMENTS:'); +fprintf(fid_out,'%11.0f%',NT*5); +fprintf(fid_out,'%s',' MATERIAL: 2 NFLAGS: 1'); +fprintf(fid_out,'\n%s',' fluid'); +fprintf(fid_out,'\n%s',' 0'); +fprintf(fid_out,'\n%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f',(1:1:NT*5)); +fprintf(fid_out,'\n%s\n','ENDOFSECTION'); + +fclose(fid_out); + +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\n\t Conversion finished successfully! (%g CPU sec)\n',cputime-t)); +disp('-----------------------------------------------------------------------------------') +disp(' ') diff --git a/legacy_scripts/gambit_hex2tet_with_BC.m b/legacy_scripts/gambit_hex2tet_with_BC.m new file mode 100644 index 0000000..5a72f2e --- /dev/null +++ b/legacy_scripts/gambit_hex2tet_with_BC.m @@ -0,0 +1,349 @@ +%% +% @file +% This file is part of SeisSol. +% +% @author Martin Kaeser (martin.kaeser AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/kaeser) +% @author Thomas Ulrich (thomas.ulrich AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/ulrich) +% +% @section LICENSE +% Copyright (c) 2005, 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 +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. + + + + +home; +disp(' ') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' %% %%') +disp(' %% GAMBIT_HEX2TET_with_BC %%') +disp(' %% Creates regular tetrahedral mesh %%') +disp(' %% on the base of an hexahedral mesh %%') +disp(' %% Then add boundary conditions %%') +disp(' %% useful for complex geometries %%') +disp(' %% i.e. when gambit cannot recreate the geometry %%') +disp(' %% from the mesh %%') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' ') + +clear, close all; + +filename = input(' Filename(suffix ".neu" assumed) : ','s'); +disp(' ') +nx = input(' Give number of elements in x-dimension : '); +ny = input(' Give number of elements in y-dimension : '); +nz = input(' Give number of elements in z-dimension : '); + +disp(' ') +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\t Reading data from: \t\t%s',[filename,'.neu'])); + +fid_in = fopen([filename ,'.neu']); +fid_out = fopen([filename ,'_tetra.neu'],'w'); +junk = fgetl(fid_in); fprintf(fid_out,'%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +param = fscanf(fid_in,'%i',[6,1]); NX = param(1); NT = param(2); +param(4)=3; +param(2) = param(2)*5; fprintf(fid_out,'\n%10.0f%10.0f%10.0f%10.0f%10.0f%10.0f',param); +junk = fgetl(fid_in); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +X = fscanf(fid_in,'%g',[4,NX]); + fprintf(fid_out,'\n%10.0f%20.10e%20.10e%20.10e',X); +X = X'; X(:,1) = []; +junk = fgetl(fid_in); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +hex = fscanf(fid_in,'%g',[11,NT]); hex = hex'; hex(:,1:3) = []; +junk = fgetl(fid_in); +junk = fgetl(fid_in); +junk = fgetl(fid_in); +fclose(fid_in); + +disp(sprintf('\t Read \t\t\t\t\t%i vertices',NX)); +disp(sprintf('\t Read \t\t\t\t\t%i hexahedral elements',NT)); + +xmin = min(X(:,1)); xmax = max(X(:,1)); +ymin = min(X(:,2)); ymax = max(X(:,2)); +zmin = min(X(:,3)); zmax = max(X(:,3)); + +dx = (xmax-xmin)/nx; +dy = (ymax-ymin)/ny; +dz = (zmax-zmin)/nz; + +disp(' ') +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\t Creating tetraheral elements ... \n')); +t = cputime; + +tet = zeros(NT*5,7); +count = 0; + +b101=[]; +b103=[]; +b105=[]; + +for i = 1:NT + + center = sum([X(hex(i,:),1), X(hex(i,:),2), X(hex(i,:),3)])/8; + + xpos = ceil( (center(1)-xmin)/dx ); + ypos = ceil( (center(2)-ymin)/dy ); + zpos = ceil( (center(3)-zmin)/dz ); + + %mx (="minus x" taken as an example) will contains + %(1) the id of the tetra elements having a face in x=xmin. + %(2) the face id of these faces (between 1 and 4). + %This information is used for setting up the boundary conditions + + mx=[]; + my=[]; + mz=[]; + px=[]; + py=[]; + pz=[]; + + %hex2 is a permutation of the index of the hexahedron in order to put + %it in a(n) (arbitrary) situation of reference + % 4---8 + % 3---7 ! + % ! ! + % 2---6 + % 1---5 + % + %then we have the face id: top=pz=3, down=mz=2, my=6, py=5 mx=4,px=2 + + + hex2=zeros(8,2); + for j = 1:8 + hex2(j,2)=j; + xj=X(hex(i,j),1); + yj=X(hex(i,j),2); + zj=X(hex(i,j),3); + + if (xj-center(1)<0) + if (yj-center(2)<0) + if (zj-center(3)<0) + hex2(j,1)=5; + else + hex2(j,1)=7; + end + else + if (zj-center(3)<0) + hex2(j,1)=1; + else + hex2(j,1)=3; + end + end + else + if (yj-center(2)<0) + if (zj-center(3)<0) + hex2(j,1)=6; + else + hex2(j,1)=8; + end + else + if (zj-center(3)<0) + hex2(j,1)=2; + else + hex2(j,1)=4; + end + end + end + end + hex2 = sortrows(hex2,1); + hex2=hex2(:,2); + + if( (mod(xpos,2)==1 & mod(ypos,2)==1 & mod(zpos,2)==1) | ... + (mod(xpos,2)==0 & mod(ypos,2)==0 & mod(zpos,2)==1) | ... + (mod(xpos,2)==1 & mod(ypos,2)==0 & mod(zpos,2)==0) | ... + (mod(xpos,2)==0 & mod(ypos,2)==1 & mod(zpos,2)==0) ) + + count = count+1; + tet(count,:) = [count,6,4,hex(i,hex2(1)),hex(i,hex2(2)),hex(i,hex2(3)),hex(i,hex2(5))]; + mx=[mx;[count,surfid('134')]]; + py=[py;[count,surfid('123')]]; + mz=[mz;[count,surfid('124')]]; + + count = count+1; + tet(count,:) = [count,6,4,hex(i,hex2(6)),hex(i,hex2(2)),hex(i,hex2(5)),hex(i,hex2(8))]; + px=[px;[count,surfid('124')]]; + my=[my;[count,surfid('134')]]; + mz=[mz;[count,surfid('123')]]; + + count = count+1; + tet(count,:) = [count,6,4,hex(i,hex2(4)),hex(i,hex2(3)),hex(i,hex2(2)),hex(i,hex2(8))]; + px=[px;[count,surfid('134')]]; + py=[py;[count,surfid('123')]]; + pz=[pz;[count,surfid('124')]]; + + count = count+1; + tet(count,:) = [count,6,4,hex(i,hex2(7)),hex(i,hex2(5)),hex(i,hex2(3)),hex(i,hex2(8))]; + mx=[mx;[count,surfid('123')]]; + my=[my;[count,surfid('124')]]; + pz=[pz;[count,surfid('134')]]; + + count = count+1; + tet(count,:) = [count,6,4,hex(i,hex2(3)),hex(i,hex2(5)),hex(i,hex2(2)),hex(i,hex2(8))]; + + else + + count = count+1; + tet(count,:) = [count,6,4,hex(i,hex2(8)),hex(i,hex2(7)),hex(i,hex2(4)),hex(i,hex2(6))]; + px=[px;[count,surfid('134')]]; + my=[my;[count,surfid('124')]]; + pz=[pz;[count,surfid('123')]]; + + count = count+1; + tet(count,:) = [count,6,4,hex(i,hex2(3)),hex(i,hex2(1)),hex(i,hex2(4)),hex(i,hex2(7))]; + mx=[mx;[count,surfid('124')]]; + py=[py;[count,surfid('123')]]; + pz=[pz;[count,surfid('134')]]; + + count = count+1; + tet(count,:) = [count,6,4,hex(i,hex2(2)),hex(i,hex2(1)),hex(i,hex2(6)),hex(i,hex2(4))]; + px=[px;[count,surfid('134')]]; + py=[py;[count,surfid('124')]]; + mz=[mz;[count,surfid('123')]]; + + count = count+1; + tet(count,:) = [count,6,4,hex(i,hex2(5)),hex(i,hex2(6)),hex(i,hex2(1)),hex(i,hex2(7))]; + mx=[mx;[count,surfid('134')]]; + my=[my;[count,surfid('124')]]; + mz=[mz;[count,surfid('123')]]; + + count = count+1; + tet(count,:) = [count,6,4,hex(i,hex2(4)),hex(i,hex2(6)),hex(i,hex2(7)),hex(i,hex2(1))]; + + end + %free surface + if (zpos==nz) + b101=vertcat(b101,pz); + end + + %rupture + %if (zpos>5)&&(xpos>5)&&(xpos<=45) + if (ypos==ny/2) + b103=vertcat(b103,py); + end + if (ypos==(ny/2+1)) + b103=vertcat(b103,my); + end + %end + %absorbing surfaces + if (ypos==1) + b105=vertcat(b105,my); + end + if (ypos==ny) + b105=vertcat(b105,py); + end + if (xpos==1) + b105=vertcat(b105,mx); + end + if (xpos==nx) + b105=vertcat(b105,px); + end + if (zpos==1) + b105=vertcat(b105,mz); + end +end + +disp(sprintf('\t Writing output-file: \t\t\t%s',[filename,'_tetra.neu'])); + +fprintf(fid_out,'\n%8.0f%3.0f%3.0f%9.0f%8.0f%8.0f%8.0f',tet'); +fprintf(fid_out,'\n%s','ENDOFSECTION'); +fprintf(fid_out,'\n%s',' ELEMENT GROUP 2.0.0'); +fprintf(fid_out,'\n%s','GROUP: 1 ELEMENTS:'); +fprintf(fid_out,'%11.0f%',NT*5); +fprintf(fid_out,'%s',' MATERIAL: 2 NFLAGS: 1'); +fprintf(fid_out,'\n%s',' fluid'); +fprintf(fid_out,'\n%s',' 0'); +fprintf(fid_out,'\n%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f',(1:1:NT*5)); +fprintf(fid_out,'\n%s\n','ENDOFSECTION'); + +fprintf(fid_out,' BOUNDARY CONDITIONS 2.3.16\n'); +aSize=size(b101); +nbound=aSize(1); +fprintf(fid_out,'%s%8.0f%8.0f%8.0f%8.0f\n',' 101',1,nbound,0,6); +for i = 1:nbound + fprintf(fid_out,'%10.0f%5.0f%5.0f\n', b101(i,1),6,b101(i,2)); +end +fprintf(fid_out,'%s\n','ENDOFSECTION'); + +fprintf(fid_out,' BOUNDARY CONDITIONS 2.3.16\n'); +aSize=size(b103); +nbound=aSize(1); +fprintf(fid_out,'%s%8.0f%8.0f%8.0f%8.0f\n',' 103',1,nbound,0,6); +for i = 1:nbound + fprintf(fid_out,'%10.0f%5.0f%5.0f\n', b103(i,1),6,b103(i,2)); +end +fprintf(fid_out,'%s\n','ENDOFSECTION'); + +fprintf(fid_out,' BOUNDARY CONDITIONS 2.3.16\n'); +aSize=size(b105); +nbound=aSize(1); +fprintf(fid_out,'%s%8.0f%8.0f%8.0f%8.0f\n',' 105',1,nbound,0,6); +for i = 1:nbound + fprintf(fid_out,'%10.0f%5.0f%5.0f\n', b105(i,1),6,b105(i,2)); +end +fprintf(fid_out,'%s\n','ENDOFSECTION'); + +fclose(fid_out); + +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\n\t Conversion finished successfully! (%g CPU sec)\n',cputime-t)); +disp('-----------------------------------------------------------------------------------') +disp(' ') + + +%here we can plot the boundary condition to be sure that everything is +%alright + +plotmesh='n'; +tet=tet(:,4:7); +s_vert(1,:) = [1,3,2]; s_vert(2,:) = [1,2,4]; s_vert(3,:) = [2,3,4]; s_vert(4,:) = [1,4,3]; + +aSize=size(b105); +NS=aSize(1); + +tri=zeros(NS,3); +for i = 1:NS + tri(i,:)=tet(b105(i,1),s_vert(b105(i,2),:))'; +end + +if plotmesh == 'y' + figure; hold on + trimesh(tri,X(:,1),X(:,2),X(:,3),X(:,3)); + axis equal,xlabel('x'),ylabel('y'),zlabel('z') +end + diff --git a/legacy_scripts/gambit_redrefine.m b/legacy_scripts/gambit_redrefine.m new file mode 100644 index 0000000..435050e --- /dev/null +++ b/legacy_scripts/gambit_redrefine.m @@ -0,0 +1,207 @@ +%% +% @file +% This file is part of SeisSol. +% +% @author Martin Kaeser (martin.kaeser AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/kaeser) +% +% @section LICENSE +% Copyright (c) 2005, 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 +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. +home; +disp(' ') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' %% %%') +disp(' %% GAMBIT_RED_REFINE (2D) %%') +disp(' %% %%') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' ') + +clear, close all; + +filename = input(' Filename(suffix ".neu" assumed): ','s'); +disp(' ') +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\t Reading data from: \t\t\t%s',[filename,'.neu'])); +disp(' ') +fid_in = fopen([filename ,'.neu']); +fid_out = fopen([filename ,'_refined.neu'],'w'); +junk = fgetl(fid_in); fprintf(fid_out,'%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +junk = fgetl(fid_in); fprintf(fid_out,'\n%s',junk); +param = fscanf(fid_in,'%i',[6,1]); +NX = param(1); NT = param(2); NG = param(3); NBSETS = param(4); NDFCD = param(5); NDFVL = param(6); +junk = fgetl(fid_in); +junk = fgetl(fid_in); +junk = fgetl(fid_in); +X = fscanf(fid_in,'%g',[3,NX]); X = X'; X(:,1) = []; +junk = fgetl(fid_in); +junk = fgetl(fid_in); +junk = fgetl(fid_in); +tri = fscanf(fid_in,'%g',[6,NT]); tri = tri'; tri(:,1:3) = []; +junk = fgetl(fid_in); +junk = fgetl(fid_in); +junk = fgetl(fid_in); +junk = fgetl(fid_in); +junk = fgetl(fid_in); +junk = fgetl(fid_in); +junk = fscanf(fid_in,'%g',[10,ceil(NT/10)]); +%junk = fgetl(fid_in); +%junk = fgetl(fid_in); +%param = fscanf(fid_in,'%i',[5,1]); +%bc = fscanf(fid_in,'%g',[3,param(3)]); bc = bc'; +%junk = fgetl(fid_in); +%junk = fgetl(fid_in); +fclose(fid_in); + +disp(sprintf('\t Read \t\t\t\t\t%i vertices',NX)); +disp(sprintf('\t Read \t\t\t\t\t%i triangular elements',NT)); + +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\t Start refinement ...')); + +subplot(121) +h = trimesh(tri,X(:,1),X(:,2),zeros(length(X),1));view(0,90),axis square +set(h,'EdgeColor','k') +xlabel('x','FontWeight','bold','FontSize',12) +ylabel('y','FontWeight','bold','FontSize',12,'Rotation',0) +axis([-1 1 -1 1 0 1]), view(0,90), axis square, drawnow; + +tin=cputime; + +MX = [(X(tri(:,1),1)+X(tri(:,2),1)),(X(tri(:,2),1)+X(tri(:,3),1)),(X(tri(:,3),1)+X(tri(:,1),1))]/2; +MY = [(X(tri(:,1),2)+X(tri(:,2),2)),(X(tri(:,2),2)+X(tri(:,3),2)),(X(tri(:,3),2)+X(tri(:,1),2))]/2; +X_ref = []; +bc_ref = []; +c=0; +for t = 1:NT + for n = 1:4 + c = c+1; + switch n + case 1 + x = [X(tri(t,1),:);MX(t,1),MY(t,1);MX(t,3),MY(t,3)]; + case 2 + x = [X(tri(t,2),:);MX(t,2),MY(t,2);MX(t,1),MY(t,1)]; + case 3 + x = [X(tri(t,3),:);MX(t,3),MY(t,3);MX(t,2),MY(t,2)]; + case 4 + x = [MX(t,1),MY(t,1);MX(t,2),MY(t,2);MX(t,3),MY(t,3)]; + end + [tmp,io,in] = intersect(X_ref,x,'rows'); + len = size(X_ref,1); + if length(in) > 0; + x(in,:) = []; + ind = [io',len+1,len+2]; + tri_ref(c,1:3) = ind(1:3); + X_ref = [X_ref;x]; + else + tri_ref(c,1:3) = (len+1:len+3); + X_ref = [X_ref;x]; + end + % compute counter-clockwise oriented edges + p = [X_ref(tri_ref(c,1),:),X_ref(tri_ref(c,2),:),X_ref(tri_ref(c,3),:)]; + % compute cross product + cr = (p(3)-p(1)).*(p(6)-p(2))-(p(4)-p(2)).*(p(5)-p(1)); + % reverse order of vertices where necessary + tri_ref(c,:) = fliplr(tri_ref(c,:)); + + end +% bcind = find(bc(:,1)==t); +% if(length(bcind)>0) +% bedge = bc(bcind,3); +% switch bedge +% case 1 +% bc_new = [c-3;c-2]; +% evec = X(tri(t,2),:)-X(tri(t,1),:); +% case 2 +% bc_new = [c-2;c-1]; +% evec = X(tri(t,3),:)-X(tri(t,2),:); +% case 3 +% bc_new = [c-1;c-3]; +% evec = X(tri(t,1),:)-X(tri(t,3),:); +% end +% +% bc_ref = [bc_ref; [bc_new,[3;3],[1;3]]]; +% end + +end + +tri = tri_ref; X = X_ref; +[NX,tmp] = size(X); +[NT,tmp] = size(tri); +disp(sprintf('\t Writing output-file: \t\t\t%s',[filename,'_refined.neu'])); +disp(sprintf('\n\t Write \t\t\t\t\t%i vertices',NX)); +disp(sprintf('\t Write \t\t\t\t\t%i triangular elements',NT)); + +fprintf(fid_out,'\n%10.0f%10.0f%10.0f%10.0f%10.0f%10.0f',NX,NT,NG,NBSETS,NDFCD,NDFVL); +fprintf(fid_out,'\n%s','ENDOFSECTION'); +fprintf(fid_out,'\n%s',' NODAL COORDINATES 2.0.0'); +fprintf(fid_out,'\n%10.0f%20.10e%20.10e',[(1:NX)',X]'); +fprintf(fid_out,'\n%s','ENDOFSECTION'); +fprintf(fid_out,'\n%s',' ELEMENTS/CELLS 2.0.0'); +fprintf(fid_out,'\n%8.0f%3.0f%3.0f%9.0f%8.0f%8.0f',[(1:NT)',ones(NT,1)*3,ones(NT,1)*3,tri]'); +fprintf(fid_out,'\n%s','ENDOFSECTION'); +fprintf(fid_out,'\n%s',' ELEMENT GROUP 2.0.0'); +fprintf(fid_out,'\n%s','GROUP: 1 ELEMENTS:'); +fprintf(fid_out,'%11.0f%',NT); +fprintf(fid_out,'%s',' MATERIAL: 2 NFLAGS: 1'); +fprintf(fid_out,'\n%s',' fluid'); +fprintf(fid_out,'\n%s',' 0'); +fprintf(fid_out,'\n%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f',(1:1:NT)); +fprintf(fid_out,'\n%s\n','ENDOFSECTION'); +%fprintf(fid_out,'\n%s',' BOUNDARY CONDITIONS 2.0.0'); +%fprintf(fid_out,'\n%32.0f%8.0f%8.0f%8.0f%8.0f',param(1),param(2),param(3)*2,param(4),param(5)); +%fprintf(fid_out,'\n%10.0f%5.0f%5.0f',bc_ref'); +%fprintf(fid_out,'\n%s\n','ENDOFSECTION'); +fclose(fid_out); +subplot(122) +h = trimesh(tri,X(:,1),X(:,2),zeros(length(X),1));view(0,90),axis square +%set(gca,'XTick',[-0.5 0 0.5]),set(gca,'YTick',[-0.5 0 0.5]) +%set(gca,'XTickLabel',[-0.5 0 0.5],'FontSize',11,'FontWeight','bold') +%set(gca,'YTickLabel',[-0.5 0 0.5],'FontSize',11,'FontWeight','bold') +%xlabel('x','FontSize',12,'FontWeight','bold') +%ylabel('y','FontSize',12,'FontWeight','bold','Rotation',0) +%save tri.mat tri; save X.mat X; +set(h,'EdgeColor','k') +xlabel('x','FontWeight','bold','FontSize',12) +ylabel('y','FontWeight','bold','FontSize',12,'Rotation',0) +%set(gca,'XTick',[-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5],... +% 'XTickLabels',[-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5]) +%set(gca,'YTick',[-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5],... +% 'YTickLabels',[-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5]) +axis([-1 1 -1 1 0 1]), view(0,90), axis square, drawnow; + +tt = cputime-tin; + +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\n\t Refinement finished successfully! (%g CPU sec)\n',tt)); +disp('-----------------------------------------------------------------------------------') +disp(' ') diff --git a/legacy_scripts/matlab2gambit.m b/legacy_scripts/matlab2gambit.m new file mode 100644 index 0000000..9b3640b --- /dev/null +++ b/legacy_scripts/matlab2gambit.m @@ -0,0 +1,105 @@ +%% +% @file +% This file is part of SeisSol. +% +% @author Martin Kaeser (martin.kaeser AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/kaeser) +% +% @section LICENSE +% Copyright (c) 2006, 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 +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. +home; +disp(' ') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' %% %%') +disp(' %% MATLAB2GAMBIT %%') +disp(' %% %%') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' ') + +clear, close all; + +filename = input(' MATLAB filename: ','s'); + +%%%% READ IN MATLAB MESH FILE +fid = fopen([filename,'.msh']); +tmp = fscanf(fid,'%i',[2,1]); n_vert = tmp(1); n_tri = tmp(2); +tri = fscanf(fid,'%i',[3,n_tri]); tri = tri'; +X = fscanf(fid,'%g',[3,n_vert]); X = X'; +fclose(fid); + +fid_out = fopen([filename,'.neu'],'w'); +fprintf(fid_out,'%s',' CONTROL INFO 2.0.0'); +fprintf(fid_out,'\n%s','** GAMBIT NEUTRAL FILE'); +fprintf(fid_out,'\n%s',filename); +fprintf(fid_out,'\n%s','PROGRAM: Gambit VERSION: 2.0.0'); +fprintf(fid_out,'\n%s',date); +fprintf(fid_out,'\n%s',' NUMNP NELEM NGRPS NBSETS NDFCD NDFVL'); +fprintf(fid_out,'\n%10.0f%10.0f%10.0f%10.0f%10.0f%10.0f',n_vert,n_tri,1,1,2,2); +fprintf(fid_out,'\n%s','ENDOFSECTION'); +fprintf(fid_out,'\n%s',' NODAL COORDINATES 2.0.0'); +counter = (1:n_vert)'; +fprintf(fid_out,'\n%10.0f%20.10e%20.10e',[counter,X(:,1:2)]'); +fprintf(fid_out,'\n%s','ENDOFSECTION'); +fprintf(fid_out,'\n%s',' ELEMENTS/CELLS 2.0.0'); +counter = (1:n_tri)'; +fprintf(fid_out,'\n%8.0f%3.0f%3.0f%9.0f%8.0f%8.0f',[counter,ones(size(counter))*3,ones(size(counter))*3,tri]'); +fprintf(fid_out,'\n%s','ENDOFSECTION'); +fprintf(fid_out,'\n%s',' ELEMENT GROUP 2.0.0'); +fprintf(fid_out,'\n%s','GROUP: 1 ELEMENTS:'); +fprintf(fid_out,'%11.0f%',n_tri); +fprintf(fid_out,'%s',' MATERIAL: 2 NFLAGS: 1'); +fprintf(fid_out,'\n%s',' solid'); +fprintf(fid_out,'\n%s',' 0'); +fprintf(fid_out,'\n%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f%8.0f',(1:1:n_tri)); +fprintf(fid_out,'\n%s','ENDOFSECTION'); +fprintf(fid_out,'\n%s',' BOUNDARY CONDITIONS 2.0.0'); + +trimesh(tri,X(:,1),X(:,2),zeros(n_vert,1)), view(0,90), axis square, hold on +b = find(X(:,3)==101); +plot(X(b,1),X(b,2),'r.') + +Btri = ismember(tri,b); +sumtri = sum(Btri,2); +b = find(sumtri==2); +sumedge(:,1) = sum(Btri(b,[1,2]),2); +sumedge(:,2) = sum(Btri(b,[2,3]),2); +sumedge(:,3) = sum(Btri(b,[1,3]),2); +[tmp,edge] = find(sumedge==2); +tmp = sortrows([tmp,edge]); +edge = tmp(:,2); +CX = ( X(tri(:,1),1)+X(tri(:,2),1)+X(tri(:,3),1) ) / 3; +CY = ( X(tri(:,1),2)+X(tri(:,2),2)+X(tri(:,3),2) ) / 3; +plot(CX(b),CY(b),'b.') + +fprintf(fid_out,'\n%32.0f%8.0f%8.0f%8.0f%8.0f',101,1,length(b),0,6); +fprintf(fid_out,'\n%10.0f%5.0f%5.0f',[b,ones(size(b))*3,edge]'); +fprintf(fid_out,'\n%s\n','ENDOFSECTION'); + + +fclose(fid_out); diff --git a/legacy_scripts/mesh_mod.f90 b/legacy_scripts/mesh_mod.f90 new file mode 100644 index 0000000..3565020 --- /dev/null +++ b/legacy_scripts/mesh_mod.f90 @@ -0,0 +1,440 @@ +!> +!! @file +!! This file is part of SeisSol. +!! +!! @author Gilbert B. Brietzke (Gilbert.Brietzke AT lrz.de, http://www.geophysik.uni-muenchen.de/Members/brietzke) +!! +!! @section LICENSE +!! Copyright (c) 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 +!! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +!! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!! POSSIBILITY OF SUCH DAMAGE. +module mesh_mod + implicit none + + integer,parameter, private :: dp=kind(1.d0) + + type stringpair + character(:), allocatable :: str,val + end type stringpair + + type sectionlabel + character(:), allocatable :: label,version,endlabel + contains + procedure :: read => read_label + procedure :: end => read_endlabel + procedure :: write => write_label + end type sectionlabel + + type header + type(sectionlabel):: seclab + type(stringpair) :: file + type(stringpair) :: soft + type(stringpair) :: vers + type(stringpair) :: date + integer :: NUMNP, NELEM, NGRPS, NBSETS, NDFCD, NDFVL + contains + procedure :: read => read_header + procedure :: write => write_header + procedure :: print => print_header + end type header + + type nodetype + integer :: numb + real(dp) :: crds(3) + end type nodetype + + type celltype + integer :: numb,dims,nvert + integer :: ngrp + integer, allocatable :: nodes(:) + end type celltype + + type grouptype + integer :: id + integer :: ncells + integer :: nprop + integer :: nflags + character(len=100) :: name + integer, dimension(:), allocatable :: cell_id + end type grouptype + + type bctype + character(len=100) :: label + integer :: n1,n,n2,n3 + integer, dimension(:), allocatable :: id + integer, dimension(:), allocatable :: dim + integer, dimension(:), allocatable :: side + end type bctype + + type meshtype + type(header) :: head + integer,pointer :: NUMNP, NELEM, NGRPS, NBSETS, NDFCD, NDFVL + type(nodetype) , allocatable :: nodes(:) + type(celltype) , allocatable :: cells(:) + type(grouptype), allocatable :: groups(:) + type(bctype) , allocatable :: bcs(:) + character(len=100) :: fileinfo + contains + procedure :: read_gambit => mesh_read_gambit + procedure :: write_gambit => mesh_write_gambit + end type meshtype + + private int2str + + type(sectionlabel):: seclab + +contains + + subroutine mesh_read_gambit(mesh,meshfile) + implicit none + class(meshtype),target :: mesh + character(len=*), intent(in) ::meshfile + type(nodetype) :: zeronode + type(celltype) :: zerocell + type(grouptype) :: zerogroup + integer :: i,j,ic=0,ibc=0,iloop,iun,iostat + write(*,'(a)')'opening meshfile: '//trim(adjustl(meshfile))//' ... : ' + open(newunit=iun,file=trim(adjustl(meshfile)),status='old',action='read') + call mesh%head%read(iun) + call mesh%head%print + mesh%NUMNP => mesh%head%NUMNP; mesh%NELEM => mesh%head%NELEM + mesh%NGRPS => mesh%head%NGRPS; mesh%NBSETS=> mesh%head%NBSETS + mesh%NDFCD => mesh%head%NDFCD; mesh%NDFVL => mesh%head%NDFVL + zeronode=nodetype(0,[0.d0,0.d0,0.d0]); + zerocell=celltype(0,0,0,0,null()) + zerogroup=grouptype(0,0,0,0,'',null()) + allocate(mesh%cells(mesh%nelem));mesh%cells=zerocell + allocate(mesh%groups(mesh%ngrps));mesh%groups=zerogroup + allocate(mesh%bcs(mesh%nbsets)); + allocate(mesh%nodes(mesh%numnp));mesh%nodes=zeronode + write(*,'(A)')'now looping sections:' + do iloop = 1,10000 + call seclab%read(iun) + if(seclab%label(1:17)=='NODAL COORDINATES') call mesh_read_nodes(mesh,iun) + if(seclab%label(1:14)=='ELEMENTS/CELLS') call mesh_read_elements(mesh,iun) + if(seclab%label(1:13)=='ELEMENT GROUP') call mesh_read_groups(mesh,iun) + if(seclab%label(1:19)=='BOUNDARY CONDITIONS')call mesh_read_bcs(mesh,iun) + cycle + enddo + close(iun) + end subroutine mesh_read_gambit + !------------------------------------------------------------------------------------- + subroutine mesh_write_gambit(mesh,meshfile,filename) + implicit none + class(meshtype) :: mesh + character(len=*), intent(in) ::meshfile + character(len=*), intent(in) ::filename + character(len=100) :: myform + integer :: i,j,iun,ibc + write(*,*)' opening meshfile: ',trim(adjustl(meshfile)) + open(newunit=iun,file=trim(adjustl(meshfile)),action='write') + call mesh%head%write(iun) + call mesh_write_nodes(mesh,iun) + call mesh_write_elements(mesh,iun) + call mesh_write_groups(mesh,iun) + call mesh_write_bcs(mesh,iun) + close(iun) + end subroutine mesh_write_gambit + + subroutine mesh_read_nodes(mesh,iun) + implicit none + type(meshtype) :: mesh + integer :: iun,i + write(*,'(" ",a,a,i10,a)',advance='no')trim(seclab%label),', reading: ',mesh%numnp,' mesh-nodes' + do i = 1,mesh%numnp + read(iun,*) mesh%nodes(i)%numb,mesh%nodes(i)%crds(1:mesh%ndfcd) + enddo + write(*,'(a,a)')', done! ',seclab%end(iun) + end subroutine mesh_read_nodes + + subroutine mesh_read_elements(mesh,iun) + implicit none + type(meshtype) :: mesh + integer :: iun,i,iostat + character(len=100) :: info + write(*,'(" ",a,a,i10,a)',advance='no')trim(seclab%label),', reading: ',mesh%nelem,' mesh-elements' + do i = 1,mesh%nelem + read(iun,'(A100)')info + read(info,*)mesh%cells(i)%numb,mesh%cells(i)%dims,mesh%cells(i)%nvert + allocate(mesh%cells(i)%nodes(mesh%cells(i)%nvert));mesh%cells(i)%nodes=0; + read(info,*,iostat=iostat) mesh%cells(i)%numb,mesh%cells(i)%dims, & + & mesh%cells(i)%nvert,mesh%cells(i)%nodes + if(.not.(is_iostat_end(iostat).or.is_iostat_eor(iostat)))cycle + read(info,*) mesh%cells(i)%numb,mesh%cells(i)%dims,mesh%cells(i)%nvert, & + & mesh%cells(i)%nodes(1:min(mesh%cells(i)%nvert,7)) + if(mesh%cells(i)%nvert>7)read(iun,*) mesh%cells(i)%nodes(8:mesh%cells(i)%nvert) + enddo + write(*,'(a,a)')', done! ',seclab%end(iun) + end subroutine mesh_read_elements + + subroutine mesh_read_groups(mesh,iun) + implicit none + type(meshtype) :: mesh + integer :: iun,j,iostat + integer, save :: ic=0 + character(len=100) :: char1 + ic=ic+1 + if(ic<=mesh%ngrps)then + write(*,'(" ",a,a,i5,a,i5,a)',advance="no")trim(seclab%label),', reading: ',ic,' of ',mesh%ngrps,' groups' + read(iun,*)char1,mesh%groups(ic)%id ,char1,mesh%groups(ic)%ncells, & + & char1,mesh%groups(ic)%nprop,char1,mesh%groups(ic)%nflags + read(iun,'(A100)')mesh%groups(ic)%name + read(iun,*) + allocate(mesh%groups(ic)%cell_id(mesh%groups(ic)%ncells)) + do j = 1,mesh%groups(ic)%ncells,10 + read(iun,*)mesh%groups(ic)%cell_id(j:min(j+9,mesh%groups(ic)%ncells)) + enddo + mesh%cells(mesh%groups(ic)%cell_id)%ngrp = mesh%groups(ic)%id; + write(*,'(a,a)')', done! ',seclab%end(iun) + endif + end subroutine mesh_read_groups + + subroutine mesh_read_bcs(mesh,iun) + implicit none + type(meshtype) :: mesh + integer :: iun,j,iostat + integer, save :: ibc=0 + ibc=ibc+1 + if(ibc<=mesh%nbsets)then + write(*,'(" ",a,a,i5,a,i5,a)',advance="no")trim(seclab%label),', reading: ',ibc,' of ',mesh%nbsets,' BCs' + read(iun,*,iostat=iostat)mesh%bcs(ibc)%label,mesh%bcs(ibc)%n1,mesh%bcs(ibc)%n,mesh%bcs(ibc)%n2,mesh%bcs(ibc)%n3 + if(is_iostat_err(iostat)) stop'error during read!' + select case (mesh%bcs(ibc)%n3) + case(6) + allocate(mesh%bcs(ibc)%id(mesh%bcs(ibc)%n)) + allocate(mesh%bcs(ibc)%dim(mesh%bcs(ibc)%n)) + allocate(mesh%bcs(ibc)%side(mesh%bcs(ibc)%n)) + do j = 1,mesh%bcs(ibc)%n + read(iun,*,iostat=iostat)mesh%bcs(ibc)%id(j),mesh%bcs(ibc)%dim(j),mesh%bcs(ibc)%side(j) + if(is_iostat_err(iostat)) stop'error during read!' + enddo + case(24) + allocate(mesh%bcs(ibc)%id(mesh%bcs(ibc)%n)) + do j = 1,mesh%bcs(ibc)%n + read(iun,*,iostat=iostat)mesh%bcs(ibc)%id(j) + if(is_iostat_err(iostat)) stop'error during read!' + enddo + end select + write(*,'(a,a)')', done! ',seclab%end(iun) + endif + end subroutine mesh_read_bcs + + subroutine mesh_write_nodes(mesh,iun) + implicit none + type(meshtype) :: mesh + integer :: iun,j + write(*,*)' writing: ',mesh%numnp,' mesh-nodes' + write(iun,'(A)')" NODAL COORDINATES "//trim(mesh%head%vers%val) + do j = 1,mesh%numnp + write(iun,'(i10,3es20.11E3)')mesh%nodes(j)%numb,mesh%nodes(j)%crds(1:mesh%ndfcd) + enddo + write(iun,'("ENDOFSECTION")') + write(*,*)' ... ',mesh%numnp,' nodes done' + end subroutine mesh_write_nodes + + subroutine mesh_write_elements(mesh,iun) + implicit none + type(meshtype) :: mesh + integer :: iun,j + character(len=100) :: myform + write(*,*)' writing: ',mesh%nelem,' mesh-elements' + write(iun,'(A)')" ELEMENTS/CELLS "//trim(mesh%head%vers%val) + do j = 1,mesh%nelem + myform='(i8,i3,i3," ",'//int2str(mesh%cells(j)%nvert,1,'(i1)')//'i8)' + write(iun,trim(adjustl(myform)))mesh%cells(j)%numb,mesh%cells(j)%dims,mesh%cells(j)%nvert,mesh%cells(j)%nodes + enddo + write(iun,'("ENDOFSECTION")') + write(*,*)' ... ',mesh%nelem,' elements done' + end subroutine mesh_write_elements + + subroutine mesh_write_groups(mesh,iun) + type(meshtype) :: mesh + integer :: iun,i,j + if(mesh%ngrps>0)then + write(*,*)' writing: ',mesh%ngrps,' groups' + do i =1,mesh%ngrps + write(iun,'(A)')" ELEMENT GROUP"//trim(mesh%head%vers%val) + write(iun,'("GROUP:",i11," ELEMENTS:",i11," MATERIAL:",i11," NFLAGS:",i11)') & + & mesh%groups(i)%id,mesh%groups(i)%ncells,mesh%groups(i)%nprop,mesh%groups(i)%nflags + write(iun,'(a40)')mesh%groups(i)%name + write(iun,*)' 0' + do j = 1,mesh%groups(i)%ncells,10 + write(iun,'(10i8)')mesh%groups(i)%cell_id(j:min(j+9,mesh%groups(i)%ncells)) + enddo + write(iun,'("ENDOFSECTION")') + enddo + write(*,*)' ... ',mesh%ngrps,' groups done' + endif + end subroutine mesh_write_groups + + subroutine mesh_write_bcs(mesh,iun) + implicit none + type(meshtype) :: mesh + integer :: iun,ibc,j + if(mesh%nbsets>0)then + do ibc =1,mesh%nbsets + write(iun,'(A)')" BOUNDARY CONDITIONS "//trim(mesh%head%vers%val) + write(*,'(" ",a,i5,a,i5,a)',advance="no")'BOUNDARY CONDITIONS, writing: ',ibc,' of ',mesh%nbsets,' BCs' + write(iun,'(" ",A8,4i8)')trim(adjustl(mesh%bcs(ibc)%label)),mesh%bcs(ibc)%n1, & + & mesh%bcs(ibc)%n,mesh%bcs(ibc)%n2,mesh%bcs(ibc)%n3 + select case (mesh%bcs(ibc)%n3) + case(6) + do j = 1,mesh%bcs(ibc)%n + write(iun,'(I10,I5,I5)')mesh%bcs(ibc)%id(j),mesh%bcs(ibc)%dim(j),mesh%bcs(ibc)%side(j) + enddo + case(24) + do j = 1,mesh%bcs(ibc)%n + write(iun,'(I10)')mesh%bcs(ibc)%id(j) + enddo + end select + write(iun,'("ENDOFSECTION")') + write(*,'(a)')', done! ' + end do + endif + end subroutine mesh_write_bcs + + subroutine read_label(lab,iun) + class(sectionlabel) :: lab + integer :: iun + character(len=20) :: dum1,dum2 + integer :: iostat + read(iun,'(A20,A)',iostat=iostat)dum1,dum2; + if(is_iostat_end(iostat))return + lab%label=trim(adjustl(dum1));lab%version=trim(adjustl(dum2)) + end subroutine read_label + + subroutine write_label(lab,iun) + class(sectionlabel) :: lab + integer :: iun + write(iun,'(A20,A,A)')adjustr(lab%label),' ',lab%version + end subroutine write_label + + function read_endlabel(lab,iun)result(endlabel) + class(sectionlabel) :: lab + character(len=12):: dum1,endlabel + integer :: iun + read(iun,'(A12)')dum1; + endlabel='ENDOFSECTION' + if(trim(adjustl(dum1))=='ENDOFSECTION')then + lab%endlabel=endlabel + return + end if + endlabel='INCONSISTENT' + end function read_endlabel + + subroutine read_header(head,iun) + class(header) :: head + integer :: iun + character(len=100) :: dum1,dum2,dum3,dum4 + call head%seclab%read(iun) + read(iun,'(A100)')dum1; read(iun,'(A)')dum2; + head%file%str=trim(adjustl(dum1));head%file%val=trim(adjustl(dum2)) + read(iun,*)dum1,dum2,dum3,dum4 + head%soft%str=trim(adjustl(dum1));head%soft%val=trim(adjustl(dum2)) + head%vers%str=trim(adjustl(dum3));head%vers%val=trim(adjustl(dum4)) + read(iun,'(A100)')dum2; + head%date%str='date';head%date%val=trim(adjustl(dum2)) + read(iun,*) + read(iun,*)head%NUMNP,head%NELEM,head%NGRPS,head%NBSETS,head%NDFCD,head%NDFVL + write(*,'(A)')'reading of section: '//head%seclab%label//': check: '//head%seclab%end(iun) + end subroutine read_header + + subroutine print_header(head) + class(header) :: head + call head%write() + end subroutine print_header + + subroutine write_header(head,unit) + use iso_fortran_env + class(header) :: head + integer, optional :: unit + integer :: iun + iun=output_unit + if(present(unit))iun=unit + call head%seclab%write(iun) + write(iun,'(A)')adjustl(head%file%str) + write(iun,'(A)')adjustl(head%file%val); + write(iun,'(A8,A,A21,A,A12,A,A)') & + & adjustr(head%soft%str),' ',adjustl(head%soft%val),' ',& + & adjustr(head%vers%str),' ',adjustl(head%vers%val); + write(iun,'(" ",A)')adjustl(head%date%val); + write(iun,'(" NUMNP NELEM NGRPS NBSETS NDFCD NDFVL")') + write(iun,'(6i10)')head%NUMNP,head%NELEM,head%NGRPS,head%NBSETS,head%NDFCD,head%NDFVL + write(iun,'(A)')'ENDOFSECTION' + end subroutine write_header + + function int2str(n,nc,myform) + implicit none + integer :: n,nc + character(len=nc)::int2str + character(len=*)::myform + write(int2str,myform)n + end function int2str + + logical function is_iostat_err(iostat) + integer :: iostat + if(is_iostat_end(iostat))is_iostat_err=.false. + if(is_iostat_eor(iostat))is_iostat_err=.false. + end function is_iostat_err + +end module mesh_mod + + +!!$program test_mesh +!!$ +!!$ use mesh_mod +!!$ +!!$ type(meshtype) :: mesh +!!$ character(len=512)::arg +!!$ character(:),allocatable :: infile,outfile +!!$ integer :: n +!!$ +!!$ call get_command_argument(1,arg) +!!$ infile=trim(adjustl(arg));n=len(infile) +!!$ outfile=infile(1:n-4)//'_new'//infile(n-3:n) +!!$ call mesh%read_gambit(infile) +!!$ call mesh%write_gambit(outfile,'newinfo') +!!$ +!!$end program test_mesh + +!!$program test +!!$ +!!$ use testo +!!$ use mesh_mod +!!$ +!!$ type(header) :: head +!!$ +!!$ open(newunit=iun,file='scec_box40x20x30_small.neu') +!!$ call head%read(iun) +!!$ close(iun) +!!$ call head%print() +!!$ open(newunit=iun,file='scec_box40x20x30_small_new.neu') +!!$ call head%write(iun) +!!$ close(iun) +!!$ +!!$end program test diff --git a/legacy_scripts/surfid.m b/legacy_scripts/surfid.m new file mode 100644 index 0000000..561af77 --- /dev/null +++ b/legacy_scripts/surfid.m @@ -0,0 +1,13 @@ +function y = surfid(x) + %determine Gambit surface id from node numbering + %called by gambit_hex2tet_with_BC + if x=='123' + y=1; + elseif x=='124' + y=2; + elseif x=='134' + y=4; + else + error('unkown surface %s',x) + end +end