diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5d0ed8e196ade597db1eac8670fb5c1fdf40538f..aa76b3065d3e51585782db40206cd4cc04827c00 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -48,7 +48,6 @@ OPTION(WAVES_USE_HEAT "Compile heat module" ON)
OPTION(WAVES_USE_MIRRORS "Compile mirrors module" ON)
OPTION(WAVES_USE_KATOPTRON "Compile katoptron module" ON)
OPTION(WAVES_USE_PARAMS "Compile params module" ON)
-OPTION(WAVES_USE_SPH "Compile sph module" ON)
OPTION(WAVES_USE_TLNOS "Compile tlnos module" ON)
OPTION(WAVES_USE_WAVES "Compile waves module" ON)
# (de)activate extensions
@@ -181,10 +180,6 @@ IF(WAVES_USE_PARAMS)
ADD_SUBDIRECTORY( params )
ENDIF()
-IF(WAVES_USE_SPH)
- ADD_SUBDIRECTORY( sph )
-ENDIF()
-
IF(WAVES_USE_TLNOS)
ADD_SUBDIRECTORY( tlnos )
ENDIF()
diff --git a/sph/.gitignore b/sph/.gitignore
deleted file mode 100644
index 3a3a98aba4f9c46384d9357b5a51c5ea6e0fe352..0000000000000000000000000000000000000000
--- a/sph/.gitignore
+++ /dev/null
@@ -1,2 +0,0 @@
-louisB
-.metadata
diff --git a/sph/CMakeLists.txt b/sph/CMakeLists.txt
deleted file mode 100644
index 156e07c4fcdc3410cc3e6e9413998ceffd8b5793..0000000000000000000000000000000000000000
--- a/sph/CMakeLists.txt
+++ /dev/null
@@ -1,20 +0,0 @@
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-# Add source dir
-ADD_SUBDIRECTORY( src )
-ADD_SUBDIRECTORY( _src )
-
-# Add test dir
-MACRO_AddTest(${CMAKE_CURRENT_SOURCE_DIR}/tests)
diff --git a/sph/__init__.py b/sph/__init__.py
deleted file mode 100644
index e56cbd377a485f7887d9b0465a4811d25a5b8348..0000000000000000000000000000000000000000
--- a/sph/__init__.py
+++ /dev/null
@@ -1,6 +0,0 @@
-# -*- coding: utf-8 -*-
-# sph MODULE initialization file
-
-import fwk
-import tbox
-from sphw import *
diff --git a/sph/_src/CMakeLists.txt b/sph/_src/CMakeLists.txt
deleted file mode 100644
index 70b86a4ce5f01699dac9bc11b1c6b96f0ecacb74..0000000000000000000000000000000000000000
--- a/sph/_src/CMakeLists.txt
+++ /dev/null
@@ -1,49 +0,0 @@
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-# CMake input file of the SWIG wrapper around "sph.so"
-
-INCLUDE(${SWIG_USE_FILE})
-
-FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
-FILE(GLOB ISRCS *.i)
-
-SET(CMAKE_SWIG_FLAGS "")
-SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
-
-SET(SWINCFLAGS
--I${PROJECT_SOURCE_DIR}/sph/src
--I${PROJECT_SOURCE_DIR}/tbox/src
--I${PROJECT_SOURCE_DIR}/tbox/_src
--I${PROJECT_SOURCE_DIR}/fwk/src
--I${PROJECT_SOURCE_DIR}/fwk/_src
--I${EIGEN_INCLUDE_DIRS}
-)
-SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES SWIG_FLAGS "${SWINCFLAGS}")
-
-if (${CMAKE_VERSION} VERSION_LESS "3.8.0")
- SWIG_ADD_MODULE(sphw python ${ISRCS} ${SRCS})
-else()
- SWIG_ADD_LIBRARY(sphw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
-endif()
-MACRO_DebugPostfix(_sphw)
-
-TARGET_INCLUDE_DIRECTORIES(_sphw PRIVATE ${PROJECT_SOURCE_DIR}/fwk/_src
- ${PROJECT_SOURCE_DIR}/tbox/_src
- ${PYTHON_INCLUDE_PATH}
-)
-
-SWIG_LINK_LIBRARIES(sphw
- sph tbox fwk ${PYTHON_LIBRARIES}
-)
diff --git a/sph/_src/sphw.i b/sph/_src/sphw.i
deleted file mode 100644
index baa298c937fb3dd5c8fef0e812aea977ffc6d285..0000000000000000000000000000000000000000
--- a/sph/_src/sphw.i
+++ /dev/null
@@ -1,106 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-// SWIG input file of the 'sph' module
-
-%feature("autodoc","1");
-
-%module(docstring=
-"'sphw' module: projet MP 2016
-(c) ULg - A&M",
-directors="1",
-threads="1"
-) sphw
-%{
-
-#include <string>
-#include <sstream>
-#include "sph.h"
-
-#include "fwkw.h"
-#include "tboxw.h"
-
-#include "wKernel.h"
-#include "wDofs.h"
-#include "wParticle.h"
-#include "wEqState.h"
-#include "wTimeIntegration.h"
-#include "wSorter.h"
-#include "wProblem.h"
-#include "wEigenTest.h"
-#include <Eigen/Core>
-
-#include "wDisplayHook.h"
-
-%}
-
-%include "std_vector.i" // should be included first otherwise SWIGPY_SLICE_ARG is undefined
-
-%include "fwkw.swg"
-
-// ----------- MODULES UTILISES ------------
-%import "fwkw.i"
-%import "tboxw.i"
-
-// ----------- SPH CLASSES ---------------
-%include "sph.h"
-
-%shared_ptr(sph::EqState);
-%shared_ptr(sph::IdealGas);
-%shared_ptr(sph::QIncFluid);
-%shared_ptr(sph::Kernel);
-%shared_ptr(sph::CubicSplineKernel);
-%shared_ptr(sph::QuadraticKernel);
-%shared_ptr(sph::QuinticSplineKernel);
-%shared_ptr(sph::Sorter);
-%shared_ptr(sph::LouisSorter);
-%shared_ptr(sph::BruteForceSorter);
-%shared_ptr(sph::Problem);
-
-%feature("director") DisplayHook;
-%include "wDisplayHook.h"
-
-%feature("director:except") {
- if ($error != NULL) {
- std::cout << "[in director:except]\n";
- //throw Swig::DirectorMethodException();
- throw std::runtime_error("Director problem");
- }
-}
-
-
-%include "wKernel.h"
-%include "wDofs.h"
-%include "wParticle.h"
-%include "wEqState.h"
-%include "wTimeIntegration.h"
-%include "wSorter.h"
-%include "wProblem.h"
-%include "wEigenTest.h"
-
-%extend sph::Dofs {
- std::string __str__() {
- std::stringstream str; str << *self;
- return str.str();
- }
-}
-
-
-
-// Instantiate some std templates
-namespace std {
- %template(std_vector_Particlep) std::vector<sph::Particle*>;
-}
diff --git a/sph/broken/cpp/waterdrop.py b/sph/broken/cpp/waterdrop.py
deleted file mode 100644
index 8f6e318caf2a9d4471fc03bd393f48bb1d7a4e4b..0000000000000000000000000000000000000000
--- a/sph/broken/cpp/waterdrop.py
+++ /dev/null
@@ -1,48 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-from sph.helpers import *
-
-if __name__=="__main__":
-
- boxL = 2.
- Lfloor = 0.7
- Lwater = 0.5
- sep = 0.05
-
- kernel = Kernel('cubic', True) # 'cubic', 'quadratic' or 'quintic'
- law = EqState('liquid') # 'gas' or 'liquid'
-
- # parameters
- model = Model()
- model.kernel = kernel
- model.law = law
-
- model.h_0 = 0.06 # initial smoothing length [m]
- model.c_0 = 35.0 # initial speed of sound [m/s]
- model.rho_0 = 1000.0 # initial density [kg/m^3]
- model.dom_dim = boxL # domain size (cube)
- model.alpha = 0.5 # artificial viscosity factor 1
- model.beta = 0.0 # artificial viscosity factor 2
- model.maxTime = 1.0 # simulation time
- model.saveInt = 0.01 # save interval
-
- # mobile particles
- cube = Cube( o=(((boxL-Lwater)/2),((boxL-Lwater)/2), ((boxL)/2)+0.5), L=(Lwater,Lwater,Lwater), rho=model.rho_0, s=sep)
- model.addMobile(cube.generate())
-
- # fixed particles
- plane = Cube( o=(((boxL-Lfloor)/2),((boxL-Lfloor)/2), (boxL/2)), L=(Lfloor,Lfloor,sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- plane = Cube( o=(0,0,0), L=(boxL,boxL,sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
-
- # run SPH model
- print(model)
- model.run()
-
- # convert to VTK
- import sph.gui as gui
- gui.ToParaview(verb=False).convertall()
-
-
\ No newline at end of file
diff --git a/sph/genp.py b/sph/genp.py
deleted file mode 100644
index 79698d948d064c6e5d2cf0d92d919504bb5225c7..0000000000000000000000000000000000000000
--- a/sph/genp.py
+++ /dev/null
@@ -1,52 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-import math
-
-class Cube:
- """ a basic "cube" defined by its origin (o), size (L), density (rho) and distance between layers (s)
- """
- def __init__(self, o=(0.0, 0.0, 0.0), L=(1.0, 1.0, 1.0), s=0.1):
- self.ox = o[0]
- self.oy = o[1]
- self.oz = o[2]
- self.Lx = L[0]
- self.Ly = L[1]
- self.Lz = L[2]
- self.s = s
-
- def generate(self):
- parts = []
- ni = int(math.ceil((self.Lx/self.s)))+1
- dx = self.Lx/ni
- nj = int(math.ceil((self.Ly/self.s)))+1
- dy = self.Ly/nj
- nk = int(math.ceil((self.Lz/self.s)))+1
- dz = self.Lz/nk
- vol = (dx*dy*dz)
-
- for i in range(ni):
- x = self.ox+ i*dx
- for j in range(nj):
- y = self.oy+j*dy
- for k in range(nk):
- z = self.oz+k*dz
- p = (x, y, z, vol)
- parts.append(p)
- return parts
-
\ No newline at end of file
diff --git a/sph/louis/.gitignore b/sph/louis/.gitignore
deleted file mode 100644
index 8d83a3c883f32bd28d02ccae1837fdd5e626f1ea..0000000000000000000000000000000000000000
--- a/sph/louis/.gitignore
+++ /dev/null
@@ -1,17 +0,0 @@
-build/
-*.mod
-*.53
-*.x3d
-
-*.fp
-*.mp
-*.prm
-
-log.txt
-paths.txt
-*.out
-*.res
-*.vtu
-
-*.kdev4
-
diff --git a/sph/louis/CMakeLists.txt b/sph/louis/CMakeLists.txt
deleted file mode 100644
index 5482be3ce8bf58322e9d47145e59e3148dd6d840..0000000000000000000000000000000000000000
--- a/sph/louis/CMakeLists.txt
+++ /dev/null
@@ -1,41 +0,0 @@
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-PROJECT(LOUIS Fortran)
-CMAKE_MINIMUM_REQUIRED(VERSION 2.6)
-
-#SET(CMAKE_VERBOSE_MAKEFILE TRUE) # affiche les lignes de commande
-
-# configure Fortran compiler
-GET_FILENAME_COMPONENT(Fortran_COMPILER_NAME ${CMAKE_Fortran_COMPILER} NAME)
-
-IF (Fortran_COMPILER_NAME MATCHES "f95" OR Fortran_COMPILER_NAME MATCHES "gfortran")
- SET(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Wall -fopenmp")
-ENDIF()
-
-IF(Fortran_COMPILER_NAME MATCHES "ifort")
- IF(WIN32)
- SET(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} /Qopenmp" )
- SET(CMAKE_Fortran_FLAGS_RELEASE "/O3 /D NDEBUG")
- #SET(CMAKE_Fortran_FLAGS_DEBUG "/Od /debug:full /dbglibs")
- ENDIF()
-ENDIF()
-
-MESSAGE(STATUS "Fortran_COMPILER_NAME : " ${Fortran_COMPILER_NAME})
-MESSAGE(STATUS "CMAKE_Fortran_FLAGS : " ${CMAKE_Fortran_FLAGS})
-MESSAGE(STATUS "CMAKE_Fortran_FLAGS_RELEASE : " ${CMAKE_Fortran_FLAGS_RELEASE})
-MESSAGE(STATUS "CMAKE_Fortran_FLAGS_DEBUG : " ${CMAKE_Fortran_FLAGS_DEBUG})
-
-ADD_EXECUTABLE(louis src/SPH_module.f90 src/SPH_simulation.f90)
-
diff --git a/sph/louis/README.md b/sph/louis/README.md
deleted file mode 100644
index 6d50fb379e33589fcbf762807d7cb8202e9c9455..0000000000000000000000000000000000000000
--- a/sph/louis/README.md
+++ /dev/null
@@ -1,21 +0,0 @@
-# SPH Louis Goffin
-
-[Reference in ORBi](http://orbi.ulg.ac.be/handle/2268/156166)
-
-
-
-## Compilation (Windows or Linux)
-
-```bash
-build.py
-```
-requires [Python](https://www.python.org/), [CMake](https://cmake.org/) and a fortran compiler (either [Intel](https://software.intel.com/en-us/fortran-compilers) or [gfortran](https://gcc.gnu.org/fortran/))
-
-## Run
-
-```bash
-run.py -k 12 tests/waterdrop.py
-```
-Results go to `workspace/test_waterdrop`. Conversion to paraview requires [VTK](http://www.vtk.org/).
-
-Load `.vtu` files in [Paraview](http://www.paraview.org/).
diff --git a/sph/louis/build.py b/sph/louis/build.py
deleted file mode 100755
index bd218164247fe4fbc9f55e2cb129a6db01f20f22..0000000000000000000000000000000000000000
--- a/sph/louis/build.py
+++ /dev/null
@@ -1,38 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-
-
-import os, sys, subprocess, platform
-
-try:
- # create build dir
- os.chdir('..')
- if not os.path.isdir('louisB'): os.mkdir('louisB')
- os.chdir('louisB')
- # cmake
- if 'Windows' in platform.uname():
- subprocess.call(r'cmake -G "Visual Studio 11 Win64" ..\louis', shell=True)
- else:
- subprocess.call('cmake -G"Eclipse CDT4 - Unix Makefiles" -DCMAKE_ECLIPSE_VERSION=4.6 -DCMAKE_ECLIPSE_GENERATE_SOURCE_PROJECT=TRUE ../louis', shell=True)
- #subprocess.call('cmake -G"Eclipse CDT4 - Unix Makefiles" -DCMAKE_ECLIPSE_GENERATE_SOURCE_PROJECT=TRUE ../louis', shell=True)
- subprocess.call('cmake --build . --config Release', shell=True)
- os.chdir('../louis')
-except Exception as e:
- print(e)
- print("<press ENTER to quit>"); input()
diff --git a/sph/louis/devenv-vs2012.bat b/sph/louis/devenv-vs2012.bat
deleted file mode 100644
index 8a275be2c9d8ce7f04a4b7b11d7d486f9197ff94..0000000000000000000000000000000000000000
--- a/sph/louis/devenv-vs2012.bat
+++ /dev/null
@@ -1,8 +0,0 @@
-@echo off
-:: - cmake -G "Visual Studio 11 2012 Win64" ..\waves
-:: - cmake --build . --config Release
-
-:: MKL + VS en ligne de commande
-call "C:\Program Files (x86)\Intel\Composer XE\mkl\bin\intel64\mklvars_intel64.bat"
-call "C:\Program Files (x86)\Intel\Composer XE\tbb\bin\tbbvars.bat" intel64 vs2012
-%comspec% /K ""C:\Program Files (x86)\Microsoft Visual Studio 11.0\VC\vcvarsall.bat" amd64"
diff --git a/sph/louis/run.py b/sph/louis/run.py
deleted file mode 100755
index 8ee193c6e3cb10db9bbc8f09a0eb4a86ce71a0c3..0000000000000000000000000000000000000000
--- a/sph/louis/run.py
+++ /dev/null
@@ -1,80 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-# runs a test as if it was installed
-# - fixes the python path in a dev environment
-# - creates a workspace folder
-
-
-class DupStream:
- def __init__(self, stream1, stream2):
- self.stream1 = stream1
- self.stream2 = stream2
- def write(self, data):
- self.stream1.write(data)
- self.stream2.write(data)
- def flush(self):
- self.stream1.flush()
- self.stream2.flush()
-
-
-class Tee:
- def __init__(self, name):
- import sys
- self.file = open(name, 'w')
- self.stdoutbak = sys.stdout
- self.stderrbak = sys.stderr
- sys.stdout = DupStream(sys.stdout, self.file)
- sys.stderr = DupStream(sys.stderr, self.file)
- def __del__(self):
- import sys
- sys.stdout = self.stdoutbak
- sys.stderr = self.stderrbak
- self.file.close()
-
-
-if __name__=="__main__":
- import sys, os
- # adds "." to the pythonpath
- thisdir = os.path.split(__file__)[0]
- sys.path.append(thisdir)
- import sph.wutils as wu
-
- # parse args
- args = wu.parseargs()
-
- # run all tests sequentially
- for testname in args.file:
- testname = os.path.abspath(testname)
- if not os.path.isfile(testname):
- raise Exception("file not found: %s" % testname)
-
- wu.setupwdir(testname)
- __file__ = testname
-
- # split streams
- tee = Tee('stdout.txt')
-
- # start test
- import time, platform
- print('-'*80)
- print("starting test", testname)
- print("time:", time.strftime("%c"))
- print("hostname:", platform.node())
- print('-'*80)
- exec(open(testname, 'r', encoding='utf8').read())
diff --git a/sph/louis/screenshot.png b/sph/louis/screenshot.png
deleted file mode 100644
index 57acc618a70984dba4acc4b5238b07e7171ff0dc..0000000000000000000000000000000000000000
Binary files a/sph/louis/screenshot.png and /dev/null differ
diff --git a/sph/louis/sph/__init__.py b/sph/louis/sph/__init__.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/sph/louis/sph/gui.py b/sph/louis/sph/gui.py
deleted file mode 100644
index a3f6702c360315dc42c52252bff814f7b86c21c0..0000000000000000000000000000000000000000
--- a/sph/louis/sph/gui.py
+++ /dev/null
@@ -1,146 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-# converts .res files to Paraview (VTK)
-
-
-
-
-
-import vtk
-version=vtk.vtkVersion().GetVTKMajorVersion()
-import glob
-
-class ToParaview:
- def __init__(self, verb = False):
- self.verb = verb
-
- def convertall(self, pattern='*MP*.res'):
- print("converting grid to VTK")
- self.convertGrid()
- print("converting %s to VTK" % pattern)
- for f in glob.glob(pattern):
- self.convertParts(f)
-
- def convertGrid(self, fname='grid.out'):
- if self.verb: print('converting', fname)
-
- # read file
- file = open(fname)
- line = file.readline()
- nx, dx = line.strip().split()
- nx = int(nx)
- dx = float(dx)
- #print 'nx=%d, dx=%f' % (nx,dx)
- file.close()
-
- # build a sgrid
- grid = vtk.vtkStructuredGrid()
- points = vtk.vtkPoints()
- grid.SetPoints(points)
-
- for k in range(nx+1):
- for j in range(nx+1):
- for i in range(nx+1):
- points.InsertNextPoint(i*dx,j*dx,k*dx)
- grid.SetDimensions(nx+1,nx+1,nx+1)
-
- writer = vtk.vtkXMLStructuredGridWriter()
- compressor = vtk.vtkZLibDataCompressor()
- writer.SetCompressor(compressor)
- writer.SetDataModeToBinary()
- if version>5:
- writer.SetInputData(grid)
- else:
- writer.SetInput(grid)
- writer.SetFileName(fname.replace('.out','.vts'))
- writer.Write()
-
- def convertParts(self, fname):
-
- #if self.verb:
- print('converting', fname)
- file = open(fname)
-
- ugrid = vtk.vtkUnstructuredGrid()
- points = vtk.vtkPoints()
- ugrid.SetPoints(points)
-
- codes = {'p': (1, "Pressure"),
- 'rho': (1, "Mass density"),
- 'v' : (3, "Velocity"),
- 'm': (1, "Mass"),
- 'c': (1, "Speed of sound"),
- 'h': (1, "Smoothing length"),
- 'mu': (1, "max(mu_ab)"),
- 'nv': (1, "Nb of neighbours") }
- results = {}
- for c in codes:
- scalars = vtk.vtkFloatArray()
- ncomp, fullname = codes[c]
- scalars.SetNumberOfComponents(ncomp)
- scalars.SetName(fullname)
- ugrid.GetPointData().AddArray(scalars)
- results[c] = scalars
-
- i=0
- for line in file:
- #x,y,z, vx,vy,vz, m, p = map(float, line.strip().split())
- try:
- x,y,z, vx,vy,vz, rho, p, m, c, h, mu, nv = list(map(float, line.strip().split()))
- except:
- print("**ERROR while reading file %s!\n\tline=\"%s\"" % (fname,line))
- break
- points.InsertPoint(i, x,y,z)
- results['v'].InsertNextTuple3(vx,vy,vz)
- results['p'].InsertNextValue(p)
- results['rho'].InsertNextValue(rho)
- results['m'].InsertNextValue(m)
- results['c'].InsertNextValue(c)
- results['h'].InsertNextValue(h)
- results['mu'].InsertNextValue(mu)
- results['nv'].InsertNextValue(nv)
- i+=1
- file.close()
-
- ntot=i
- if self.verb: print("\t%d lines read. Converting grid to VTK..." % ntot)
-
- for i in range(ntot):
- vertex = vtk.vtkVertex()
- ids = vertex.GetPointIds()
- ids.SetId(0, i)
- ugrid.InsertNextCell(vertex.GetCellType(), ids)
-
- if self.verb:
- print("\t...", ugrid.GetNumberOfPoints(), 'points and', \
- ugrid.GetNumberOfCells(), 'cells converted')
-
- writer = vtk.vtkXMLUnstructuredGridWriter()
- compressor = vtk.vtkZLibDataCompressor()
- writer.SetCompressor(compressor)
- writer.SetDataModeToBinary()
- if version>5:
- writer.SetInputData(ugrid)
- else:
- writer.SetInput(ugrid)
- writer.SetFileName(fname.replace('.res','.vtu'))
- writer.Write()
-
-if __name__=="__main__":
- ToParaview(verb=True).convertall()
diff --git a/sph/louis/sph/helpers.py b/sph/louis/sph/helpers.py
deleted file mode 100644
index c05b7713da85d180265ebbb9b7e42539341c1c9e..0000000000000000000000000000000000000000
--- a/sph/louis/sph/helpers.py
+++ /dev/null
@@ -1,253 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-# Python helper classes for Louis' code
-
-import math, subprocess, os, platform, sys, glob
-import sph.wutils as wu
-
-class Kernel:
- names = { 'cubic':1, 'quadratic':2, 'quintic':3 }
- def __init__(self, name='cubic', corr=True):
- if name in Kernel.names:
- self.name = name
- else:
- raise Exception('unknown kernel kind "%s"' % name)
- self.corrected = corr
-
- def __str__(self):
- return '%s kernel (%d)' % (self.name, self.val())
- def val(self):
- return Kernel.names[self.name]
-
-class EqState:
- names = { 'gas':1, 'liquid':2 }
- def __init__(self, name='liquid'):
- if name in EqState.names:
- self.name = name
- else:
- raise Exception('unknown equ of state "%s"' % name)
- self.gamma = 7
- self.molMass = 28.97e-3
- def __str__(self):
- return '%s (%d)' % (self.name, self.val())
- def val(self):
- return EqState.names[self.name]
-
-
-class Model:
- def __init__(self):
- # default values of the parameters
- self.h_0 = 0.1 # 3: [double] initial smoothing length
- self.c_0 = 1480.0 # 4: [double] initial speed of sound [m/s]
- self.rho_0 = 1000.0 # 5: [double] initial density [kg/m^3]
- self.dom_dim = 0.0 # 6: [double] domain size (cube)
- self.kernel = Kernel() # 7: [integer] (1:'cubic'/2:'quadratic'/3:'quintic')
- self.alpha = 0.5 # 8: [double] artificial viscosity factor 1
- self.beta = 0.0 # 9: [double] artificial viscosity factor 2
- self.law = EqState() # 10: [integer] (1:'gas'/2:'fluid')
- self.maxTime = 1.0 # 14: [double] simulation time [s]
- self.saveInt = 0.01 # 15: [double] save interval [s]
-
- # sets of particles
- self.fparts = []
- self.mparts = []
-
- def writeprm(self):
- """ writes parameters (input.prm)
- """
- file = open('input.prm','w')
- file.write('%d\n' % len(self.fparts) )
- file.write('%d\n' % len(self.mparts) )
- file.write('%e\n' % self.h_0)
- file.write('%e\n' % self.c_0)
- file.write('%e\n' % self.rho_0)
- file.write('%e\n' % self.dom_dim)
- file.write('%d\n' % self.kernel.val())
- file.write('%e\n' % self.alpha)
- file.write('%e\n' % self.beta)
- file.write('%d\n' % self.law.val())
- file.write('%d\n' % self.law.gamma)
- file.write('%e\n' % self.law.molMass)
- file.write('%d\n' % self.kernel.corrected)
- file.write('%e\n' % self.maxTime)
- file.write('%e\n' % self.saveInt)
- file.close()
-
- def addFixed(self, parts):
- """ adds fixed particles to the model
- """
- self.fparts += parts
-
- def addMobile(self, parts):
- """ adds mobile particles to the model
- """
- self.mparts += parts
-
- def clean(self):
- for p in ['res*.*', 'input.*', 'grid.*']:
- for f in glob.glob(p):
- print('rm %s' % f)
- os.remove(f)
-
- def run(self):
-
- # clean prev results
- self.clean()
-
- # write parameters
- self.writeprm()
-
- # input.mp - contains mobile particles
- file = open('input.mp','w')
- for p in self.mparts:
- p.write(file)
- file.close()
-
- # input.fp - contains fixed particles
- file = open('input.fp','w')
- for p in self.fparts:
- p.write(file)
- file.close()
-
- # paths.txt - contains path to the input files
- file = open('paths.txt','w')
- file.write('input.prm\n')
- file.write('input.fp\n')
- file.write('input.mp\n')
- file.close()
-
- # set nb of threads
- args = wu.parseargs()
- os.environ['OMP_NUM_THREADS'] = str(args.k)
-
- # run program
- exename = self.getexe()
- print("running %s using %s threads" % (exename, os.environ['OMP_NUM_THREADS']))
-
- #iop = subprocess.call(self.getexe(), shell=True, stdout=sys.stdout, stderr=sys.stderr)
- #print 'louis returned iop=%d', iop
-
- #http://stackoverflow.com/questions/2715847/python-read-streaming-input-from-subprocess-communicate/17698359#17698359
-
- #file=open('pipo.txt','w')
- try:
- proc = subprocess.Popen(exename, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, bufsize=1)
- with proc.stdout:
- for line in iter(proc.stdout.readline, b''):
- line = line.rstrip('\n').rstrip('\r')
- print('[F]%s' % line)
- #sys.stdout.write(line)
- #file.write("line=\"%s\"\n" % line)
- proc.wait()
- except KeyboardInterrupt:
- print('Ignoring CTRL-C')
- pass
- #file.close()
-
- def getexe(self):
- """ looks for fortran exe
- """
- dir1=os.path.abspath(os.path.dirname(__file__)+os.sep+".."+os.sep+"..")+os.sep+"louisB"
- if 'Windows' in platform.uname():
- exename = os.path.join(dir1, "Release/louis.exe")
- else:
- exename = os.path.join(dir1, "louis")
- if not os.path.isfile(exename):
- raise Exception ("%s NOT found" %exename)
- return exename
-
- def __str__(self):
- txt = "SPH Model:\n"
- txt += "\tnb fixed particiles = %d\n" % len(self.fparts)
- txt += "\tnb mobile particules = %d\n" % len(self.mparts)
- txt += "\tinitial smoothing length = %f\n" % self.h_0
- txt += "\tinitial speed of sound [m/s] = %f\n" % self.c_0
- txt += "\tinitial density [kg/m^3] = %f\n" % self.rho_0
- txt += "\tdomain size (cube) = %f\n" % self.dom_dim
- txt += "\tkernel kind = %s\n" % self.kernel
- txt += "\tartificial viscosity factor 1 = %f\n" % self.alpha
- txt += "\tartificial viscosity factor 2 = %f\n" % self.beta
- txt += "\tequ of state (1:'gas'/2:'fluid') = %s\n" % self.law
- txt += "\tfluid prm = %d\n" % self.law.gamma
- txt += "\tgas prm [kg/mol] = %f\n" % self.law.molMass
- txt += "\tkernel correction (0:no 1:yes) = %s\n" % self.kernel.corrected
- txt += "\tsimulation time [s] = %f\n" % self.maxTime
- txt += "\tsave interval [s] = %f\n" % self.saveInt
- return txt
-
-class Particle:
- def __init__(self, x,y,z, vx,vy,vz, rho0, m0):
- self.x = x
- self.y = y
- self.z = z
- self.vx = vx
- self.vy = vz
- self.vz = vz
- self.rho0 = rho0
- self.m0 = m0
- def write(self, file):
- file.write('%e\t%e\t%e\t' % (self.x,self.y,self.z))
- file.write('%e\t%e\t%e\t' % (self.vx,self.vy,self.vz))
- file.write('%e\t%e\n' % (self.rho0, self.m0))
- def __str__(self):
- return "pos=(%f, %f, %f) v=(%f, %f, %f) rho=%f m=%f" % \
- (self.x, self.y, self.z, self.vx, self.vy, self.vz, self.rho0, self.m0)
-
-class Cube:
- """ a basic "cube" defined by its origin (o), size (L), density (rho) and distance between layers (s)
- """
- def __init__(self, o=(0.0, 0.0, 0.0), L=(1.0, 1.0, 1.0), rho=1.0, s=0.1):
- self.ox = o[0]
- self.oy = o[1]
- self.oz = o[2]
- self.Lx = L[0]
- self.Ly = L[1]
- self.Lz = L[2]
- self.rho = rho
- self.s = s
-
- def generate(self):
- parts = []
- ni = int(math.ceil((self.Lx/self.s)))+1
- dx = self.Lx/(ni-1)
- nj = int(math.ceil((self.Ly/self.s)))+1
- dy = self.Ly/(nj-1)
- nk = int(math.ceil((self.Lz/self.s)))+1
- dz = self.Lz/(nk-1)
- vx = vy = vz = 0.0
- m0 = (dx*dy*dz)*self.rho
- rho0 = self.rho
-
- for i in range(ni):
- x = self.ox+ i*dx
- for j in range(nj):
- y = self.oy+j*dy
- for k in range(nk):
- z = self.oz+k*dz
- p = Particle(x, y, z, vx, vy, vz, rho0, m0)
- parts.append(p)
- return parts
-
-
-
-
-
-
-
-
diff --git a/sph/louis/sph/wutils.py b/sph/louis/sph/wutils.py
deleted file mode 100644
index 2125cdc6c965e754c5a351c7e8318d8cfc7f04bd..0000000000000000000000000000000000000000
--- a/sph/louis/sph/wutils.py
+++ /dev/null
@@ -1,57 +0,0 @@
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-# Python utilities
-
-
-
-def setupwdir(testname):
- """
- creates a working folder for each test
- """
- import os, os.path
- #print "__file__=",__file__
- dir1=os.path.abspath(os.path.dirname(__file__)+os.sep+"..")+os.sep
- print("dir1=",dir1)
- print("testname=",testname)
- common = os.path.commonprefix( (testname, dir1) )
- #print "common=", common
- resdir = testname[len(common):].replace(os.sep,"_")
- resdir = os.path.splitext(resdir)[0] # remove ".py"
- #print "resdir=", resdir
- wdir=os.path.join('workspace', resdir)
- if not os.path.isdir(wdir):
- print("creating", wdir)
- os.makedirs(wdir)
- os.chdir(wdir)
-
-# ------------------------------------------------------------------------------
-
-def parseargs():
- """
- parses command line arguments
- """
- import argparse
- parser = argparse.ArgumentParser()
- parser.add_argument("-v", "--verb", help="increase output verbosity", action="count", default=0)
- parser.add_argument("--nogui", help="disable any graphical output",
- action="store_true")
- parser.add_argument("-k", help="nb of threads", type=int, default=1)
- #parser.add_argument("-p", help="misc parameters")
- parser.add_argument('file', nargs='*', help='python files')
- args = parser.parse_args()
- return args
diff --git a/sph/louis/src/SPH_module.f90 b/sph/louis/src/SPH_module.f90
deleted file mode 100644
index b6c1f72afcf04480964dfa2696cbdaf28ab55b85..0000000000000000000000000000000000000000
--- a/sph/louis/src/SPH_module.f90
+++ /dev/null
@@ -1,1241 +0,0 @@
-!>SPH_module
-!! @n Contains all the classes necessary to run a SPH simulation.
-!! @brief Group of classes (types) and procedures definitions.
-!! @author Louis Goffin
-!! @date 2013-05-26
-!! @version 1.0.0
-
-module SPH_module
-
- implicit none
-
- integer, parameter :: DP = KIND(1.0D0) !< double precision
- !integer, parameter :: DP = KIND(1.0) !< single precision
- real(DP), parameter :: pi = 3.141592653589793238462643383279502884197_dp !< \f$ \pi \f$
-
-
- !> Minimal link class
- !! @n This class is a minimal pointer, i.e. the object is composed of only a pointer.
- !! No other information is included. The pointers point toward particles
- type min_link
- class(fixed_particle), pointer :: ptr => null() !< pointer toward a particle
- end type min_link
-
-
- !> Link class
- !! @n This class contains a pointer that points toward an object
- !! and the distance between 2 particles. This class is used to
- !! build vectors of pointers toward objects.
-
- type link
- class(fixed_particle), pointer :: ptr => null() !< pointer toward a particle
- real(DP) :: r = 0.0_dp !< distance between neighbours
- end type link
-
-
- !> List class
- !! @n This class is a list that contains pointers toward objects (+distance).
- !! The only problem of this list class is that only link objects (ptr+r) can be added
-
- type list
- integer :: nbr = 0 !< number of elements in the list
- integer :: max_nbr = 0 !< maxnumber of elements in the list (size of the array)
- integer :: incr = 35 !< increment: number of spaces to add when the list is full
- type(link), dimension(:), allocatable :: lst !< list containing elements
-
- contains
- procedure :: initList => list_initList
- procedure :: addElement => list_addElement
- procedure :: resetList => list_resetList
- end type list
-
-
- !> Fixed particle class
- !! @n This class contains a certain number of parameters describing
- !! the state of a fixed particle (boundary particle).It also
- !! includes the needed procedures to calculate the continuity
- !! and some other equations.
-
- type fixed_particle
- real(DP), dimension(1:3, 1:3) :: coord !< 3x3 array containing the coordinates of a particle.
- !! column 1 = currentTime
- !! column 2 = RKstep
- !! column 3 = nextTime
- real(DP), dimension(1:3, 1:3) :: speed !< 3x3 array containing the velocity of a particle.
- !! column 1 = currentTime
- !! column 2 = RKstep
- !! column 3 = nextTime
- real(DP), dimension(1:3) :: rho !< 3x1 array containing the density of a particle.
- !! element 1 = currentTime
- !! element 2 = RKstep
- !! element 3 = nextTime
- real(DP) :: m !< mass of the particle
- real(DP), dimension(1:3) :: p !< pressure of the particle
- !! element 1 = currentTime
- !! element 2 = RKstep
- !! element 3 = nextTime
- real(DP), dimension(1:3) :: c !< speed of sound of a particle
- !! element 1 = currentTime
- !! element 2 = RKstep
- !! element 3 = nextTime
- real(DP) :: h !< smoothing length
- type(list) :: neighbours !< list of neighbours
- integer :: numOfNeighbours !< number of neighbours
- real(DP), dimension(1:3, 1:150) :: vec_gradW !< array that contains the gradient for every
- !! neighbours; initially set to 150 elements to
- !! increase the computational efficiency
- real(DP), dimension(1:3, 1:150):: vec_gradW_mod !< corrected vec_gradW if asked
- class(particle_manager), pointer :: manager !< pointer toward the object particle_manager
- real(DP) :: max_mu_ab !< maximum mu_ab of a particle (used for the timestep calculation)
-
- contains
- procedure :: save2disk => particle_save2disk
- procedure :: loadfromdisk => particle_loadfromdisk
- procedure :: getNeighbours
- procedure :: calcPressure
- procedure :: calcCelerity
- procedure :: gradW
- procedure :: kernel_corr
- procedure :: varUpdate => varUpdate_fixed
- end type fixed_particle
-
-
- !> Mobile_particle class
- !! @n This is an extension of the fixed_particle class.
- !! The procedure varUpdtate is overwritten to include the update of u and x.
-
- type, extends(fixed_particle) :: mobile_particle
- contains
- procedure :: varUpdate => varUpdate_mobile
- procedure :: ArtificialViscosity
- end type mobile_particle
-
-
- !> Particle_sort class
- !! @n This class is able to sort the particles. A grid is generated
- !! and the particles are sorted in each cell.
- type particle_sort
- real(DP) :: h_max !< maximum smoothing length
- real(DP) :: cellSize !< length of a side of a cube
- integer :: nCells = 0 !< number of cells in the domain
- integer :: nCellsSide !< number of cells on a row
- logical :: init !< true if the cells must be initialised
- type(list), dimension(:), allocatable :: storage !< vector of lists that contain
- !! the particles in a cell
- class(particle_manager), pointer :: manager !< pointer toward the object particle_manager
-
- contains
- procedure :: get_h_max
- procedure :: setCells
- procedure :: particlesSort
- end type particle_sort
-
- !> "enums"
-
- integer, parameter :: K_CUBIC_SPLINE = 1
- integer, parameter :: K_QUADRATIC = 2
- integer, parameter :: K_QUINTIC_SPLINE = 3
-
- integer, parameter :: KCORR_OFF = 0
- integer, parameter :: KCORR_ON = 1
-
- integer, parameter :: LAW_IDEAL_GAS = 1
- integer, parameter :: LAW_QINC_FLUID = 2
-
-
- !> Particle_manager class
- !! @n This class is used to manage all the particles,
- !! i.e. it contains a reference to every particles,
- !! it contains a number of parameters useful for the problem
- !! (variable smoothing length or not, ...), it has a solver, etc.
-
- type particle_manager
- type(particle_sort) :: sorting !< sorting machine
- type(min_link), dimension(:), allocatable :: part !< array of pointers toward particles
-
- integer :: numFP !< number of fixed particles
- integer :: numMP !< number of mobile particles
- integer :: numPart !< number of particles (FP+MP)
- integer :: kernelKind !< kind of kernel
- integer :: kappa !< kappa linked to the eqn state
- real(DP) :: alpha !< weighting factor in the artificial viscosity formulation
- real(DP) :: beta !< weighting factor in the artificial viscosity formulation
- integer :: eqnState !< equation of state
- !! 1 = ideal gas law
- !! 2 = quasi-incompressible fluid
- integer :: state_gamma !< power in eqn State 2.
- !! often taken around 7
- real(DP) :: molMass !< Molar mass of the fluid for the prefect gas law
-
- integer :: kernelCorrection !< correction of the kernel
- !! 0 = no correction
- !! 1 = correction enabled
-
- real(DP) :: maxTime !< simulation time in seconds
- real(DP) :: saveInt !< saving interval
- real(DP) :: h_0 !< initial smoothing length
- real(DP) :: rho_0 !< density of the fluid at free surface
- real(DP) :: c_0 !< speed of sound in normal conditions
- real(DP) :: timeStep !< timestep (not constant)
- real(DP) :: currentTime !< current time
- integer :: RKstep !< used to know in which RK iteration we are [RB] (1 or 2)
-
- real(DP) :: dom_dim !< length of a side of the domain
- !! (exterior particle to exterior particle).
- !! the domain is assumed to be cubic
-
- contains
- procedure :: readPRM
- procedure :: initialisation
- procedure :: solver
- procedure :: timeStepUpdate
- procedure :: slUpdate
- procedure :: savePartSet
- end type particle_manager
-
-
- contains
-
- !> calculates the distance between two particles
- function eval_r(xyz, xyz2)
- real(DP), dimension(3) :: xyz
- real(DP), dimension(3) :: xyz2
- real(DP) :: eval_r
-
- eval_r = sqrt(sum( (xyz(:)-xyz2(:))*(xyz(:)-xyz2(:)) ))
- end function eval_r
-
- ! ------------------------------------------------------------------------
- ! list
- ! ------------------------------------------------------------------------
-
- !> list/initList : initialise a list when it is first created
- subroutine list_initList(this)
- class(list) :: this
-
- if(this%incr < 1) then ! [RB] me semble inutile (incr=35 et ne varie pas)
- this%incr = 1
- end if
- this%max_nbr = this%incr
- allocate( this%lst(1:this%max_nbr) )
- end subroutine list_initList
-
-
- !> list/addElement : adds an element to a list.
- !! If the list is full, then it increases the size of the list.
- !! @param ptr : pointer of type link to be added in the list
- !! @param r : distance between two particles
-
- subroutine list_addElement(this, ptr, r)
- class(list) :: this
-
- class(fixed_particle), pointer :: ptr
- real(DP) :: r
- type(link), dimension(:), allocatable :: temp_lst !< temporary list used when it is necessary
- !! to increase the size of the existing list
-
- ! if the list is empty, it must be initialised
- if(this%max_nbr == 0) then
- call this%initList
-
- ! if the list is full it must be resized
- elseif(this%nbr == this%max_nbr) then
- allocate(temp_lst(1:this%max_nbr+this%incr))
- temp_lst(1:this%max_nbr) = this%lst(1:this%max_nbr)
- call move_alloc(temp_lst, this%lst) ! [RB] c'est une fct fortran
- this%max_nbr = this%max_nbr + this%incr
- end if
- this%nbr = this%nbr + 1
- this%lst(this%nbr)%ptr => ptr
- this%lst(this%nbr)%r = r
- end subroutine list_addElement
-
-
- !> list/resetList : resets the list.
- !! It only sets the number of elements to 0
- !! but it keeps the maximal size to max_nbr
-
- subroutine list_resetList(this)
- class(list) this
- this%nbr = 0
- end subroutine list_resetList
-
- ! ------------------------------------------------------------------------
- ! fixed_particle
- ! ------------------------------------------------------------------------
-
- subroutine particle_save2disk(this, ufile)
- class(fixed_particle), target :: this
- integer, intent(in) :: ufile !< unit
-
- write(unit=ufile, fmt="(E15.7,E15.7,E15.7,E15.7,E15.7,E15.7,E15.7,E15.7,E15.7,E15.7,E15.7,E15.7, I8)") &
- this%coord(1, 1), &
- this%coord(2, 1), &
- this%coord(3, 1), &
- this%speed(1, 1), &
- this%speed(2, 1), &
- this%speed(3, 1), &
- this%rho(1), &
- this%p(1), &
- this%m, &
- this%c(1), &
- this%h, &
- this%max_mu_ab, &
- this%numOfNeighbours
- end subroutine particle_save2disk
-
-
- subroutine particle_loadfromdisk(this, ufile, h_0)
- class(fixed_particle), target :: this
- integer, intent(in) :: ufile !< unit
- real(DP) :: x, y, z, u_x, u_y, u_z, rho, m
- real(DP) :: h_0
-
- read(ufile, *) x, y, z, u_x, u_y, u_z, rho, m
- this%coord = 0.d0
- this%coord(1, 1) = x
- this%coord(2, 1) = y
- this%coord(3, 1) = z
- this%speed = 0.d0
- this%speed(1, 1) = u_x
- this%speed(2, 1) = u_y
- this%speed(3, 1) = u_z
- this%rho = 0.d0
- this%rho(1) = rho
- this%m = m
- this%h = h_0
- this%p = 0.d0
- this%p(1) = this%calcPressure(this%rho(1))
- this%c = 0.d0
- this%c(1) = this%calcCelerity(this%rho(1))
- end subroutine particle_loadfromdisk
-
- !> fixed_particle/getNeighbours is a routine that questions the "particle_sort" object
- !! to get the particles in the neighbouring cells. The distance is
- !! calculated and the neighbours selected.
-
- subroutine getNeighbours(this)
- class(fixed_particle), target :: this
-
- real(DP), dimension(:), pointer, contiguous :: xyz !< position of the particle
- integer :: xCell, yCell, zCell !< number of the cell according to x, y and z
- integer :: nCellsSide !< number of cells on a row of the domain
- integer :: i, j, k !< loop counters
- integer, dimension(1:27) :: cellsToCheck !< number of the cells to check for the neighbours
- real(DP), dimension(:), pointer, contiguous :: neighXYZ !< coordinates of a neighbour
-
- real(DP) :: r !< distance between two particles
- class(fixed_particle), pointer :: cur_ptr !< current pointer toward a particle
- type(link), pointer :: cur_neigh !< current list of neighbours
- type(particle_sort), pointer :: srt !< pointer toward the sorting machine
- type(list), pointer :: storage !< pointer toward the storage
-
- integer :: cur_RKstep !< RKstep
-
- ! [RB] --
- integer :: twice
-
- ! pointer initialisation
- srt => this%manager%sorting ! [RB] chaque particule pointe vers le manager
- cur_RKstep = this%manager%RKstep
-
- ! coordinates of the particle
- xyz => this%coord(:, cur_RKstep)
-
- if(cur_RKstep == 1) then
-
- ! calculates the number of the cell in which the particle is
- nCellsSide = srt%nCellsSide
-
- xCell = nint( (xyz(1)-mod(xyz(1), srt%cellSize)) / srt%cellSize ) + 1
- yCell = nint( (xyz(2)-mod(xyz(2), srt%cellSize)) / srt%cellSize ) + 1
- zCell = nint( (xyz(3)-mod(xyz(3), srt%cellSize)) / srt%cellSize ) + 1
-
- if(xCell<1) then
- xCell = 1
- end if
- if(xCell>nCellsSide) then
- xCell = nCellsSide
- end if
- if(yCell<1) then
- yCell = 1
- end if
- if(yCell>nCellsSide) then
- yCell = nCellsSide
- end if
- if(zCell<1) then
- zCell = 1
- end if
- if(zCell>nCellsSide) then
- zCell = nCellsSide
- end if
-
- ! calculates the number of the neighbouring cells
- do i = -1, 1
- do j = -1, 1
- do k = -1, 1
- if((xCell+i>0).and.(yCell+j>0).and.(zCell+k>0).and. &
- (xCell+i <= nCellsSide).and. &
- (yCell+j <= nCellsSide).and. &
- (zCell+k <= nCellsSide) ) then
- cellsToCheck( (i+1)*9 + (j+1)*3 + (k+2) ) = &
- (xCell+i-1)*nCellsSide**2 + (yCell+j-1)*nCellsSide + (zCell+k)
- else
- cellsToCheck((i+1)*9 + (j+1)*3 + (k+2)) = 0
- end if
- end do
- end do
- end do
- end if
-
-
- ! stores the neighbours of the particle in the neighbours list.
- ! First, the list is reset, then the neighbouring cells are scanned.
- ! In each cell, the distance between the two particles is calculated.
- ! If it is lower than the support domain and it is not the particle we
- ! are working with (r>0 but here r>1E-12 for numerical errors), an element
- ! link(ptr+r) is added to the neighbours list.
- ! For the second RK step, only the distances r are recalculated. It is assumed that
- ! the neighbours remain the same between 2 RK step.
-
- if(cur_RKstep == 1) then
- call this%neighbours%resetList()
- twice = 0 ! [RB]
- do i = 1, 27
- if(cellsToCheck(i)>0) then
- storage => srt%storage(cellsToCheck(i))
- do j = 1, srt%storage(cellsToCheck(i))%nbr
- cur_ptr => storage%lst(j)%ptr
- neighXYZ => cur_ptr%coord(:, cur_RKstep)
- r = eval_r(xyz, neighXYZ)
- if(r<= this%manager%kappa*this%h) then
- if(r>1E-12) then ! [RB] why not cur_ptr /= this?
- call this%neighbours%addElement(cur_ptr, r)
- else
- twice = twice + 1 ! [RB]
- end if
- end if
- end do
- end if
- end do
- ! [RB] safeguard
- if (twice.ne. 1) then
- print *, 'safeguard activated!'
- print *, ' one particle has been taken into account', twice, ' times'
- print *, ' xCell =', xCell
- print *, ' yCell =', yCell
- print *, ' zCell =', zCell
- print *, ' cellsToCheck =', cellsToCheck
- stop
- end if
- this%numOfNeighbours = this%neighbours%nbr
- else
- ! RK step2 - same neighbours and r is updated
- do i = 1, this%numOfNeighbours
- cur_neigh => this%neighbours%lst(i)
- neighXYZ => cur_neigh%ptr%coord(:, cur_RKstep)
- cur_neigh%r = eval_r(xyz, neighXYZ)
- end do
- end if
- end subroutine getNeighbours
-
-
- !> fixed_particle/calcPressure is a function that calculates the pressure according
- !! to the equation of state chosen.
- !! @param rho : actual density
-
- function calcPressure(this, rho)
- class(fixed_particle) this
-
- real(DP) :: rho
- real(DP) :: calcPressure
- real(DP), parameter :: idealGasCst = 8.3144621d0
- real(DP) :: B
-
- select case(this%manager%eqnState)
- case ( LAW_IDEAL_GAS)
- calcPressure = (rho/this%manager%rho_0-1) &
- * idealGasCst*293.15d0/this%manager%molMass ! eq (3.24)
- case ( LAW_QINC_FLUID )
- B = this%manager%c_0**2 &
- * this%manager%rho_0 / this%manager%state_gamma ! eq (3.27)
- calcPressure = B*((rho/this%manager%rho_0)**this%manager%state_gamma - 1) ! eq (3.25)
- case default
- print *, 'Bad Equ of state (1,2)'
- stop
- end select
- end function calcPressure
-
-
- !> fixed_particle/calcCelerity : calculates the celerity according
- !! to the equation of state chosen.
- !! The equation used is @f[c = \sqrt{\frac{dp}{d\rho}}@f]
- !! @param rho : actual density
-
- function calcCelerity(this, rho)
- class(fixed_particle) :: this
-
- real(DP), intent(in) :: rho
- real(DP) :: calcCelerity
-
- select case(this%manager%eqnState)
- case( LAW_IDEAL_GAS )
- ! 1 = considering the ideal gas law at 20 degrees C
- calcCelerity = this%manager%c_0 ! eq (3.36)
- case ( LAW_QINC_FLUID )
- ! 2 = considering a quasi-incompressible fluid
- calcCelerity = this%manager%c_0 &
- * sqrt((rho/this%manager%rho_0)**(this%manager%state_gamma-1)) ! eq (3.37)
- case default
- print *, 'Bad Equ of state (1,2)'
- stop
- end select
- end function calcCelerity
-
-
- !> gradW : creates a vector that contains the values
- !! of the gradient for each neighbour.
-
- subroutine gradW(this)
- class(fixed_particle), target :: this
- integer :: i
- real(DP) :: alpha_d !< normalisation coefficient
- real(DP) :: r !< distance between a particle and a neighbour
- class(fixed_particle), pointer :: cur_neigh !< pointer toward a neighbour
- real(DP) :: cur_h !< value of h
- integer :: cur_RKstep !< pointer toward the current RK step
-
- if(this%numOfNeighbours>150) then
- print *, 'Error: Number of neighbours greater than expected (max 150 for vec_gradW): ', this%numOfNeighbours
- stop
- end if
-
- cur_h = this%h
- cur_RKstep = this%manager%RKstep
-
- select case(this%manager%kernelKind)
- case ( K_CUBIC_SPLINE)
-
- alpha_d = 3.d0/(2.d0*pi*cur_h**3) ! [RB] efficiency of x**3.0d0 vs x**3 vs x*x*x ??
- ! values of alpha_d in table 2.1 p 23
- do i = 1, this%numOfNeighbours
- r = this%neighbours%lst(i)%r ! [RB] a pointer is useless here!
- cur_neigh => this%neighbours%lst(i)%ptr
- if((r/cur_h>= 0.d0).and.(r/cur_h<1.d0)) then ! [RB] could "r/cur_h" be negative??
- this%vec_gradW(:, i) = alpha_d/cur_h &
- * (3.d0/2.d0*(r/cur_h)**2 - 2.d0*(r/cur_h) ) &
- * (this%coord(:, cur_RKstep) - cur_neigh%coord(:, cur_RKstep))/r ! eq (2.26)
- else if((r/cur_h>= 1.d0).and.(r/cur_h<2.d0)) then
- this%vec_gradW(:, i) = alpha_d/cur_h &
- * (-0.5d0*(2.d0-r/cur_h)**2) &
- * (this%coord(:, cur_RKstep) - cur_neigh%coord(:, cur_RKstep))/r
- else
- this%vec_gradW(:, i) = 0.d0
- end if
- end do
-
- case ( K_QUADRATIC )
-
- alpha_d = 5.d0/(4.d0*pi*cur_h**3)
-
- do i = 1, this%numOfNeighbours
- r = this%neighbours%lst(i)%r ! [RB] a pointer is useless here!
- cur_neigh => this%neighbours%lst(i)%ptr
-
- if((r/cur_h>= 0.d0).and.(r/cur_h<= 2.d0)) then ! [RB] could "r/cur_h" be negative??
- this%vec_gradW(:, i) = alpha_d/cur_h &
- * (3.d0/8.d0*r/cur_h-3.d0/4.d0) &
- * (this%coord(:, cur_RKstep)-cur_neigh%coord(:, cur_RKstep))/r ! eq (2.28)
- else
- this%vec_gradW(:, i) = 0.d0
- end if
- end do
-
- case ( K_QUINTIC_SPLINE)
-
- alpha_d = 3.d0/(359.d0*pi*cur_h**3)
-
- do i = 1, this%numOfNeighbours
- r = this%neighbours%lst(i)%r ! [RB] a pointer is useless here!
- cur_neigh => this%neighbours%lst(i)%ptr
- if((r/cur_h>= 0.d0).and.(r/cur_h<1.d0)) then
- this%vec_gradW(:, i) = alpha_d/cur_h &
- * ( -5.d0*(3.d0-r/cur_h)**4 &
- + 30.d0*(2.d0-r/cur_h)**4 &
- - 75.d0*(1.d0-r/cur_h)**4 ) &
- * (this%coord(:, cur_RKstep)-cur_neigh%coord(:, cur_RKstep))/r ! eq (2.32)
- elseif((r/cur_h>= 1.d0).and.(r/cur_h<2.d0)) then
- this%vec_gradW(:, i) = alpha_d/cur_h &
- * ( -5.d0*(3.d0-r/cur_h)**4 &
- +30.d0*(2.d0-r/cur_h)**4 ) &
- * (this%coord(:, cur_RKstep)-cur_neigh%coord(:, cur_RKstep))/r
- elseif((r/cur_h>= 2.d0).and.(r/cur_h<3.d0)) then
- this%vec_gradW(:, i) = alpha_d/cur_h &
- * (-5.d0*(3.d0-r/cur_h)**4) &
- * (this%coord(:, cur_RKstep)-cur_neigh%coord(:, cur_RKstep))/r
- else
- this%vec_gradW(:, i) = 0.d0
- end if
- end do
-
- case default
- print *, 'Bad value for kernel kind (1,2,3)'
- stop
- end select
-
- if(this%manager%kernelCorrection == KCORR_ON) then
- this%vec_gradW_mod = this%vec_gradW
- end if
- end subroutine gradW
-
-
- !> kernel_corr is a routine that takes into account the fact that the kernel may be truncated.
- !! It corrects the gradient of the kernel
-
- subroutine kernel_corr(this)
- class(fixed_particle) :: this
- real(DP), dimension(3, 3) :: M !< matrix used to correct the kernel gradient
- real(DP), dimension(3, 3) :: L !< inverse of the matrix used to correct the kernel gradient
- real(DP) :: detM !< determinant of M
- integer :: i !< loop counter
- real(DP) :: MDivRho !< m_b/rho_b
- class(fixed_particle), pointer :: cur_neigh !< pointer toward the current neighbour
- integer :: cur_RKstep !< current RK step
-
- cur_RKstep = this%manager%RKstep
- M(1, 1) = 0.d0
- M(2, 2) = 0.d0
- M(3, 3) = 0.d0
- M(1, 2) = 0.d0
- M(1, 3) = 0.d0
- M(2, 3) = 0.d0
-
- do i = 1, this%numOfNeighbours
- cur_neigh => this%neighbours%lst(i)%ptr
- MDivRho = cur_neigh%m/cur_neigh%rho(cur_RKstep)
- M(1, 1) = M(1, 1) + MDivRho * (cur_neigh%coord(1, cur_RKstep) - this%coord(1, cur_RKstep)) * this%vec_gradW(1, i)
- M(2, 2) = M(2, 2) + MDivRho * (cur_neigh%coord(2, cur_RKstep) - this%coord(2, cur_RKstep)) * this%vec_gradW(2, i)
- M(3, 3) = M(3, 3) + MDivRho * (cur_neigh%coord(3, cur_RKstep) - this%coord(3, cur_RKstep)) * this%vec_gradW(3, i)
- M(1, 2) = M(1, 2) + MDivRho * (cur_neigh%coord(1, cur_RKstep) - this%coord(1, cur_RKstep)) * this%vec_gradW(2, i)
- M(1, 3) = M(1, 3) + MDivRho * (cur_neigh%coord(1, cur_RKstep) - this%coord(1, cur_RKstep)) * this%vec_gradW(3, i)
- M(2, 3) = M(2, 3) + MDivRho * (cur_neigh%coord(2, cur_RKstep) - this%coord(2, cur_RKstep)) * this%vec_gradW(3, i)
- end do
- M(2, 1) = M(1, 2) !< M is symmetric
- M(3, 1) = M(1, 3) !< M is symmetric
- M(3, 2) = M(2, 3) !< M is symmetric
-
- detM = M(1, 1) * ( M(2, 2)*M(3, 3)-M(3, 2)*M(2, 3) ) &
- - M(1, 2) * ( M(2, 1)*M(3, 3)-M(3, 1)*M(2, 3) ) &
- + M(1, 3) * ( M(2, 1)*M(3, 2)-M(3, 1)*M(2, 2) )
- if(detM.eq.0.0) then
- print *, "detM==0!"
- stop
- end if
-
- L(1, 1) = M(2, 2)*M(3, 3) - M(3, 2)*M(2, 3)
- L(2, 2) = M(1, 1)*M(3, 3) - M(3, 1)*M(1, 3)
- L(3, 3) = M(1, 1)*M(2, 2) - M(2, 1)*M(1, 2)
- L(1, 2) = M(3, 1)*M(2, 3) - M(2, 1)*M(3, 3)
- L(2, 1) = L(1, 2) ! the inverse of a symmetric matrix is symmetric
- L(1, 3) = M(2, 1)*M(3, 2) - M(3, 1)*M(2, 2)
-
- L(3, 1) = L(1, 3) ! the inverse of a symmetric matrix is symmetric
- L(2, 3) = M(3, 1)*M(1, 2) - M(1, 1)*M(3, 2)
- L(3, 2) = L(2, 3) ! the inverse of a symmetric matrix is symmetric
- L = (1.d0/detM)*L
-
- do i = 1, this%numOfNeighbours
- !this%vec_gradW_mod(1, i) = L(1, 1)*this%vec_gradW(1, i) + L(1, 2)*this%vec_gradW(1, i) + L(1, 3)*this%vec_gradW(3, i)
- this%vec_gradW_mod(1, i) = L(1, 1)*this%vec_gradW(1, i) + L(1, 2)*this%vec_gradW(2, i) + L(1, 3)*this%vec_gradW(3, i) ! [RB] 1=>2 BUG?
- this%vec_gradW_mod(2, i) = L(2, 1)*this%vec_gradW(1, i) + L(2, 2)*this%vec_gradW(2, i) + L(2, 3)*this%vec_gradW(3, i)
- this%vec_gradW_mod(3, i) = L(3, 1)*this%vec_gradW(1, i) + L(3, 2)*this%vec_gradW(2, i) + L(3, 3)*this%vec_gradW(3, i)
- end do
-
- end subroutine kernel_corr
-
-
- !> fixed_particle/varUpdate_fixed : update the density of a fixed particle.
- !! The update of the velocity and the position are not performed.
- !! The integration scheme is a RK22 scheme.
-
- subroutine varUpdate_fixed(this)
- class(fixed_particle) :: this
- real(DP) :: Delta_rho !< \f$d\rho/dt\f$
- real(DP), dimension(1:3) :: u_ab !< relative velocity between the particle and a neighbour
- integer :: i !< loop counter
- integer :: cur_RKstep !< pointer toward the value of the current RK step
- class(fixed_particle), pointer :: cur_neigh !< current neighbour
- real(DP) :: dt
-
- Delta_rho = 0.d0
- cur_RKstep = this%manager%RKstep
-
- call this%getNeighbours()
- call this%gradW()
-
- do i = 1, this%numOfNeighbours
- cur_neigh => this%neighbours%lst(i)%ptr
- u_ab(:) = this%speed(:, cur_RKstep) - cur_neigh%speed(:, cur_RKstep)
- Delta_rho = Delta_rho + this%m * dot_product(u_ab, this%vec_gradW(:, i))
- end do
-
- dt = this%manager%timeStep
-
- if(cur_RKstep == 1) then !< 1st RK step
- this%rho(2) = this%rho(1) + Delta_rho * dt
- this%rho(3) = this%rho(1) + Delta_rho * dt/2.d0
- this%speed(:, 2) = this%speed(:, 1)
- this%coord(:, 2) = this%coord(:, 1)
- this%p(2) = this%calcPressure(this%rho(2))
- this%c(2) = this%calcCelerity(this%rho(2))
- else !< 2nd RK step
- this%rho(3) = this%rho(3) + Delta_rho * dt/2.d0
- this%speed(:, 3) = this%speed(:, 2)
- this%coord(:, 3) = this%coord(:, 2)
- this%p(3) = this%calcPressure(this%rho(3))
- this%c(3) = this%calcCelerity(this%rho(3))
- end if
- end subroutine varUpdate_fixed
-
- ! ------------------------------------------------------------------------
- ! mobile_particle
- ! ------------------------------------------------------------------------
-
- !> mobile_particle/varUpdate_mobile is a routine that updates the density, the velocity
- !! and the position of a mobile particle.
- !! The integration scheme is a RK22 scheme.
-
- subroutine varUpdate_mobile(this)
- class(mobile_particle) :: this
- real(DP) :: Delta_rho !< \f$d\rho/dt\f$
- real(DP), dimension(1:3) :: Delta_u !< \f$du/dt\f$
- real(DP), dimension(1:3) :: Delta_x !< \f$dx/dt\f$
- real(DP), dimension(1:3) :: u_ab !< relative velocity between the particle and a neighbour
- real(DP), dimension(1:3) :: F !< Volume forces
- real(DP) :: pi_ab !< Artificial viscosity
- integer :: i !< loopcounter
- integer :: cur_RKstep !< pointer toward the value of the current RK step
- class(fixed_particle), pointer :: cur_neigh !< current neighbour
- real(DP) :: dt
-
- Delta_rho = 0.d0
- Delta_u = (/ 0.d0, 0.d0, 0.d0 /)
- Delta_x = (/ 0.d0, 0.d0, 0.d0 /)
- F = (/ real(0, DP), real(0, DP), real(-9.81, DP) /)
-
- cur_RKstep = this%manager%RKstep
-
- call this%getNeighbours()
- call this%gradW()
-
- if(this%manager%kernelCorrection == KCORR_ON) then
- call this%kernel_corr()
- end if
-
- ! reset max_mu_ab
- if(cur_RKstep == 1) then
- this%max_mu_ab = 0.d0
- end if
-
-
-
- select case(this%manager%kernelCorrection)
- case( KCORR_ON )
- do i = 1, this%numOfNeighbours
- cur_neigh => this%neighbours%lst(i)%ptr
- u_ab(:) = this%speed(:, cur_RKstep) - cur_neigh%speed(:, cur_RKstep)
- pi_ab = this%ArtificialViscosity(cur_neigh, this%manager%alpha, this%manager%beta)
- Delta_rho = Delta_rho &
- + this%m * dot_product(u_ab, this%vec_gradW(:, i))
- Delta_u = Delta_u &
- + this%m * (cur_neigh%p(cur_RKstep) / cur_neigh%rho(cur_RKstep)**2 &
- + this%p(cur_RKstep) / this%rho(cur_RKstep)**2 &
- + pi_ab ) &
- * this%vec_gradW_mod(:, i) ! <= difference here
- end do
- case ( KCORR_OFF )
- do i = 1, this%numOfNeighbours
- cur_neigh => this%neighbours%lst(i)%ptr
- u_ab(:) = this%speed(:, cur_RKstep) - cur_neigh%speed(:, cur_RKstep)
- pi_ab = this%ArtificialViscosity(cur_neigh, this%manager%alpha, this%manager%beta)
- Delta_rho = Delta_rho &
- + this%m * dot_product(u_ab, this%vec_gradW(:, i))
- Delta_u = Delta_u &
- + this%m * (cur_neigh%p(cur_RKstep)/cur_neigh%rho(cur_RKstep)**2 &
- + this%p(cur_RKstep) / this%rho(cur_RKstep)**2 &
- + pi_ab ) &
- * this%vec_gradW(:, i) ! <= difference here
- end do
- case default
- print *,'Bad value of kernel correction'
- stop
- end select
-
- Delta_u = -Delta_u + F
- dt = this%manager%timeStep
-
- if(cur_RKstep == 1) then ! 1st RK step
- this%rho(2) = this%rho(1) + Delta_rho * dt
- this%rho(3) = this%rho(1) + Delta_rho * dt/2.d0
- this%speed(:, 2) = this%speed(:, 1) + Delta_u * dt
- this%speed(:, 3) = this%speed(:, 1) + Delta_u * dt/2.d0
- this%coord(:, 2) = this%coord(:, 1) + this%speed(:, 1) * dt
- this%coord(:, 3) = this%coord(:, 1) + this%speed(:, 1) * dt/2.d0
- this%p(2) = this%calcPressure(this%rho(2))
- this%c(2) = this%calcCelerity(this%rho(2))
- else ! 2nd RK step
- this%rho(3) = this%rho(3) + Delta_rho * dt/2.d0
- this%speed(:, 3) = this%speed(:, 3) + Delta_u * dt/2.d0
- this%coord(:, 3) = this%coord(:, 3) + this%speed(:, 2) * dt/2.d0
- this%p(3) = this%calcPressure(this%rho(3))
- this%c(3) = this%calcCelerity(this%rho(3))
- end if
- end subroutine varUpdate_mobile
-
-
- !> mobile_particle/ArtificialViscosity calculates the viscosity term
- !! in the momentum equation.
- !! @param neighObj : neighbouring object
- !! @param alpha : coefficicent in the artificial viscosity formulation
- !! @param beta : coefficicent in the artificial viscosity formulation
-
- function ArtificialViscosity(this, neighObj, alpha, beta)
- class(mobile_particle) :: this
- class(fixed_particle) :: neighObj
- real(DP) :: alpha
- real(DP) :: beta
- real(DP) :: ArtificialViscosity
- real(DP) :: mu_ab = 0.d0 !< this term represents a kind of viscosity
- real(DP), dimension(1:3) :: u_ab !< relative velocity of a in comparison with b
- real(DP), dimension(1:3) :: x_ab !< the distance between a and b
- real(DP) :: c_ab !< mean speed of sound
- real(DP) :: rho_ab !< mean density
-
- u_ab = this%speed(:, this%manager%RKstep) - neighObj%speed(:, this%manager%RKstep)
- x_ab = this%coord(:, this%manager%RKstep) - neighObj%coord(:, this%manager%RKstep)
-
- if(dot_product(u_ab, x_ab)<0.0d0) then
- ! mu_ab is calculated using \f[\mu_{ab} = \frac{h\vec{u}_{ab}\cdot\vec{x}_{ab}}{\vec{x}_{ab}^2+\eta^2}\f]
- mu_ab = this%h * dot_product(u_ab, x_ab) / ( dot_product(x_ab, x_ab) + 0.01d0*this%h**2.d0 )
- c_ab = 0.5d0 * (this%c(this%manager%RKstep) + neighObj%c(this%manager%RKstep))
- rho_ab = 0.5d0 * (this%rho(this%manager%RKstep) + neighObj%rho(this%manager%RKstep))
-
- ArtificialViscosity = (-alpha*c_ab*mu_ab + beta*mu_ab**2.d0)/rho_ab
- else
- ArtificialViscosity = 0.d0
- end if
-
- ! update of max_mu_ab for the calculation of the timestep
- if( (this%manager%RKstep == 1).and.(mu_ab>this%max_mu_ab) ) then
- this%max_mu_ab = mu_ab
- end if
- end function ArtificialViscosity
-
- ! ------------------------------------------------------------------------
- ! particle_manager
- ! ------------------------------------------------------------------------
-
- !> Reading and storing of the data in the parameter files
- subroutine readPRM(this, param_path)
- class(particle_manager), target :: this
- character(len=*) :: param_path
-
- open(unit = 1, file = trim(param_path))
- read(1, *) this%numFP
- read(1, *) this%numMP
- read(1, *) this%h_0
- read(1, *) this%c_0
- read(1, *) this%rho_0
- read(1, *) this%dom_dim
- read(1, *) this%kernelKind
- read(1, *) this%alpha
- read(1, *) this%beta
- read(1, *) this%eqnState
- read(1, *) this%state_gamma
- read(1, *) this%molMass
- read(1, *) this%kernelCorrection
- read(1, *) this%maxTime
- read(1, *) this%saveInt
- close(unit = 1)
- end subroutine readPRM
-
-
- !> particle_manager/initialisation : initialises the particle manager
- !! and all the particles. The routine takes the data from external files.
- !! The files paths are given in a file saved in the same directory
- !! as the program.
- !! The external files contains:
- !! the parameters,
- !! the fixed particles properties and the mobile particles.
-
- subroutine initialisation(this)
- class(particle_manager), target :: this
-
- character(250) :: param_path !< path of the parameters file
- character(250) :: fp_path !< path of the fixed particle(fp) file
- character(250) :: mp_path !< path of the mobile particle(mp) file
-
- integer :: i !< loop counter
- class(fixed_particle), pointer :: cur_ptr
-
- this%timeStep = 1.0d-15 !< initial time step
- this%currentTime = 0.d0 !< current time initialisation
- this%RKstep = 1 !< RK step counter initialisation
-
- ! Reading of the paths of the input files
- open(unit = 1, file = 'paths.txt')
- read(1, *) param_path
- read(1, *) fp_path
- read(1, *) mp_path
- close(unit = 1)
-
- call this%readPRM(param_path)
-
- ! allocation of the particles array
- this%numPart = this%numFP + this%numMP
- allocate(this%part(1:this%numPart))
-
- select case(this%kernelKind)
- case( K_CUBIC_SPLINE )
- this%kappa = 2
- case( K_QUADRATIC )
- this%kappa = 2
- case( K_QUINTIC_SPLINE )
- this%kappa = 3
- end select
-
- ! Reading and storing of the data for the fixed particles
- open(unit = 1, file = trim(fp_path))
- do i = 1, this%numFP
- allocate(fixed_particle :: this%part(i)%ptr)
- cur_ptr => this%part(i)%ptr
- cur_ptr%manager => this
- call cur_ptr%loadfromdisk(1, this%h_0)
- call cur_ptr%neighbours%initList()
- end do
- close(unit = 1)
-
- ! Reading and storing of the data for the mobile particles
- open(unit = 1, file = trim(mp_path))
- do i = 1, this%numMP
- allocate(mobile_particle :: this%part(this%numFP+i)%ptr)
- cur_ptr => this%part(this%numFP+i)%ptr
- cur_ptr%manager => this
- call cur_ptr%loadfromdisk(1, this%h_0)
- call cur_ptr%neighbours%initList()
- end do
- close(unit = 1)
-
- ! Particle sort
- this%sorting%manager => this
- this%sorting%init = .true.
-
- print *, 'Initialisation finished.'
- end subroutine initialisation
-
-
- !> save a particle set onto disk
-
- subroutine savePartSet(this, fname, ite, nbegin, nend)
- class(particle_manager) :: this
- character (len=*), intent(in) :: fname !< file name prefix
- integer, intent(in) :: ite
- integer, intent(in) :: nbegin, nend
-
- integer :: i, ios
- character (len=250) :: filename
-
- write(filename , fmt="(A,'_',I0.8,'.res')") fname,ite
- open (unit=1, file=filename, action="write", iostat=ios)
- if ( ios /= 0 ) then
- print "('Cannot open result file')"
- stop
- end if
- do i = nbegin, nend
- call this%part(i)%ptr%save2disk(1)
- end do
- close(unit=1)
- end subroutine savePartSet
-
-
- !> particle_manager/solver: solves the problem. It loops over time and uses
- !! a RK22 time integration scheme.
-
- subroutine solver(this)
- class(particle_manager) :: this
-
- integer :: i, j
- integer :: ite !< iteration counter
- logical :: to_save !< saving flag. If true a saving is done
- class(fixed_particle), pointer :: p
-
- ite = 0
-
- do while(this%currentTime <= this%maxTime)
-
- ! Time increment and saving status
- if((floor(this%currentTime/this%saveInt) /= &
- floor((this%currentTime+this%timeStep)/this%saveInt)) .or. ite == 0) then
- to_save = .true.
- end if
- this%currentTime = this%currentTime + this%timeStep
-
- ! Runge-Kutta loop
- do j = 1, 2
- this%RKstep = j
- !if (j.eq.1) then ! [RB]
- call this%sorting%particlesSort()
- !end if ! [RB]
- ! Loop over the particles
- !$OMP PARALLEL DO PRIVATE(i) SCHEDULE(DYNAMIC)
- do i = 1, this%numPart
- call this%part(i)%ptr%varUpdate()
- end do
- !$OMP END PARALLEL DO
- end do
-
- ! Update of the current time variables (currentTime = nextTime)
- do i = 1, this%numPart
- p => this%part(i)%ptr
- p%rho(1) = p%rho(3)
- p%p(1) = p%p(3)
- p%c(1) = p%c(3)
- p%speed(:, 1) = p%speed(:, 3)
- p%coord(:, 1) = p%coord(:, 3)
- end do
-
- ! Test for the data saving
- if(to_save) then
- call this%savePartSet('resMP', ite, this%numFP+1, this%numFP+this%numMP)
- call this%savePartSet('resFP', ite, 1, this%numFP)
-
- print *, 'Iteration nb ', ite
- print *, ' Time (s) = ', this%currentTime
- print *, ' Time step (s) = ', this%timeStep
- to_save = .false.
- end if
-
- call this%timeStepUpdate()
- call this%slUpdate()
-
- ite = ite + 1
- end do
-
- end subroutine solver
-
-
- !> particle_manager/timeStepUpdate: computes the next timestep using
- !! the properties of the particles.
-
- subroutine timeStepUpdate(this)
- class(particle_manager) :: this
- real(DP) :: dTf, dTftemp !< time step relative to the body forces
- real(DP) :: dTcv, dTcvtemp !< time step relative to the viscous forces and Courrant number
- integer :: i
- class(fixed_particle), pointer :: cur_ptr
-
- ! computes the timestep relative to the body forces
- dTf = sqrt( this%part(this%numFP+1)%ptr%h / 9.81d0 )
- do i = this%numFP+2, this%numPart
- cur_ptr => this%part(i)%ptr
- dTftemp = sqrt(cur_ptr%h / 9.81d0)
- if(dTftemp<dTf) then
- dTf = dTftemp
- end if
- end do
-
- ! computes the timestep relative to the CN and the viscous forces
- dTcv = this%part(this%numFP+1)%ptr%h &
- / (this%part(this%numFP+1)%ptr%c(1) &
- + 0.6d0 * (this%alpha * this%part(this%numFP+1)%ptr%c(1) &
- + this%beta * this%part(this%numFP+1)%ptr%max_mu_ab ) )
- do i = this%numFP+2, this%numPart
- cur_ptr => this%part(i)%ptr
- dTcvtemp = cur_ptr%h / &
- (cur_ptr%c(1) &
- + 0.6d0 * (this%alpha * cur_ptr%c(1) + this%beta * cur_ptr%max_mu_ab))
- if(dTcvtemp<dTcv) then
- dTcv = dTcvtemp
- end if
- end do
-
- ! computes the final timestep
- if(0.4d0*dTf > 0.25d0*dTcv) then
- this%timeStep = 0.25d0 * dTcv
- else
- this%timeStep = 0.4d0 * dTf
- end if
-
- !> possibility to change the timestep if we use the ideal gas law
- if(this%eqnState == 1) then
- this%timeStep = 5.d0 * this%timeStep
- end if
- end subroutine timeStepUpdate
-
-
- !> particle_manager/slUpdate is a routine that updates the smoothing length at each timestep.
- !! It is written to provide the same smoothing length for every particle.
-
- subroutine slUpdate(this)
- class(particle_manager) :: this
- real(DP) :: mean_rho !< mean value of the densities of the mobile particles
- real(DP) :: new_h !< new smoothing length
- integer :: i
- mean_rho = 0.d0
-
- ! calculation of the average density
- do i = 1, this%numPart
- mean_rho = mean_rho+this%part(i)%ptr%rho(1)
- end do
- mean_rho = mean_rho/this%numPart
-
- ! calculation of the new smoothing length
- new_h = this%h_0 * (this%rho_0/mean_rho)**(1.d0/3.d0)
-
- ! if the smoothing length is greater than 0.5 the size of a cell, it is limited
- if(new_h > 0.5d0 * this%sorting%cellSize) then
- new_h = 0.5d0 * this%sorting%cellSize
- print *, 'Warning: the smoothing has been limited'
- end if
-
- ! update of the smoothing length
- do i = 1, this%numPart
- this%part(i)%ptr%h = new_h
- end do
- end subroutine slUpdate
-
- ! ------------------------------------------------------------------------
- ! particle_sort
- ! ------------------------------------------------------------------------
-
- !> particle_sort/get_h_max: finds the largest smoothing length of the particles.
- !! This is useful when it is not constant over the particles.
-
- subroutine get_h_max(this)
- class(particle_sort) :: this
- integer :: i
-
- this%h_max = 0
- do i = 1, this%manager%numPart
- if(this%manager%part(i)%ptr%h > this%h_max) then
- this%h_max = this%manager%part(i)%ptr%h
- end if
- end do
-
- ! increase of h_max in order to have a security if h changes.
- ! This is done according to the equation of state used.
- if(this%manager%eqnState == LAW_IDEAL_GAS) then
- this%h_max = 1.1d0 * this%h_max
- else
- this%h_max = 1.02d0 * this%h_max
- end if
- end subroutine get_h_max
-
-
- !> particle_sort/setCells : sets the size of the cells in which the particles
- !! will be sorted. A cell must be cubic. The domain is assumed to be cubic.
- !! This routine also sets the number of cells and allocates the vector which contains
- !! the lists of particles.
- !! In ordrer to be as efficient as possible, the storage vector is not deallocated and
- !! reallocated at each iteration.
-
- subroutine setCells(this)
- class(particle_sort) :: this
- integer :: ios
-
- ! calculates the necessary number of cells on a side
- call this%get_h_max()
-
- this%nCellsSide = 0
- do while(this%manager%dom_dim / (this%nCellsSide+1) > this%manager%kappa * this%h_max)
- this%nCellsSide = this%nCellsSide + 1
- end do
-
- this%nCells = this%nCellsSide**3
-
- ! allocated with the necessary number of cells.
- allocate(this%storage(1:this%nCells))
- this%cellSize = this%manager%dom_dim / this%nCellsSide
-
- ! [RB] info
- print *, 'INFO particle_sort/setCells()'
- print *, ' .nCellsSide =', this%nCellsSide
- print *, ' .nCells =', this%nCells
- print *, ' .cellSize =', this%cellSize
-
- ! save info 2 disk
- open (unit=1, file="grid.out", action="write", iostat=ios)
- if ( ios /= 0 ) then
- print "('Cannot open grid file')"
- stop
- end if
- write(unit=1, fmt="(I8,E15.7)") &
- this%nCellsSide, &
- this%cellSize
- close(unit=1)
-
- end subroutine setCells
-
-
- !> particle_sort/particlesSort: sorts every particle in a cell.
- !! This will be useful in order to find the neighbours.
-
- subroutine particlesSort(this)
- class(particle_sort), target :: this
-
- integer :: i
- integer :: xCell, yCell, zCell !< number of the cell in the x, y and z direction
- integer :: part_pos !< absolute position of a particle
- real(DP), dimension(:), pointer, contiguous :: xyz !< position of a particle
- integer, pointer :: nCellsSide !< number of cells on a row
-
- class(fixed_particle), pointer :: prt
-
- !print *,"sorting particles"
-
- if(this%init) then
- call this%setCells()
- this%init = .false.
- end if
-
- ! the lists of every cell are reset
- do i = 1, this%nCells
- call this%storage(i)%resetList()
- end do
-
- nCellsSide => this%nCellsSide
- do i = 1, this%manager%numPart
- prt => this%manager%part(i)%ptr
- xyz => prt%coord(:, this%manager%RKstep)
-
- xCell = nint( (xyz(1)-mod(xyz(1), this%cellSize)) / this%cellSize ) + 1
- yCell = nint( (xyz(2)-mod(xyz(2), this%cellSize)) / this%cellSize ) + 1
- zCell = nint( (xyz(3)-mod(xyz(3), this%cellSize)) / this%cellSize ) + 1
-
- if(xCell<1) then
- xCell = 1
- end if
- if(xCell>nCellsSide) then
- xCell = nCellsSide
- end if
- if(yCell<1) then
- yCell = 1
- end if
- if(yCell>nCellsSide) then
- yCell = nCellsSide
- end if
- if(zCell<1) then
- zCell = 1
- end if
- if(zCell>nCellsSide) then
- zCell = nCellsSide
- end if
-
- part_pos = (xCell-1)*nCellsSide**2 + (yCell-1)*nCellsSide + zCell
-
- call this%storage(part_pos)%addElement( prt, 0.0_dp )
- end do
-
- end subroutine particlesSort
-
-end module SPH_module
diff --git a/sph/louis/src/SPH_simulation.f90 b/sph/louis/src/SPH_simulation.f90
deleted file mode 100644
index f16b20b4eebe6eda7035443cd9e1015e0482bee5..0000000000000000000000000000000000000000
--- a/sph/louis/src/SPH_simulation.f90
+++ /dev/null
@@ -1,31 +0,0 @@
-!> SPH simulation
-!! @n This program is used to solve the Navier-Stokes equations
-!! using the SPH method. A number of files must be given:
-!! paths.txt, *.prm, *.fp and *.mp.
-!! @warning The domain must be cubic!
-!! @brief Main program to launch a SPH simulation
-!! @author Louis Goffin
-!! @date 2013-05-26
-!! @version 1.0.0
-
-program SPH_simulation
-
- use SPH_module
- implicit none
-
- integer :: t1, t2, clock_rate, clock_max
- type(particle_manager) :: manager
-
- print *,'============= SPH_simulation (L.Goffin)'
-
- call system_clock(t1, clock_rate, clock_max)
-
- call manager%initialisation()
- call manager%solver()
-
- call system_clock(t2, clock_rate, clock_max)
-
- print *,'Elapsed real time = ', (t2-t1)/clock_rate
-
-end program SPH_simulation
-
diff --git a/sph/louis/tests/waterdrop.py b/sph/louis/tests/waterdrop.py
deleted file mode 100644
index 19f84b5da57d8505517645c50ac401e04685a6b9..0000000000000000000000000000000000000000
--- a/sph/louis/tests/waterdrop.py
+++ /dev/null
@@ -1,63 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-from sph.helpers import *
-
-if __name__=="__main__":
-
- boxL = 2.
- Lfloor = 0.7
- Lwater = 0.5
- sep = 0.05
-
- kernel = Kernel('cubic', True) # 'cubic', 'quadratic' or 'quintic'
- law = EqState('liquid') # 'gas' or 'liquid'
-
- # parameters
- model = Model()
- model.kernel = kernel
- model.law = law
-
- model.h_0 = 0.06 # initial smoothing length [m]
- model.c_0 = 35.0 # initial speed of sound [m/s]
- model.rho_0 = 1000.0 # initial density [kg/m^3]
- model.dom_dim = boxL # domain size (cube)
- model.alpha = 0.5 # artificial viscosity factor 1
- model.beta = 0.0 # artificial viscosity factor 2
- model.maxTime = 1.0 # simulation time
- model.saveInt = 0.01 # save interval
-
- # mobile particles
- cube = Cube( o=(((boxL-Lwater)/2),((boxL-Lwater)/2), ((boxL)/2)+0.5), L=(Lwater,Lwater,Lwater), rho=model.rho_0, s=sep)
- model.addMobile(cube.generate())
-
- # fixed particles
- plane = Cube( o=(((boxL-Lfloor)/2),((boxL-Lfloor)/2), (boxL/2)), L=(Lfloor,Lfloor,sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- plane = Cube( o=(0,0,0), L=(boxL,boxL,sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
-
- # run SPH model
- print(model)
- model.run()
-
- # convert to VTK
- import sph.gui as gui
- gui.ToParaview(verb=False).convertall()
-
-
\ No newline at end of file
diff --git a/sph/louis/tests/waterdrop2.py b/sph/louis/tests/waterdrop2.py
deleted file mode 100644
index 954fef2261c49f45d9bf8838a009b20326094fd3..0000000000000000000000000000000000000000
--- a/sph/louis/tests/waterdrop2.py
+++ /dev/null
@@ -1,77 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-from sph.helpers import *
-
-if __name__=="__main__":
-
- boxL = 2.
- Lfloor = 0.7
- Lwater = 0.5
- sep = 0.05
-
- kernel = Kernel('cubic', False) # 'cubic', 'quadratic' or 'quintic'
- law = EqState('liquid') # 'gas' or 'liquid'
-
- # parameters
- model = Model()
- model.kernel = kernel
- model.law = law
-
- model.h_0 = 0.06 # initial smoothing length [m]
- model.c_0 = 35.0 # initial speed of sound [m/s]
- model.rho_0 = 1000.0 # initial density [kg/m^3]
- model.dom_dim = boxL # domain size (cube)
- model.alpha = 0.5 # artificial viscosity factor 1
- model.beta = 0.0 # artificial viscosity factor 2
- model.maxTime = 3.0 # simulation time
- model.saveInt = 0.01 # save interval
-
- # mobile particles
- cube = Cube( o=(((boxL-Lwater)/2),((boxL-Lwater)/2), ((boxL)/2)+0.5), L=(Lwater,Lwater,Lwater), rho=model.rho_0, s=sep)
- model.addMobile(cube.generate())
-
- # fixed particles
- #obstacle
- plane = Cube( o=(((boxL-Lfloor)/2),((boxL-Lfloor)/2), (boxL/2)), L=(Lfloor,Lfloor,sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #floor
- plane = Cube( o=(0,0,0), L=(boxL,boxL,sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #x=0
- plane = Cube( o=(0,0,2*sep), L=(sep,boxL,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #y=0
- plane = Cube( o=(2*sep,0,2*sep), L=(boxL-4*sep,sep,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #x=L
- plane = Cube( o=(boxL-sep,0,2*sep), L=(sep,boxL,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #y=L
- plane = Cube( o=(2*sep,boxL-sep,2*sep), L=(boxL-4*sep,sep,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
-
- # run SPH model
- print(model)
- model.run()
-
- # convert to VTK
- import sph.gui as gui
- gui.ToParaview(verb=False).convertall()
-
-
diff --git a/sph/louis/tests/waterdrop3.py b/sph/louis/tests/waterdrop3.py
deleted file mode 100644
index ce6b149727250d51430e68d5e40fede6c29759de..0000000000000000000000000000000000000000
--- a/sph/louis/tests/waterdrop3.py
+++ /dev/null
@@ -1,77 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-from sph.helpers import *
-
-if __name__=="__main__":
-
- boxL = 2.
- Lfloor = 0.7
- Lwater = 0.5
- sep = 0.05/2
-
- kernel = Kernel('cubic', False) # 'cubic', 'quadratic' or 'quintic'
- law = EqState('liquid') # 'gas' or 'liquid'
-
- # parameters
- model = Model()
- model.kernel = kernel
- model.law = law
-
- model.h_0 = 0.06/2 # initial smoothing length [m]
- model.c_0 = 35.0 # initial speed of sound [m/s]
- model.rho_0 = 1000.0 # initial density [kg/m^3]
- model.dom_dim = boxL # domain size (cube)
- model.alpha = 0.5 # artificial viscosity factor 1
- model.beta = 0.0 # artificial viscosity factor 2
- model.maxTime = 3.0 # simulation time
- model.saveInt = 0.01/2 # save interval
-
- # mobile particles
- cube = Cube( o=(((boxL-Lwater)/2),((boxL-Lwater)/2), ((boxL)/2)+0.5), L=(Lwater,Lwater,Lwater), rho=model.rho_0, s=sep)
- model.addMobile(cube.generate())
-
- # fixed particles
- #obstacle
- plane = Cube( o=(((boxL-Lfloor)/2),((boxL-Lfloor)/2), (boxL/2)), L=(Lfloor,Lfloor,sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #floor
- plane = Cube( o=(0,0,0), L=(boxL,boxL,sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #x=0
- plane = Cube( o=(0,0,2*sep), L=(sep,boxL,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #y=0
- plane = Cube( o=(2*sep,0,2*sep), L=(boxL-4*sep,sep,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #x=L
- plane = Cube( o=(boxL-sep,0,2*sep), L=(sep,boxL,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #y=L
- plane = Cube( o=(2*sep,boxL-sep,2*sep), L=(boxL-4*sep,sep,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
-
- # run SPH model
- print(model)
- model.run()
-
- # convert to VTK
- import sph.gui as gui
- gui.ToParaview(verb=False).convertall()
-
-
diff --git a/sph/louis/tests/waterdrop4.py b/sph/louis/tests/waterdrop4.py
deleted file mode 100644
index df08ea33eb252be91830be46a65598a0bc9e4fbc..0000000000000000000000000000000000000000
--- a/sph/louis/tests/waterdrop4.py
+++ /dev/null
@@ -1,77 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-from sph.helpers import *
-
-if __name__=="__main__":
-
- boxL = 2.
- Lfloor = 0.7
- Lwater = 0.5
- sep = 0.05/4
-
- kernel = Kernel('cubic', False) # 'cubic', 'quadratic' or 'quintic'
- law = EqState('liquid') # 'gas' or 'liquid'
-
- # parameters
- model = Model()
- model.kernel = kernel
- model.law = law
-
- model.h_0 = 0.06/4 # initial smoothing length [m]
- model.c_0 = 35.0 # initial speed of sound [m/s]
- model.rho_0 = 1000.0 # initial density [kg/m^3]
- model.dom_dim = boxL # domain size (cube)
- model.alpha = 0.5 # artificial viscosity factor 1
- model.beta = 0.0 # artificial viscosity factor 2
- model.maxTime = 3.0 # simulation time
- model.saveInt = 0.01/2 # save interval
-
- # mobile particles
- cube = Cube( o=(((boxL-Lwater)/2),((boxL-Lwater)/2), ((boxL)/2)+0.5), L=(Lwater,Lwater,Lwater), rho=model.rho_0, s=sep)
- model.addMobile(cube.generate())
-
- # fixed particles
- #obstacle
- plane = Cube( o=(((boxL-Lfloor)/2),((boxL-Lfloor)/2), (boxL/2)), L=(Lfloor,Lfloor,sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #floor
- plane = Cube( o=(0,0,0), L=(boxL,boxL,sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #x=0
- plane = Cube( o=(0,0,2*sep), L=(sep,boxL,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #y=0
- plane = Cube( o=(2*sep,0,2*sep), L=(boxL-4*sep,sep,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #x=L
- plane = Cube( o=(boxL-sep,0,2*sep), L=(sep,boxL,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
- #y=L
- plane = Cube( o=(2*sep,boxL-sep,2*sep), L=(boxL-4*sep,sep,boxL-2*sep), rho=model.rho_0, s=sep)
- model.addFixed(plane.generate())
-
- # run SPH model
- print(model)
- model.run()
-
- # convert to VTK
- import sph.gui as gui
- gui.ToParaview(verb=False).convertall()
-
-
diff --git a/sph/louis/tests/waterdrop_post.py b/sph/louis/tests/waterdrop_post.py
deleted file mode 100644
index 78d57fa3cf6ab30f57ab46e397b8e6014ad54a74..0000000000000000000000000000000000000000
--- a/sph/louis/tests/waterdrop_post.py
+++ /dev/null
@@ -1,188 +0,0 @@
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-#### import the simple module from the paraview
-from paraview.simple import *
-import glob, os.path
-
-#### disable automatic camera reset on 'Show'
-paraview.simple._DisableFirstRenderCameraReset()
-
-print("cwd=", os.getcwd())
-dirname = '/home/boman/dev/Projet_MP/waves/sph/louis/workspace/tests_waterdrop'
-
-# create a new 'XML Structured Grid Reader'
-gridvts = XMLStructuredGridReader(FileName=os.path.join(dirname,'grid.vts'))
-RenameSource('Grid', gridvts)
-
-# create a new 'XML Unstructured Grid Reader'
-inpfiles = sorted(glob.glob(os.path.join(dirname,'resFP_*.vtu')))
-#print inpfiles
-resFP_00000 = XMLUnstructuredGridReader(FileName=inpfiles)
-resFP_00000.PointArrayStatus = ['max(mu_ab)', 'Pressure', 'Speed of sound', 'Mass density', 'Velocity', 'Smoothing length', 'Mass', 'Nb of neighbours']
-RenameSource('Fixed particles', resFP_00000)
-
-# get animation scene
-animationScene1 = GetAnimationScene()
-
-# update animation scene based on data timesteps
-animationScene1.UpdateAnimationUsingDataTimeSteps()
-
-# create a new 'XML Unstructured Grid Reader'
-inpfiles = sorted(glob.glob(os.path.join(dirname,'resMP_*.vtu')))
-resMP_00000 = XMLUnstructuredGridReader(FileName=inpfiles)
-resMP_00000.PointArrayStatus = ['max(mu_ab)', 'Pressure', 'Speed of sound', 'Mass density', 'Velocity', 'Smoothing length', 'Mass', 'Nb of neighbours']
-RenameSource('Mobile particles', resMP_00000)
-
-# get active view
-renderView1 = GetActiveViewOrCreate('RenderView')
-# uncomment following to set a specific view size
-# renderView1.ViewSize = [997, 849]
-
-# show data in view
-gridvtsDisplay = Show(gridvts, renderView1)
-# trace defaults for the display properties.
-gridvtsDisplay.Representation = 'Outline'
-gridvtsDisplay.ColorArrayName = ['POINTS', '']
-gridvtsDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
-gridvtsDisplay.SelectOrientationVectors = 'None'
-gridvtsDisplay.ScaleFactor = 0.2
-gridvtsDisplay.SelectScaleArray = 'None'
-gridvtsDisplay.GlyphType = 'Arrow'
-gridvtsDisplay.ScalarOpacityUnitDistance = 0.21650635094610968
-
-# reset view to fit data
-renderView1.ResetCamera()
-
-# show data in view
-resFP_00000Display = Show(resFP_00000, renderView1)
-# trace defaults for the display properties.
-resFP_00000Display.ColorArrayName = [None, '']
-resFP_00000Display.OSPRayScaleArray = 'Mass'
-resFP_00000Display.OSPRayScaleFunction = 'PiecewiseFunction'
-resFP_00000Display.SelectOrientationVectors = 'Mass'
-resFP_00000Display.ScaleFactor = 0.19512200355529785
-resFP_00000Display.SelectScaleArray = 'Mass'
-resFP_00000Display.GlyphType = 'Arrow'
-resFP_00000Display.ScalarOpacityUnitDistance = 0.18843877254666028
-resFP_00000Display.GaussianRadius = 0.09756100177764893
-resFP_00000Display.SetScaleArray = ['POINTS', 'Mass']
-resFP_00000Display.ScaleTransferFunction = 'PiecewiseFunction'
-resFP_00000Display.OpacityArray = ['POINTS', 'Mass']
-resFP_00000Display.OpacityTransferFunction = 'PiecewiseFunction'
-
-# show data in view
-resMP_00000Display = Show(resMP_00000, renderView1)
-# trace defaults for the display properties.
-resMP_00000Display.ColorArrayName = [None, '']
-resMP_00000Display.OSPRayScaleArray = 'Mass'
-resMP_00000Display.OSPRayScaleFunction = 'PiecewiseFunction'
-resMP_00000Display.SelectOrientationVectors = 'Mass'
-resMP_00000Display.ScaleFactor = 0.04545450210571289
-resMP_00000Display.SelectScaleArray = 'Mass'
-resMP_00000Display.GlyphType = 'Arrow'
-resMP_00000Display.ScalarOpacityUnitDistance = 0.07157227916349206
-resMP_00000Display.GaussianRadius = 0.022727251052856445
-resMP_00000Display.SetScaleArray = ['POINTS', 'Mass']
-resMP_00000Display.ScaleTransferFunction = 'PiecewiseFunction'
-resMP_00000Display.OpacityArray = ['POINTS', 'Mass']
-resMP_00000Display.OpacityTransferFunction = 'PiecewiseFunction'
-
-# set active source
-SetActiveSource(gridvts)
-
-# change representation type
-gridvtsDisplay.SetRepresentationType('Wireframe')
-
-
-# Properties modified on gridvtsDisplay
-gridvtsDisplay.Opacity = 0.2
-
-# set active source
-SetActiveSource(resFP_00000)
-
-# change representation type
-resFP_00000Display.SetRepresentationType('Points')
-
-# Properties modified on resFP_00000Display
-resFP_00000Display.PointSize = 5.0
-
-# set active source
-SetActiveSource(resMP_00000)
-
-# create a new 'Glyph'
-glyph1 = Glyph(Input=resMP_00000, GlyphType='Sphere')
-glyph1.Scalars = ['POINTS', 'Smoothing length']
-glyph1.Vectors = ['POINTS', 'None']
-glyph1.ScaleMode = 'scalar'
-glyph1.ScaleFactor = 1.0
-glyph1.GlyphTransform = 'Transform2'
-glyph1.GlyphMode = 'All Points'
-
-
-
-# get color transfer function/color map for 'Smoothinglength'
-smoothinglengthLUT = GetColorTransferFunction('Smoothinglength')
-
-# show data in view
-glyph1Display = Show(glyph1, renderView1)
-# trace defaults for the display properties.
-glyph1Display.ColorArrayName = ['POINTS', 'Smoothing length']
-glyph1Display.LookupTable = smoothinglengthLUT
-glyph1Display.OSPRayScaleArray = 'Smoothing length'
-glyph1Display.OSPRayScaleFunction = 'PiecewiseFunction'
-glyph1Display.SelectOrientationVectors = 'Mass'
-glyph1Display.ScaleFactor = 0.051454496383666996
-glyph1Display.SelectScaleArray = 'Smoothing length'
-glyph1Display.GlyphType = 'Arrow'
-glyph1Display.GaussianRadius = 0.025727248191833498
-glyph1Display.SetScaleArray = ['POINTS', 'Smoothing length']
-glyph1Display.ScaleTransferFunction = 'PiecewiseFunction'
-glyph1Display.OpacityArray = ['POINTS', 'Smoothing length']
-glyph1Display.OpacityTransferFunction = 'PiecewiseFunction'
-
-# show color bar/color legend
-glyph1Display.SetScalarBarVisibility(renderView1, True)
-
-# get opacity transfer function/opacity map for 'Smoothinglength'
-smoothinglengthPWF = GetOpacityTransferFunction('Smoothinglength')
-
-# hide data in view
-Hide(resMP_00000, renderView1)
-
-# set scalar coloring
-ColorBy(glyph1Display, ('POINTS', 'Velocity'))
-
-# Hide the scalar bar for this color map if no visible data is colored by it.
-HideScalarBarIfNotNeeded(smoothinglengthLUT, renderView1)
-
-# rescale color and/or opacity maps used to include current data range
-glyph1Display.RescaleTransferFunctionToDataRange(True, False)
-
-# show color bar/color legend
-glyph1Display.SetScalarBarVisibility(renderView1, True)
-
-
-#### saving camera placements for all active views
-
-# current camera placement for renderView1
-renderView1.CameraPosition = [2.242394659003803, -5.320849062052243, 2.8132656553463735]
-renderView1.CameraFocalPoint = [1.0000000000000002, 0.9999999999999994, 0.9999999999999989]
-renderView1.CameraViewUp = [-0.2489518711661434, 0.22153968976360192, 0.9428378077391271]
-renderView1.CameraParallelScale = 1.7320508075688772
-
-#### uncomment the following to render all views
-RenderAllViews()
-# alternatively, if you want to write images, you can use SaveScreenshot(...).
\ No newline at end of file
diff --git a/sph/matlab/interp.m b/sph/matlab/interp.m
deleted file mode 100644
index 4e1b3f06517d75498caf9127359ef0424029fcb6..0000000000000000000000000000000000000000
--- a/sph/matlab/interp.m
+++ /dev/null
@@ -1,72 +0,0 @@
-
-function interp()
- close all; clear all;
-
- % plot M4 spline kernel
- x=linspace(-3,3,100);
- figure;
- h=0.5;
- plot(x,m4spline(x,h));
- hold on
- plot(x,gausskernel(x,h));
- grid;
- xlabel('x')
- ylabel('kernel')
- legend('M4 spline', 'gaussian')
-
- % verify that int W = 1
- %integral(@(x)m4spline(x,h),-3*h,3*h)
- integral(@(x)gausskernel(x,h),-Inf,Inf)
-
- % plot a sin and the approximate fct
- x=linspace(1,10,200);
- y=myfun(x);
-
- figure;
- plot(x,y);
- xlabel('x')
- grid;
- hold on;
- h=2;
- plot(x,approx(x,h))
- h=0.5;
- plot(x,approx(x,h))
- legend('exact','h=2','h=0.5');
-
-end
-
-function y=myfun(x)
- %y=(x-5).^2/10;
- y=sin(x);
-end
-
-
-function v=m4spline(x,h)
- r=abs(x)/h;
- i1=find(r<1);
- i2=find(r>=1 & r<2);
- v=zeros(size(r));
- v(i1) = (((2-r(i1)).^3)-4*((1-r(i1)).^3))/6;
- v(i2) = ((2-r(i2)).^3)/6;
- v=v/h;
-end
-
-function v=gausskernel(x,h)
- v = exp(-(x.*x)/(h^2))/(pi^(1/2)*h); % 1D
-end
-
-
-function v=approx(x,h)
- f = @(xp,c) m4spline(c-xp, h).*myfun(xp);
- %f = @(xp,c) gausskernel(c-xp, h).*myfun(xp);
-
-
- for i=1:length(x)
- v(i) = integral(@(xp)f(xp,x(i)),x(i)-3*h,x(i)+3*h);
- end
-end
-
-
-
-
-
\ No newline at end of file
diff --git a/sph/src/CMakeLists.txt b/sph/src/CMakeLists.txt
deleted file mode 100644
index becbe3e1677242c3b05c79497b5fe872a8cc236f..0000000000000000000000000000000000000000
--- a/sph/src/CMakeLists.txt
+++ /dev/null
@@ -1,26 +0,0 @@
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-# CMake input file of sph.so
-
-FILE(GLOB SRCS *.h *.cpp *.inl *.hpp)
-
-ADD_LIBRARY(sph SHARED ${SRCS})
-MACRO_DebugPostfix(sph)
-
-TARGET_INCLUDE_DIRECTORIES(sph PUBLIC ${PROJECT_SOURCE_DIR}/sph/src)
-
-TARGET_LINK_LIBRARIES(sph tbox fwk)
-
-SOURCE_GROUP(base REGULAR_EXPRESSION ".*\\.(cpp|inl|hpp|h)")
diff --git a/sph/src/sph.h b/sph/src/sph.h
deleted file mode 100644
index 567ecf581a6f4f4271a8cb35fa4b4292314ec988..0000000000000000000000000000000000000000
--- a/sph/src/sph.h
+++ /dev/null
@@ -1,52 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-// global header of the sph module
-
-#ifndef SPH_H
-#define SPH_H
-
-#if defined(WIN32)
-#ifdef sph_EXPORTS
-#define SPH_API __declspec(dllexport)
-#else
-#define SPH_API __declspec(dllimport)
-#endif
-#else
-#define SPH_API
-#endif
-
-#include "tbox.h"
-#include <math.h>
-
-namespace sph
-{
-class Dofs;
-class Problem;
-class Particle;
-class FixedParticle;
-class MobileParticle;
-class EqState;
-class IdealGas;
-class QIncFluid;
-class Kernel;
-class CubicSplineKernel;
-class QuadraticKernel;
-class QuinticSplineKernel;
-class DisplayHook;
-} // namespace sph
-
-#endif //SPH_H
diff --git a/sph/src/wDisplayHook.cpp b/sph/src/wDisplayHook.cpp
deleted file mode 100644
index c27f6b9be0814210fd1d33d532202de2e1453b0c..0000000000000000000000000000000000000000
--- a/sph/src/wDisplayHook.cpp
+++ /dev/null
@@ -1,32 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wDisplayHook.h"
-using namespace sph;
-
-DisplayHook::DisplayHook()
-{
-}
-
-void DisplayHook::display(int nt, double t, std::vector<double> &u)
-{
- std::cout << "DisplayHook::display()\n";
-}
-
-void DisplayHook::refresh()
-{
- std::cout << "DisplayHook::refresh()\n";
-}
diff --git a/sph/src/wDisplayHook.h b/sph/src/wDisplayHook.h
deleted file mode 100644
index 7489d09e18b15e0bb8a4349c375bdcbfd89994d4..0000000000000000000000000000000000000000
--- a/sph/src/wDisplayHook.h
+++ /dev/null
@@ -1,41 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WDISPLAYHOOK_H
-#define WDISPLAYHOOK_H
-
-#include "sph.h"
-#include "wObject.h"
-#include <vector>
-
-namespace sph
-{
-
-/**
- * @brief Dispay Hook - links with the python GUI
- */
-
-class SPH_API DisplayHook : public fwk::wObject
-{
-public:
- DisplayHook();
- virtual void display(int nt, double t, std::vector<double> &u);
- virtual void refresh();
-};
-
-} // namespace sph
-
-#endif //WDISPLAYHOOK_H
diff --git a/sph/src/wDofs.cpp b/sph/src/wDofs.cpp
deleted file mode 100644
index 65e2fd0e8a1ed39235b45eb3dd4abec05119dd64..0000000000000000000000000000000000000000
--- a/sph/src/wDofs.cpp
+++ /dev/null
@@ -1,28 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wDofs.h"
-using namespace sph;
-
-namespace sph
-{
-SPH_API std::ostream &
-operator<<(std::ostream &out, Dofs const &obj)
-{
- out << "x=" << obj.x.transpose() << " v=" << obj.u.transpose() << " rho=" << obj.rho;
- return out;
-}
-} // namespace sph
diff --git a/sph/src/wDofs.h b/sph/src/wDofs.h
deleted file mode 100644
index d86af3e42915bee46449e540061216cbfd5dd8fc..0000000000000000000000000000000000000000
--- a/sph/src/wDofs.h
+++ /dev/null
@@ -1,47 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WDOFS_H
-#define WDOFS_H
-
-#include "sph.h"
-#include <iostream>
-#include <vector>
-#include <Eigen/Core>
-
-namespace sph
-{
-
-/**
- * @brief independant variables of the problem
- */
-
-class SPH_API Dofs
-{
-public:
- Eigen::Vector3d x; //< position
- Eigen::Vector3d u; //< velocity
- double rho; //< density
- Dofs(Eigen::Vector3d const &_x, Eigen::Vector3d const &_u, double _rho) : x(_x), u(_u), rho(_rho) {}
-
-#ifndef SWIG
- friend SPH_API std::ostream &operator<<(std::ostream &out, Dofs const &obj);
-#endif
-};
-
-} // namespace sph
-
-#endif //WDOFS_H
diff --git a/sph/src/wEigenTest.cpp b/sph/src/wEigenTest.cpp
deleted file mode 100644
index 0ec8bcece263f9c8e30b536915aad39a455c731b..0000000000000000000000000000000000000000
--- a/sph/src/wEigenTest.cpp
+++ /dev/null
@@ -1,27 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wEigenTest.h"
-#include <Eigen/Eigen>
-#include <iostream>
-
-using namespace sph;
-using namespace Eigen;
-void EigenTest::test1()
-{
- Matrix4f A = Matrix4f::Random();
- std::cout << A << std::endl;
-}
diff --git a/sph/src/wEigenTest.h b/sph/src/wEigenTest.h
deleted file mode 100644
index c8bf5d3f5e034f9db17413a35255898b23e48331..0000000000000000000000000000000000000000
--- a/sph/src/wEigenTest.h
+++ /dev/null
@@ -1,39 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WEIGENTEST_H
-#define WEIGENTEST_H
-
-#include "sph.h"
-
-namespace sph
-{
-
-/**
- * @brief test of the Eigen library
- */
-
-class SPH_API EigenTest
-{
-public:
- EigenTest() {}
-
- void test1();
-};
-
-} // namespace sph
-
-#endif //WEIGENTEST_H
diff --git a/sph/src/wEqState.cpp b/sph/src/wEqState.cpp
deleted file mode 100644
index 0fb78e336ce9126b21319baa332f979cc29fdb57..0000000000000000000000000000000000000000
--- a/sph/src/wEqState.cpp
+++ /dev/null
@@ -1,63 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wEqState.h"
-using namespace sph;
-
-EqState::EqState()
-{
-}
-
-EqState::~EqState()
-{
-}
-
-void EqState::write(std::ostream &out) const
-{
- out << "sph::EqState";
-}
-
-// --------------------------------------------------------
-
-IdealGas::IdealGas()
-{
-}
-
-IdealGas::~IdealGas()
-{
- std::cout << "~IdealGas()\n";
-}
-
-void IdealGas::write(std::ostream &out) const
-{
- out << "sph::IdealGas";
-}
-
-// --------------------------------------------------------
-
-QIncFluid::QIncFluid()
-{
-}
-
-QIncFluid::~QIncFluid()
-{
- std::cout << "~QIncFluid()\n";
-}
-
-void QIncFluid::write(std::ostream &out) const
-{
- out << "sph::QIncFluid";
-}
diff --git a/sph/src/wEqState.h b/sph/src/wEqState.h
deleted file mode 100644
index a4a19b24946ce28f81d23615086452fb5216e81d..0000000000000000000000000000000000000000
--- a/sph/src/wEqState.h
+++ /dev/null
@@ -1,76 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WEQSTATE_H
-#define WEQSTATE_H
-
-#include "sph.h"
-#include "wObject.h"
-#include <iostream>
-#include <vector>
-#include <memory>
-
-namespace sph
-{
-
-/**
- * @brief Virtual class for Equations of state
- */
-
-class SPH_API EqState : public fwk::wSharedObject
-{
-public:
- EqState();
- virtual ~EqState();
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-/**
- * @brief Ideal gas
- */
-
-class SPH_API IdealGas : public EqState
-{
-public:
- IdealGas();
- virtual ~IdealGas();
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-/**
- * @brief Quasi incompressible fluid
- */
-
-class SPH_API QIncFluid : public EqState
-{
-public:
- QIncFluid();
- virtual ~QIncFluid();
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-} // namespace sph
-
-#endif //WEQSTATE_H
diff --git a/sph/src/wKernel.cpp b/sph/src/wKernel.cpp
deleted file mode 100644
index fc8129c247f8fc5348fc22a431c54a917534aa04..0000000000000000000000000000000000000000
--- a/sph/src/wKernel.cpp
+++ /dev/null
@@ -1,175 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wKernel.h"
-using namespace sph;
-
-Kernel::Kernel()
-{
-}
-
-Kernel::~Kernel()
-{
-}
-
-void Kernel::write(std::ostream &out) const
-{
- out << "sph::Kernel";
-}
-
-// -----------------------------------------------------------
-
-CubicSplineKernel::CubicSplineKernel() : Kernel()
-{
-}
-
-CubicSplineKernel::~CubicSplineKernel()
-{
- std::cout << "~CubicSplineKernel()\n";
-}
-
-void CubicSplineKernel::write(std::ostream &out) const
-{
- out << "sph::CubicSplineKernel";
-}
-
-double
-CubicSplineKernel::value(double r, double h) const
-{
- double h2 = h * h;
- double h3 = h2 * h;
- double alphad = 3.0 / (2.0 * M_PI * h3);
- double rh = r / h;
-
- if (rh < 1.0)
- return alphad * (3.0 / 2.0 - rh * rh + 0.5 * rh * rh * rh);
- else if (rh < 2.0)
- return alphad * (1.0 / 6.0 * pow(1.0 - rh, 3));
- else
- return 0.0;
-}
-
-double
-CubicSplineKernel::grad(double r, double h) const
-{
- double h2 = h * h;
- double h3 = h2 * h;
- double alphad = 3.0 / (2.0 * M_PI * h3);
- double rh = r / h;
-
- if (rh < 1.0)
- return alphad / h * (3.0 / 2.0 * rh * rh - 2.0 * rh);
- else if (rh < 2.0)
- return alphad / h * (-0.5 * (2.0 - rh) * (2.0 - rh));
- else
- return 0.0;
-}
-
-// -----------------------------------------------------------
-
-QuadraticKernel::QuadraticKernel() : Kernel()
-{
-}
-
-QuadraticKernel::~QuadraticKernel()
-{
- std::cout << "~QuadraticKernel()\n";
-}
-
-void QuadraticKernel::write(std::ostream &out) const
-{
- out << "sph::QuadraticKernel";
-}
-
-double
-QuadraticKernel::value(double r, double h) const
-{
- double h2 = h * h;
- double h3 = h2 * h;
- double alphad = 5.0 / (4.0 * M_PI * h3);
- double rh = r / h;
-
- if (rh < 2.0)
- return alphad * (3.0 / 16.0 * rh * rh - 3.0 / 4.0 * rh + 3.0 / 4.0);
- else
- return 0.0;
-}
-
-double
-QuadraticKernel::grad(double r, double h) const
-{
- double h2 = h * h;
- double h3 = h2 * h;
- double alphad = 5.0 / (4.0 * M_PI * h3);
- double rh = r / h;
-
- if (rh < 2.0)
- return alphad / h * (3.0 / 8.0 * rh - 3.0 / 4.0);
- else
- return 0.0;
-}
-
-// -----------------------------------------------------------
-
-QuinticSplineKernel::QuinticSplineKernel() : Kernel()
-{
-}
-
-QuinticSplineKernel::~QuinticSplineKernel()
-{
- std::cout << "~QuinticSplineKernel()\n";
-}
-
-void QuinticSplineKernel::write(std::ostream &out) const
-{
- out << "sph::QuinticSplineKernel";
-}
-
-double
-QuinticSplineKernel::value(double r, double h) const
-{
- double h2 = h * h;
- double h3 = h2 * h;
- double alphad = 3.0 / (359.0 * M_PI * h3);
- double rh = r / h;
-
- if (rh < 1.0)
- return alphad * (pow(3.0 - rh, 5) - 6.0 * pow(2.0 - rh, 5) + 15.0 * pow(1.0 - rh, 5));
- else if (rh < 2.0)
- return alphad * (pow(3.0 - rh, 5) - 6.0 * pow(2.0 - rh, 5));
- else if (rh < 3.0)
- return alphad * (pow(3.0 - rh, 5));
- else
- return 0.0;
-}
-
-double
-QuinticSplineKernel::grad(double r, double h) const
-{
- double h2 = h * h;
- double h3 = h2 * h;
- double alphad = 3.0 / (359.0 * M_PI * h3);
- double rh = r / h;
-
- if (rh < 1.0)
- return alphad / h * (-5.0 * pow(3.0 - rh, 4) + 30.0 * pow(2.0 - rh, 4) - 75.0 * pow(1.0 - rh, 4));
- else if (rh < 2.0)
- return alphad / h * (-5.0 * pow(3.0 - rh, 4) + 30.0 * pow(2.0 - rh, 4));
- else if (rh < 3.0)
- return alphad / h * (-5.0 * pow(3.0 - rh, 4));
- else
- return 0.0;
-}
diff --git a/sph/src/wKernel.h b/sph/src/wKernel.h
deleted file mode 100644
index aed4e0b9661a3ecba047798fcddb958e6d38f1fd..0000000000000000000000000000000000000000
--- a/sph/src/wKernel.h
+++ /dev/null
@@ -1,104 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WKERNEL_H
-#define WKERNEL_H
-
-#include "sph.h"
-#include "wObject.h"
-#include <iostream>
-#include <vector>
-#include <memory>
-
-namespace sph
-{
-
-/**
- * @brief virtual Kernel class
- */
-
-class SPH_API Kernel : public fwk::wSharedObject
-{
-public:
- Kernel();
- virtual ~Kernel();
- virtual double kappa() const = 0;
-
- virtual double value(double r, double h) const = 0;
- virtual double grad(double r, double h) const = 0;
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-/**
- * @brief Cubic spline kernel
- */
-
-class SPH_API CubicSplineKernel : public Kernel
-{
-public:
- CubicSplineKernel();
- virtual ~CubicSplineKernel();
- virtual double kappa() const override { return 2.0; }
- virtual double value(double r, double h) const override;
- virtual double grad(double r, double h) const override;
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-/**
- * @brief Quadratic kernel
- */
-
-class SPH_API QuadraticKernel : public Kernel
-{
-public:
- QuadraticKernel();
- virtual ~QuadraticKernel();
- virtual double kappa() const override { return 2.0; }
- virtual double value(double r, double h) const override;
- virtual double grad(double r, double h) const override;
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-/**
- * @brief Quintic spline kernel
- */
-
-class SPH_API QuinticSplineKernel : public Kernel
-{
-public:
- QuinticSplineKernel();
- virtual ~QuinticSplineKernel();
- virtual double kappa() const override { return 3.0; }
- virtual double value(double r, double h) const override;
- virtual double grad(double r, double h) const override;
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-} // namespace sph
-
-#endif //WKERNEL_H
diff --git a/sph/src/wParticle.cpp b/sph/src/wParticle.cpp
deleted file mode 100644
index 20841401ba33f9fdfb6b489519f6e92699e2c5a6..0000000000000000000000000000000000000000
--- a/sph/src/wParticle.cpp
+++ /dev/null
@@ -1,66 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wParticle.h"
-#include "wDofs.h"
-using namespace sph;
-
-Particle::Particle(int _no, Dofs const &_dofs) : no(_no)
-{
- dofs.push_back(new Dofs(_dofs));
-}
-
-Particle::~Particle()
-{
- std::cout << "~Particle()\n";
- for (auto d : dofs)
- delete d;
-}
-
-void Particle::write(std::ostream &out) const
-{
- out << "p#" << no;
-}
-
-// ------------------------------------------------------------------------
-
-FixedParticle::FixedParticle(int _no, Dofs const &_dofs) : Particle(_no, _dofs)
-{
-}
-
-FixedParticle::~FixedParticle()
-{
-}
-
-void FixedParticle::write(std::ostream &out) const
-{
- out << "pF#" << no;
-}
-
-// ------------------------------------------------------------------------
-
-MobileParticle::MobileParticle(int _no, Dofs const &_dofs) : Particle(_no, _dofs)
-{
-}
-
-MobileParticle::~MobileParticle()
-{
-}
-
-void MobileParticle::write(std::ostream &out) const
-{
- out << "pM#" << no;
-}
diff --git a/sph/src/wParticle.h b/sph/src/wParticle.h
deleted file mode 100644
index f0c0ba113034e07a20a6be8be959799e4a139123..0000000000000000000000000000000000000000
--- a/sph/src/wParticle.h
+++ /dev/null
@@ -1,85 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WPARTICLE_H
-#define WPARTICLE_H
-
-#include "sph.h"
-#include "wObject.h"
-#include <iostream>
-#include <vector>
-#include <memory>
-
-using namespace tbox;
-
-namespace sph
-{
-
-/**
- * @brief a particle
- */
-
-class SPH_API Particle : public fwk::wObject
-{
-public:
- int no; //< label (debug/visu)
- std::vector<Dofs *> dofs; //< variables that are time integrated
- double p; //< pressure
- double c; //< speed of sound
- std::vector<Particle *> neighs; //< list of neighbours
-
-public:
- Particle(int _no, Dofs const &_dofs);
- virtual ~Particle();
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-/**
- * @brief Fixed particle
- */
-
-class SPH_API FixedParticle : public Particle
-{
-public:
- FixedParticle(int _no, Dofs const &_dofs);
- virtual ~FixedParticle();
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-/**
- * @brief Mobile particle
- */
-
-class SPH_API MobileParticle : public Particle
-{
-public:
- MobileParticle(int _no, Dofs const &_dofs);
- virtual ~MobileParticle();
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-} // namespace sph
-
-#endif //WPARTICLE_H
diff --git a/sph/src/wProblem.cpp b/sph/src/wProblem.cpp
deleted file mode 100644
index 96f54c162bb85eb945fd5a1d5254c3a38a85829b..0000000000000000000000000000000000000000
--- a/sph/src/wProblem.cpp
+++ /dev/null
@@ -1,44 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wProblem.h"
-#include "wParticle.h"
-using namespace sph;
-
-Problem::Problem(double _h) : h(_h)
-{
-}
-
-Problem::~Problem()
-{
- std::cout << "~Problem()\n";
-}
-
-void Problem::write(std::ostream &out) const
-{
- out << "sph::Problem:\n";
- out << "\t" << prts.size() << " particles";
-}
-
-void Problem::addFixedP(int no, Dofs const &dof)
-{
- prts.push_back(new FixedParticle(no, dof));
-}
-
-void Problem::addMobileP(int no, Dofs const &dof)
-{
- prts.push_back(new MobileParticle(no, dof));
-}
diff --git a/sph/src/wProblem.h b/sph/src/wProblem.h
deleted file mode 100644
index d57a982358f2767c5bafe2f2db6b802c781dbdfb..0000000000000000000000000000000000000000
--- a/sph/src/wProblem.h
+++ /dev/null
@@ -1,53 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WPROBLEM_H
-#define WPROBLEM_H
-
-#include "sph.h"
-#include "wObject.h"
-#include <iostream>
-#include <vector>
-#include <memory>
-
-namespace sph
-{
-
-/**
- * @brief main sph Problem object
- */
-
-class SPH_API Problem : public fwk::wSharedObject
-{
-public:
- std::vector<Particle *> prts; //< set containing all the particles
- double h; //< smoothing length (same for all particles) but may vary with time
-
- Problem(double _h);
- virtual ~Problem();
-
- // -- add particles to the problem
- void addFixedP(int no, Dofs const &dof);
- void addMobileP(int no, Dofs const &dof);
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-} // namespace sph
-
-#endif //WPROBLEM_H
diff --git a/sph/src/wSorter.cpp b/sph/src/wSorter.cpp
deleted file mode 100644
index 6de12df947115dd233d4b06e437e4eadb5483779..0000000000000000000000000000000000000000
--- a/sph/src/wSorter.cpp
+++ /dev/null
@@ -1,205 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wSorter.h"
-#include "wKernel.h"
-#include "wParticle.h"
-#include "wDofs.h"
-#include <algorithm>
-#include <sstream>
-#include <exception>
-#include <tbb/parallel_for.h>
-using namespace sph;
-
-Sorter::Sorter() : nthreads(1)
-{
-}
-
-Sorter::~Sorter()
-{
-}
-
-void Sorter::write(std::ostream &out) const
-{
- out << "sph::Sorter";
-}
-
-// ------------------------------------
-
-LouisSorter::LouisSorter(double _edgeL, double _hmax, Kernel const &kernel) : edgeL(_edgeL), hmax(_hmax)
-{
- kappa = kernel.kappa();
- ncells = std::max(int(edgeL / (hmax * kappa)), 1);
-}
-
-LouisSorter::~LouisSorter()
-{
- std::cout << "~LouisSorter()\n";
-}
-
-void LouisSorter::write(std::ostream &out) const
-{
- out << "sph::LouisSorter:\n";
- out << "\tedgeL = " << edgeL << '\n';
- out << "\thmax = " << hmax << '\n';
- out << "\tncells = " << ncells << 'x' << ncells << 'x' << ncells;
-}
-
-void LouisSorter::execute(std::vector<Particle *> &prts, int step)
-{
- // alloc if required
- if (sgrid.empty())
- sgrid.resize(ncells * ncells * ncells);
- else
- for (auto &l : sgrid)
- l.clear();
-
- // fills the grid
- for (auto p : prts)
- {
- Eigen::Vector3d &x = p->dofs[step]->x;
- int ii[3];
- for (int i = 0; i < 3; ++i)
- {
- int id = int((x(i) / edgeL) * ncells);
- if (id < 0)
- id = 0;
- if (id >= ncells)
- id = ncells - 1;
- ii[i] = id;
- }
- sgrid[idx(ii[0], ii[1], ii[2])].push_back(p);
- }
-
- // fills neighbours
- double Rmax = (hmax * kappa);
- double R2max = Rmax * Rmax;
-
- for (auto p : prts)
- {
- p->neighs.clear();
- Eigen::Vector3d &x = p->dofs[step]->x;
- int binf[3], bsup[3];
- for (int i = 0; i < 3; ++i)
- {
- int id = int((x(i) / edgeL) * ncells);
- if (id < 0)
- id = 0;
- if (id >= ncells)
- id = ncells - 1;
- binf[i] = decidx(id);
- bsup[i] = incidx(id);
- }
- for (int i = binf[0]; i <= bsup[0]; ++i)
- for (int j = binf[1]; j <= bsup[1]; ++j)
- for (int k = binf[2]; k <= bsup[2]; ++k)
- for (auto pn : sgrid[idx(i, j, k)])
- if (pn != p)
- {
- Eigen::Vector3d &xn = pn->dofs[step]->x;
- double R2 = (xn - x).squaredNorm();
- if (R2 <= R2max)
- p->neighs.push_back(pn);
- }
- }
-}
-
-void LouisSorter::chkidx(int i, int j, int k) const
-{
- if (i < 0 || i >= ncells || j < 0 || j >= ncells || k < 0 || k >= ncells)
- {
- std::stringstream str;
- str << "LouisSorter:: bad index (" << i << "," << j << "," << k << ") [ncells=" << ncells << "]!";
- throw std::runtime_error(str.str());
- }
-}
-
-// ------------------------------------
-
-BruteForceSorter::BruteForceSorter(double _hmax, Kernel const &kernel) : hmax(_hmax)
-{
- kappa = kernel.kappa();
- method = 1;
-}
-
-BruteForceSorter::~BruteForceSorter()
-{
- std::cout << "~BruteForceSorter()\n";
-}
-
-void BruteForceSorter::write(std::ostream &out) const
-{
- out << "sph::BruteForceSorter:\n";
- out << "\thmax = " << hmax << '\n';
-}
-
-void BruteForceSorter::execute(std::vector<Particle *> &prts, int step)
-{
- for (auto p : prts)
- p->neighs.clear();
- //p->neighs.reserve(100); // change rien
-
- double Rmax = (hmax * kappa);
- double R2max = Rmax * Rmax;
-
- if (method == 1)
- {
- std::cout << "using serial code\n";
- for (auto it = prts.begin(); (it + 1) != prts.end(); ++it)
- {
- Particle *p = *it;
- Eigen::Vector3d &x = p->dofs[step]->x;
- for (auto it2 = it + 1; it2 != prts.end(); ++it2)
- {
- Particle *pn = *it2;
- Eigen::Vector3d &xn = pn->dofs[step]->x;
- double R2 = (xn - x).squaredNorm();
- if (R2 <= R2max)
- {
- p->neighs.push_back(pn);
- pn->neighs.push_back(p);
- }
- }
- }
- }
- else //if(method==2)
- {
- std::cout << "using parallel code with " << nthreads << " threads\n";
- // (NAIVE) PARALLEL CODE
-
- int grainsize = 200;
- tbb::parallel_for(
- tbb::blocked_range<size_t>(0, prts.size(), grainsize),
- [=](tbb::blocked_range<size_t> &r) {
- //std::cout << r.size() << '\n';
- for (size_t i = r.begin(); i < r.end(); ++i)
- {
- Particle *p = prts[i];
- Eigen::Vector3d &x = p->dofs[step]->x;
- //Eigen::Vector3d x = p->dofs[step]->x;
- for (auto it2 = prts.begin(); it2 != prts.end(); ++it2)
- {
- Particle *pn = *it2;
- Eigen::Vector3d &xn = pn->dofs[step]->x;
- double R2 = (xn - x).squaredNorm();
- if (R2 <= R2max)
- p->neighs.push_back(pn);
- }
- }
- },
- tbb::simple_partitioner());
- }
-}
diff --git a/sph/src/wSorter.h b/sph/src/wSorter.h
deleted file mode 100644
index c249c6f922519f2d782bc63330bbd7f7803d5886..0000000000000000000000000000000000000000
--- a/sph/src/wSorter.h
+++ /dev/null
@@ -1,123 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WSORTER_H
-#define WSORTER_H
-
-#include "sph.h"
-#include "wObject.h"
-#include <iostream>
-#include <vector>
-#include <memory>
-
-namespace sph
-{
-
-/**
- * @brief Nearest neighbour search
- */
-
-class SPH_API Sorter : public fwk::wSharedObject
-{
-public:
- int nthreads; //< nb of threads
-
- Sorter();
- virtual ~Sorter();
-
- virtual void execute(std::vector<Particle *> &parts, int step = 0) = 0;
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-/**
- * @brief Louis'method: particles are sorted in a cubic structured grid
- * between (x,y,z)=(0,0,0) and (x,y,z)=(edgeL,edgeL,edgeL)
- */
-
-class SPH_API LouisSorter : public Sorter
-{
- std::vector<std::vector<Particle *>> sgrid; //< structured grid
- double kappa;
-
-public:
- double edgeL; //< length of the cube edge
- double hmax; //< expected max smothing length used to calculate the number of elements in the sgrid
- int ncells; //< number of cells along 1 edge
-
- LouisSorter(double _edgeL, double _hmax, Kernel const &kernel);
- virtual ~LouisSorter();
- virtual void execute(std::vector<Particle *> &prts, int step = 0) override;
-
- // debug (only)
- size_t nbprts(int i, int j, int k) const
- {
- chkidx(i, j, k);
- return sgrid[idx(i, j, k)].size();
- }
- std::vector<Particle *> &getParticles(int i, int j, int k)
- {
- chkidx(i, j, k);
- return sgrid[idx(i, j, k)];
- }
- int size() const { return ncells; }
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-private:
- int idx(int i, int j, int k) const { return ncells * (ncells * i + j) + k; }
- void chkidx(int i, int j, int k) const;
- int incidx(int id) const
- {
- ++id;
- if (id >= ncells)
- id = ncells - 1;
- return id;
- }
- int decidx(int id) const
- {
- --id;
- if (id < 0)
- id = 0;
- return id;
- }
-};
-
-/**
- * @brief Brute force neighbour search (used for comparison)
- */
-
-class SPH_API BruteForceSorter : public Sorter
-{
- double kappa;
-
-public:
- double hmax; //< expected max smothing length used to calculate the number of elements in the sgrid
- int method;
-
- BruteForceSorter(double _hmax, Kernel const &kernel);
- virtual ~BruteForceSorter();
- virtual void execute(std::vector<Particle *> &prts, int step = 0) override;
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-} // namespace sph
-
-#endif //WSORTER_H
diff --git a/sph/src/wTimeIntegration.cpp b/sph/src/wTimeIntegration.cpp
deleted file mode 100644
index 39723352b2b5bf455b76f995c64a917bfbaf0bea..0000000000000000000000000000000000000000
--- a/sph/src/wTimeIntegration.cpp
+++ /dev/null
@@ -1,32 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "wTimeIntegration.h"
-using namespace sph;
-
-TimeIntegration::TimeIntegration()
-{
-}
-
-TimeIntegration::~TimeIntegration()
-{
- std::cout << "~TimeIntegration()\n";
-}
-
-void TimeIntegration::write(std::ostream &out) const
-{
- out << "sph::TimeIntegration";
-}
diff --git a/sph/src/wTimeIntegration.h b/sph/src/wTimeIntegration.h
deleted file mode 100644
index a07b7cfd7176ac2392659e4f36c94ebf5893be2a..0000000000000000000000000000000000000000
--- a/sph/src/wTimeIntegration.h
+++ /dev/null
@@ -1,46 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#ifndef WTIMEINTEGRATION_H
-#define WTIMEINTEGRATION_H
-
-#include "sph.h"
-#include "wObject.h"
-#include <iostream>
-#include <vector>
-#include <memory>
-
-namespace sph
-{
-
-/**
- * @brief main sph TimeIntegration object
- */
-
-class SPH_API TimeIntegration : public fwk::wObject
-{
-public:
- TimeIntegration();
- virtual ~TimeIntegration();
-
-#ifndef SWIG
- virtual void write(std::ostream &out) const override;
-#endif
-};
-
-} // namespace sph
-
-#endif //WTIMEINTEGRATION_H
diff --git a/sph/tests/eigentest.py b/sph/tests/eigentest.py
deleted file mode 100644
index bf5f33b8a13e0a953cd25add6e120e882bee9a9f..0000000000000000000000000000000000000000
--- a/sph/tests/eigentest.py
+++ /dev/null
@@ -1,22 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-import sph
-
-
-sph.EigenTest().test1()
diff --git a/sph/tests/neighbours.py b/sph/tests/neighbours.py
deleted file mode 100644
index c95e6399e5ec466d381873162107f22903873d33..0000000000000000000000000000000000000000
--- a/sph/tests/neighbours.py
+++ /dev/null
@@ -1,92 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-import sph
-from tbox import Vector3d
-from fwk import Timer
-
-
-def saveNeighbours(pbl, fname):
- print('saving results to %s...' % fname)
- f = open(fname,'w')
- for p in pbl.prts:
- nos=[]
- for pn in p.neighs:
- nos.append(pn.no)
- f.write( ' p#%d (%d): ' % (p.no, len(p.neighs)))
- for no in sorted(nos):
- f.write( '%d,' % no)
- f.write('\n')
- f.close()
-
-
-def printLouisGrid(sorter):
- for i in range(sorter.size()):
- for j in range(sorter.size()):
- for k in range(sorter.size()):
- prts = sorter.getParticles(i,j,k)
- if len(prts)!=0:
- print("(i,l,k)=(%d,%d,%d) nbprts=%d" % (i,j,k,len(prts)))
-#
-
-
-boxL = 2.
-Lfloor = 0.7
-Lwater = 0.5
-sep = 0.015
-rho = 1000.0
-h0 = 1.2*sep
-
-
-pbl = sph.Problem(h0)
-
-# generate a set of particles
-
-import sph.genp as genp
-no=0
-
-cube = genp.Cube( o=(((boxL-Lwater)/2),((boxL-Lwater)/2), ((boxL)/2)+0.5), L=(Lwater,Lwater,Lwater), s=sep)
-for p in cube.generate():
- no+=1
- pbl.addMobileP(no, sph.Dofs(Vector3d(p[0], p[1], p[2]), Vector3d(0.0, 0.0, 0.0), rho))
-
-print(pbl)
-
-# kernel
-
-kernel = sph.CubicSplineKernel()
-print(kernel)
-
-
-# sorter
-
-timer = Timer(); timer.start()
-sorter = sph.LouisSorter(boxL, h0, kernel)
-print(sorter)
-sorter.execute(pbl.prts)
-timer.stop(); print(timer)
-#saveNeighbours(pbl, 'louissorter.txt')
-
-timer = Timer(); timer.start()
-sorter = sph.BruteForceSorter(h0, kernel)
-sorter.method = 2
-sorter.nthreads = 6
-print(sorter)
-sorter.execute(pbl.prts)
-timer.stop(); print(timer)
-#saveNeighbours(pbl, 'bruteforcesorter.txt')
diff --git a/sph/tests/sandbox.py b/sph/tests/sandbox.py
deleted file mode 100644
index da1853ea0d1fd91ee7dc775aadfaaa799d650459..0000000000000000000000000000000000000000
--- a/sph/tests/sandbox.py
+++ /dev/null
@@ -1,116 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-import sph
-from tbox import Vector3d
-from fwk import Timer
-
-
-# test dofs --
-d1 = sph.Dofs(Vector3d(0,0,0), Vector3d(0,0,0), 0.0)
-d2 = sph.Dofs(Vector3d(1,1,1), Vector3d(0,0,0), 0.0)
-print("dof1:", d1)
-print("dof2:", d2)
-
-#
-
-
-boxL = 2.
-Lfloor = 0.7
-Lwater = 0.5
-sep = 0.05
-rho = 1000.0
-h0 = 1.2*sep
-
-
-pbl = sph.Problem(h0)
-print(pbl)
-
-# generate a set of particles
-
-import sph.genp as genp
-
-no=0
-# mobile particles
-cube = genp.Cube( o=(((boxL-Lwater)/2),((boxL-Lwater)/2), ((boxL)/2)+0.5), L=(Lwater,Lwater,Lwater), s=sep)
-for p in cube.generate():
- no+=1
- pbl.addMobileP(no, sph.Dofs(Vector3d(p[0], p[1], p[2]), Vector3d(0.0, 0.0, 0.0), rho))
-
-
-# fixed particles
-
-plane = genp.Cube( o=(((boxL-Lfloor)/2),((boxL-Lfloor)/2), (boxL/2)), L=(Lfloor,Lfloor,sep), s=sep)
-for p in plane.generate():
- no+=1
- pbl.addMobileP(no, sph.Dofs(Vector3d(p[0], p[1], p[2]), Vector3d(0.0, 0.0, 0.0), rho))
-
-plane = genp.Cube( o=(0,0,0), L=(boxL,boxL,sep), s=sep)
-for p in plane.generate():
- no+=1
- pbl.addMobileP(no, sph.Dofs(Vector3d(p[0], p[1], p[2]), Vector3d(0.0, 0.0, 0.0), rho))
-
-print(pbl)
-
-# kernel
-
-kernel = sph.CubicSplineKernel()
-print(kernel)
-
-
-# sorter
-
-timer = Timer(); timer.start()
-sorter = sph.LouisSorter(boxL, h0, kernel)
-sorter.execute(pbl.prts)
-timer.stop(); print(timer)
-
-print(sorter)
-
-if 0:
- for i in range(sorter.size()):
- for j in range(sorter.size()):
- for k in range(sorter.size()):
- prts = sorter.getParticles(i,j,k)
- if len(prts)!=0:
- print("(i,l,k)=(%d,%d,%d) nbprts=%d" % (i,j,k,len(prts)))
-
-def saveNeighbours(pbl, fname):
- f = open(fname,'w')
- for p in pbl.prts:
- nos=[]
- for pn in p.neighs:
- nos.append(pn.no)
- f.write( ' p#%d (%d): ' % (p.no, len(p.neighs)))
- for no in sorted(nos):
- f.write( '%d,' % no)
- f.write('\n')
- f.close()
- print('%s saved to disk' % fname)
-
-saveNeighbours(pbl, 'louissorter.txt')
-
-timer = Timer(); timer.start()
-sorter = sph.BruteForceSorter(h0, kernel)
-sorter.execute(pbl.prts)
-timer.stop(); print(timer)
-
-
-print(sorter)
-
-saveNeighbours(pbl, 'bruteforcesorter.txt')