diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index c8ad725..8d190b4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -21,7 +21,7 @@ jobs: strategy: fail-fast: false matrix: - os: [macos-latest, ubuntu-latest] + os: [ubuntu-latest] steps: - uses: actions/checkout@master with: diff --git a/examples/partitioned-heat-conduction/README.md b/examples/partitioned-heat-conduction/README.md index 56378ac..20c908c 100644 --- a/examples/partitioned-heat-conduction/README.md +++ b/examples/partitioned-heat-conduction/README.md @@ -24,15 +24,9 @@ In the G+Smo build folder, build the tutorial file. ``` make -j <#threads> +make install ``` -Go to the gismo-executable folder and link the compiled executable to the gismo_executable. - -``` -cd gismo-executable -chmod +x create_symlink.sh -./create_symlink.sh -``` You can find the corresponding `run.sh` script for running the case in the folders corresponding to the participant you want to use: diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/run.sh b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/run.sh index d201a7f..6748129 100755 --- a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/run.sh +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/run.sh @@ -1,4 +1,4 @@ #!/bin/bash set -e -u -../gismo-executable/gismo-executable-dirichlet -c ../precice-config.xml +partitioned-heat-conduction-IGA-Dirichlet -c ../precice-config.xml diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/gismo-executable/creat_symlink.sh b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/gismo-executable/creat_symlink.sh deleted file mode 100755 index 8c1c480..0000000 --- a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/gismo-executable/creat_symlink.sh +++ /dev/null @@ -1,48 +0,0 @@ -#!/bin/bash - -# Path to the source file (relative to the script's directory) -SOURCE_FLUID="../../../../../../build/bin/partitioned-heat-conduction-IGA-Dirichlet" -SOURCE_SOLID="../../../../../../build/bin/partitioned-heat-conduction-IGA-Neumann" - -# Path to the symbolic link -LINK_FLUID="./gismo-executable-dirichlet" -LINK_SOLID="./gismo-executable-neumann" - -# Check if the source file exists -if [ -e "$SOURCE_FLUID" ]; then - # Check if the symbolic link already exists - if [ -L "$LINK_FLUID" ]; then - echo "Symbolic link already exists. Updating the link." - rm "$LINK_FLUID" # Remove the existing symbolic link - elif [ -e "$LINK_FLUID" ]; then - echo "A file or directory named '$LINK_FLUID' already exists. Please remove it manually." - exit 1 - fi - - # Create the symbolic link - ln -s "$SOURCE_FLUID" "$LINK_FLUID" - echo "Symbolic link created: $LINK_FLUID -> $SOURCE_FLUID" -else - echo "Source file '$SOURCE_FLUID' does not exist. Please check the path." - exit 1 -fi - - -# Check if the source file exists -if [ -e "$SOURCE_SOLID" ]; then - # Check if the symbolic link already exists - if [ -L "$LINK_SOLID" ]; then - echo "Symbolic link already exists. Updating the link." - rm "$LINK_SOLID" # Remove the existing symbolic link - elif [ -e "$LINK_SOLID" ]; then - echo "A file or directory named '$LINK_SOLID' already exists. Please remove it manually." - exit 1 - fi - - # Create the symbolic link - ln -s "$SOURCE_SOLID" "$LINK_SOLID" - echo "Symbolic link created: $LINK_SOLID -> $SOURCE_SOLID" -else - echo "Source file '$SOURCE_SOLID' does not exist. Please check the path." - exit 1 -fi diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/run.sh b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/run.sh index 2a424ec..0c643af 100755 --- a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/run.sh +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/run.sh @@ -1,4 +1,4 @@ #!/bin/bash set -e -u -../gismo-executable/gismo-executable-neumann -c ../precice-config.xml +partitioned-heat-conduction-IGA-Neumann -c ../precice-config.xml diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/dirichlet-gismo/run.sh b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/dirichlet-gismo/run.sh index 2612a8c..b381cf8 100755 --- a/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/dirichlet-gismo/run.sh +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/dirichlet-gismo/run.sh @@ -1,4 +1,4 @@ #!/bin/bash set -e -u -../gismo-executable/gismo-executable -c ../precice-config.xml --plot -s 0 +partitioned-heat-conduction -c ../precice-config.xml --plot -s 0 diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/gismo-executable/creat_symlink.sh b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/gismo-executable/creat_symlink.sh deleted file mode 100755 index 5fbe541..0000000 --- a/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/gismo-executable/creat_symlink.sh +++ /dev/null @@ -1,26 +0,0 @@ -#!/bin/bash - -# Path to the source file (relative to the script's directory) -SOURCE="../../../../../../build/bin/partitioned-heat-conduction" - -# Path to the symbolic link -LINK="./gismo-executable" - -# Check if the source file exists -if [ -e "$SOURCE" ]; then - # Check if the symbolic link already exists - if [ -L "$LINK" ]; then - echo "Symbolic link already exists. Updating the link." - rm "$LINK" # Remove the existing symbolic link - elif [ -e "$LINK" ]; then - echo "A file or directory named '$LINK' already exists. Please remove it manually." - exit 1 - fi - - # Create the symbolic link - ln -s "$SOURCE" "$LINK" - echo "Symbolic link created: $LINK -> $SOURCE" -else - echo "Source file '$SOURCE' does not exist. Please check the path." - exit 1 -fi diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/neumann-gismo/run.sh b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/neumann-gismo/run.sh index 6e6bddd..6ecf01e 100755 --- a/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/neumann-gismo/run.sh +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/neumann-gismo/run.sh @@ -1,4 +1,4 @@ #!/bin/bash set -e -u -../gismo-executable/gismo-executable -c ../precice-config.xml --plot -s 1 +partitioned-heat-conduction -c ../precice-config.xml --plot -s 1 diff --git a/examples/perpendicular-flap-IGA-IGA-fluid.cpp b/examples/perpendicular-flap-IGA-IGA-fluid.cpp deleted file mode 100644 index 67b2ab2..0000000 --- a/examples/perpendicular-flap-IGA-IGA-fluid.cpp +++ /dev/null @@ -1,282 +0,0 @@ -/** @file test-communication-IGA-IGA-fluid.cpp - * - * @brief test communicate spline through preCICE. To see if the spline force function can be reconstructed - * - * This file is a part of G+Smo library. - * - * This Source Code Form is subject to the terms of the Mozilla Public - * License, v. 2.0. If a copy of the MPL was not distributed with this - * file, You can obtain one at http://mozilla.org/MPL/2.0/. - * - Author(s): J.Li (TU Delft, 2023-...) - * - **/ - -#include -#include -#include -#include -// #include - -#include -#include - -#ifdef gsStructuralAnalysis_ENABLED -#include -#include -#include -#include -#include -#include -#include -#endif -using namespace gismo; - -int main(int argc, char *argv[]) -{ - - //! [Parse command line] - bool plot = false; - bool write = false; - bool get_readTime = false; - bool get_writeTime = false; - index_t plotmod = 1; - index_t loadCase = 0; - std::string precice_config("../precice_config.xml"); - - gsCmdLine cmd("Coupled heat equation using PreCICE."); - cmd.addString( "c", "config", "PreCICE config file", precice_config ); - cmd.addInt("l","loadCase", "Load case: 0=constant load, 1='spring' load", loadCase); - cmd.addSwitch("write", "Create a file with point data", write); - cmd.addSwitch("readTime", "Get the read time", get_readTime); - cmd.addSwitch("writeTime", "Get the write time", get_writeTime); - try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } - - GISMO_ASSERT(gsFileManager::fileExists(precice_config),"No precice config file has been defined"); - //! [Read input file] - std::string participantName = "Fluid"; - gsPreCICE participant(participantName, precice_config); - - - /* - * Data initialization - * - * This participant receives mesh information (knot vector, control points) from the solid, - * and it creates its own mesh based on the information. - * And writes pressure, reads displacements from the solid participant. - * The follow meshes and data are made available: - * - * - Meshes: - * + GeometryControlPointMesh: This mesh contains the control points as mesh vertices. - * + GeometryKnotMesh: This mesh contains the knots as mesh vertices. - * + ForceKnotMesh: This mesh contains the knots as mesh vertices for force mesh. - * + ForceControlPointMesh: This mesh contains the control points as mesh vertices for force mesh. - * - * - * - Data: - * + GeometryControlPointData: This data is defined on the ControlPointMesh and stores the displacement of the control points - * + ForceControlPointData: This data is defined on the ForceControlPointMesh and stores the change of pressure/stress - */ - - - // Meshes - std::string GeometryControlPointMesh = "GeometryControlPointMesh"; - std::string GeometryKnotMesh = "GeometryKnotMesh"; - - - std::string ForceKnotMesh = "ForceKnotMesh"; - std::string ForceControlPointMesh = "ForceControlPointMesh"; - - - // Data - std::string GeometryControlPointData = "GeometryControlPointData"; - std::string ForceControlPointData = "ForceControlPointData"; - // std::string GeometryKnotData = "GeometryKnotData"; - // std::string ForceKnotData = "ForceKnotMeshData"; - - gsMatrix<> bbox(3,2); - bbox << -1e300, 1e300, // X dimension limits - -1e300, 1e300, // Y dimension limits - -1e300, 1e300; // Z dimension limits - bbox.transposeInPlace(); - bbox.resize(1,bbox.rows()*bbox.cols()); - - // Step 2: Add force mesh - - // Step 2.1: Define force knot mesh - gsVector forceKnotIDs; - gsKnotVector<> kv(0,1,10,3,1); //begin, end, # interiors, mult, by default =1 - gsTensorBSplineBasis<2, real_t> forceBasis(kv,kv); - gsMatrix<> forceControlPoints(forceBasis.size(), 3); - forceControlPoints.setZero(); - forceControlPoints.col(2).setConstant(-5e3); - forceControlPoints.transposeInPlace(); - - gsMatrix<> forceKnots = knotsToMatrix(forceBasis); - - - - - // gsMultiBasis<> bases(forceMesh); - // gsMatrix<> forceKnots = knotsToMatrix(bases.basis(0)); - // gsDebugVar(forceKnots.dim()); - // - - - - // Regenerate the geometry in another solver - // forceMesh.addPatch(give(forceBasis.makeGeometry(forceControlPoints))); - - // Step 2.1: Define force control point mesh - gsVector forceControlPointIDs; - participant.addMesh(ForceKnotMesh,forceKnots,forceKnotIDs); - participant.addMesh(ForceControlPointMesh,forceControlPoints,forceControlPointIDs); - participant.setMeshAccessRegion(GeometryControlPointMesh, bbox); - - real_t precice_dt = participant.initialize(); - - // Step 1: Collect geometry from PreCICE - gsVector geometryKnotIDs; - gsMatrix<> geometryKnots; - participant.getMeshVertexIDsAndCoordinates(GeometryKnotMesh, geometryKnotIDs, geometryKnots); - gsDebugVar(geometryKnots); - - - gsBasis<> * geometryKnotBasis = knotMatrixToBasis(geometryKnots).get(); - - gsVector geometryControlPointIDs; - gsMatrix<> geometryControlPoints; - participant.getMeshVertexIDsAndCoordinates(GeometryControlPointMesh, geometryControlPointIDs, geometryControlPoints); - gsDebugVar(geometryControlPoints.dim()); - - gsMultiPatch<> mp, deformation; - gsDebugVar(*geometryKnotBasis); - gsDebugVar(geometryKnotBasis->size()); - gsDebugVar(geometryControlPoints.rows()); - mp.addPatch(give(geometryKnotBasis->makeGeometry(geometryControlPoints.transpose()))); - - deformation = mp; - - - real_t t=0; - real_t dt = precice_dt; - real_t t_read = 0; - real_t t_write = 0; - - index_t timestep = 0; - - // Define the solution collection for Paraview - gsParaviewCollection collection("./output/solution"); - - gsMatrix<> points(2,1); - points.col(0)<<0.5,1; - - gsStructuralAnalysisOutput writer("./output/pointData.csv",points); - writer.init({"x","y","z"},{"time"}); - - // Time integration loop - while(participant.isCouplingOngoing()) - { - - // Read control point displacements - participant.readData(GeometryControlPointMesh, GeometryControlPointData, geometryControlPointIDs, geometryControlPoints); - deformation.patch(0).coefs() = geometryControlPoints.transpose(); - if (get_readTime) - t_read += participant.readTime(); - - /* - * Two projection approaches: - * - gsQuasiInterpolate::localIntpl(forceBasis, deformation.patch(0), coefs); - * - gsL2Projection::projectFunction(forceBasis, deformation, mp, coefs) - * - * Check if forceBasis + coefs gives the same as deformation.patch(0) - */ - if (loadCase == 0) - { - gsInfo << "Load case 0: Constant load\n"; - } - else if (loadCase == 1) - { - gsInfo << "Load case 1: Spring load\n"; - if(timestep > 0) - { - gsMatrix<> coefs; - gsL2Projection::projectFunction(forceBasis, deformation, mp, coefs); - coefs.resize(forceControlPoints.rows(),forceControlPoints.cols()); - forceControlPoints.row(2) = -coefs.row(2); - } - } - else - { - GISMO_ERROR("Load case "< pointDataMatrix; - gsMatrix<> otherDataMatrix(1,1); - if (participant.requiresWritingCheckpoint()) - { - gsInfo<<"Reading Checkpoint\n"; - } - - // do the coupling - precice_dt = participant.advance(dt); - - dt = std::min(precice_dt, dt); - - if (participant.requiresReadingCheckpoint()) - { - gsInfo<<"Writing Checkpoint \n"; - } - else - { - t += dt; - timestep++; - - gsField<> solution(mp,deformation); - if (timestep % plotmod==0 && plot) - { - // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here - std::string fileName = "./output/solution" + util::to_string(timestep); - gsWriteParaview<>(solution, fileName, 500); - fileName = "solution" + util::to_string(timestep) + "0"; - collection.addTimestep(fileName,t,".vts"); - } - if (write) - { - solution.patch(0).eval_into(points,pointDataMatrix); - otherDataMatrix< -#include -#include -#include -// #include - -#include -#include - - -#ifdef gsStructuralAnalysis_ENABLED -#include -#include -#include -#include -#include -#include -#include -#endif - -#ifdef gsKLShell_ENABLED -#include -#include -#include -#include -#include -#endif - - -using namespace gismo; - -// template -// T ProjectForceMesh( const gsMatrix<> & forceControlPoints, -// const gsMultiBasis<> & forceBasis, -// const gsMultiBasis<> & geometryBasis, -// const index_t & targetDim, -// gsMatrix & results) -// { -// gsExprAssembler<> A(1,1); - -// gsMatrix solVector; - -// // Set the discretization space -// space u = A.getSpace(forceBasis, targetDim); - -// solution sol = A.getSolution(u, solVector); -// auto f = A.getCoordinates(u); - -// u.setup(0); -// A.initSystem(); - - - -// gsMatrix<> matrixForce, matrixGeometry; -// } - - -int main(int argc, char *argv[]) -{ - bool plot = false; - bool get_readTime = false; - bool get_writeTime = false; - index_t plotmod = 1; - index_t numRefine = 1; - index_t numElevate = 0; - int method = 3; // 1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson, 6 RK4 - - real_t rho = 3000; - real_t E = 4e6; - real_t nu = 0.3; - - //! [Parse command line] - std::string precice_config("../precice_config.xml"); - - gsCmdLine cmd("Coupled heat equation using PreCICE."); - cmd.addInt( "e", "degreeElevation", - "Number of degree elevation steps to perform before solving (0: equalize degree in all directions)", numElevate ); - cmd.addInt( "r", "uniformRefine", "Number of Uniform h-refinement loops", numRefine ); - cmd.addString( "c", "config", "PreCICE config file", precice_config ); - cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot); - cmd.addInt("m", "method","1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson",method); - cmd.addSwitch("readTime", "Get the read time", get_readTime); - cmd.addSwitch("writeTime", "Get the write time", get_writeTime); - try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } - - GISMO_ASSERT(gsFileManager::fileExists(precice_config),"No precice config file has been defined"); - //! [Read input file] - std::string participantName = "Solid"; - gsPreCICE participant(participantName, precice_config); - - //Generate domain - gsMultiPatch<> patches; - patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,0.0,0.5,1.0)); - - //Embed the 2D geometry to 3D - gsMultiPatch<> solutions; - patches.addAutoBoundaries(); - patches.embed(3); - - // Create bases - // p-refine - if (numElevate!=0) - patches.degreeElevate(numElevate); - - // h-refine - for (int r =0; r < numRefine; ++r) - patches.uniformRefine(); - - // Create bases - gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) - - gsInfo << "Patches: "<< patches.nPatches() <<", degree: "<< bases.minCwiseDegree() <<"\n"; - - - /* - * Data initialization - * - * This participant receives mesh information (knot vector, control points) from the solid, - * and it creates its own mesh based on the information. - * And writes pressure, reads displacements from the solid participant. - * The follow meshes and data are made available: - * - * - Meshes: - * + GeometryControlPointMesh: This mesh contains the control points as mesh vertices. - * + GeometryKnotMesh: This mesh contains the knots as mesh vertices. - * + ForceKnotMesh: This mesh contains the knots as mesh vertices for force mesh. - * + ForceControlPointMesh: This mesh contains the control points as mesh vertices for force mesh. - * - * - * - Data: - * + GeometryControlPointData: This data is defined on the GeometryControlPointMesh and stores the displacement of the control points - * + GeometryKnotData - * + ForceControlPointData: This data is defined on the ForceControlPointMesh and stores the change of pressure/stress - * + ForceKnotData - */ - - std::string GeometryKnotMesh = "GeometryKnotMesh"; - std::string GeometryControlPointMesh = "GeometryControlPointMesh"; - std::string ForceKnotMesh = "ForceKnotMesh"; - std::string ForceControlPointMesh = "ForceControlPointMesh"; - - // std::string StressData = "StressData"; - std::string GeometryControlPointData = "GeometryControlPointData"; - // std::string GeometryKnotData = "GeometryKnotData"; - // std::string ForceKnotData = "ForceKnotMeshData"; - std::string ForceControlPointData = "ForceControlPointData"; - - gsMatrix<> bbox(3,2); - bbox << -1e300, 1e300, // X dimension limits - -1e300, 1e300, // Y dimension limits - -1e300, 1e300; // Z dimension limits - bbox.transposeInPlace(); - bbox.resize(1,bbox.rows()*bbox.cols()); - participant.setMeshAccessRegion(ForceControlPointMesh, bbox); - - // Setup geometry control point mesh - gsVector geometryControlPointIDs; - gsMatrix<> geometryControlPoints = patches.patch(0).coefs().transpose(); - participant.addMesh(GeometryControlPointMesh, geometryControlPoints, geometryControlPointIDs); - - - // Setup the geometry knot mesh - gsMatrix<> geometryKnots = knotsToMatrix(bases.basis(0)); - - - participant.addMesh(GeometryKnotMesh,geometryKnots); - - real_t precice_dt = participant.initialize(); - - //G[et the force mesh from the coupling interface - gsVector forceKnotIDs; - gsMatrix<> forceKnots; - participant.getMeshVertexIDsAndCoordinates(ForceKnotMesh,forceKnotIDs,forceKnots); - - gsBasis<> * basis = knotMatrixToBasis(forceKnots).get(); - - gsVector forceControlPointIDs; - gsMatrix<> forceControlPoints; - participant.getMeshVertexIDsAndCoordinates(ForceControlPointMesh, forceControlPointIDs,forceControlPoints); - - - // // Step 2: Regenerate the geometry - gsMultiPatch<> forceMesh; //Geometry object belongs to gsFunctionSet - forceMesh.addPatch(give(basis->makeGeometry(forceControlPoints.transpose()))); - - - // Define boundary condition for solid mesh - gsBoundaryConditions<> bcInfo; - //Bottom Side - bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, nullptr, -1); - bcInfo.addCondition(0, boundary::south, condition_type::clamped, nullptr, 2); - - // Assign geometry map - bcInfo.setGeoMap(patches); - - // Surface force equals to the force mesh - -//------------------------------------------ gsKLShell Module ------------------------------------------------ - //Setup the material properties - gsFunctionExpr<> E_modulus(std::to_string(E),3); - gsFunctionExpr<> PoissonRatio(std::to_string(nu),3); - gsFunctionExpr<> Density(std::to_string(rho),3); - - gsOptionList options; - - - //Define thickness - real_t thickness = 0.1; - - gsFunctionExpr<> t(std::to_string(thickness),3); - - std::vector*> parameters(2); - parameters[0] = &E_modulus; - parameters[1] = &PoissonRatio; - - gsMaterialMatrixBase* materialMatrix; - options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",0); - options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",1); - materialMatrix = getMaterialMatrix<3, real_t>(patches, t, parameters, Density, options); - - // REMOVE IT LATER - - // forceMesh.embed(3); - - gsThinShellAssembler<3, real_t, true> assembler(patches, bases, bcInfo, forceMesh, materialMatrix); - gsOptionList assemblerOptions = options.wrapIntoGroup("Assembler"); - - assembler.assemble(); - // forceMesh.patch(0).coefs().setZero(); - // assembler.assemble(); - // gsDebugVar(assembler.rhs()); - - assembler.setOptions(assemblerOptions); - - index_t timestep = 0; - index_t timestep_checkpoint = 0; - - - // Compute the mass matrix (since it is constant over time) - assembler.assembleMass(); - gsSparseMatrix<> M = assembler.massMatrix(); - assembler.assemble(); - gsSparseMatrix<> K = assembler.matrix(); - - - // Define the solution collection for Paraview - std::string dirname = "./output"; - - gsFileManager::mkdir(dirname); - gsParaviewCollection collection(dirname + "/solution"); - - // Time step - real_t dt = precice_dt; - real_t t_read = 0; - real_t t_write = 0; - - gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&solutions](gsMatrix const &x, gsSparseMatrix & m) - { - // to do: add time dependency of forcing - // For the shell - ThinShellAssemblerStatus status; - assembler.constructSolution(x, solutions); - status = assembler.assembleMatrix(solutions); - m = assembler.matrix(); - return true; - }; - - - // Function for the Residual - gsStructuralAnalysisOps::TResidual_t Residual = [&assembler,&solutions](gsMatrix const &x, real_t t, gsVector & result) - { - //Add assemble vector JL - ThinShellAssemblerStatus status; - assembler.constructSolution(x,solutions); - status = assembler.assembleVector(solutions); - result = assembler.rhs(); - return true; - }; - - -gsSparseMatrix<> C = gsSparseMatrix<>(assembler.numDofs(),assembler.numDofs()); - gsStructuralAnalysisOps::Damping_t Damping = [&C](const gsVector &, gsSparseMatrix & m) { m = C; return true; }; - gsStructuralAnalysisOps::Mass_t Mass = [&M]( gsSparseMatrix & m) { m = M; return true; }; - - gsDynamicBase * timeIntegrator; - if (method==1) - timeIntegrator = new gsDynamicExplicitEuler(Mass,Damping,Jacobian,Residual); - else if (method==2) - timeIntegrator = new gsDynamicImplicitEuler(Mass,Damping,Jacobian,Residual); - else if (method==3) - timeIntegrator = new gsDynamicNewmark(Mass,Damping,Jacobian,Residual); - else if (method==4) - timeIntegrator = new gsDynamicBathe(Mass,Damping,Jacobian,Residual); - else if (method==5) - { - timeIntegrator = new gsDynamicWilson(Mass,Damping,Jacobian,Residual); - timeIntegrator->options().setReal("gamma",1.4); - } - else if (method==6) - timeIntegrator = new gsDynamicRK4(Mass,Damping,Jacobian,Residual); - else - GISMO_ERROR("Method "<options().setReal("DT",dt); - timeIntegrator->options().setReal("TolU",1e-3); - timeIntegrator->options().setSwitch("Verbose",true); - - - - // Project u_wall as initial condition (violates Dirichlet side on precice interface) - // RHS of the projection - gsMatrix<> solVector; - solVector.setZero(assembler.numDofs(),1); - - // Assemble the RHS - gsVector<> F = assembler.rhs(); - - gsVector<> F_checkpoint, U_checkpoint, V_checkpoint, A_checkpoint, U, V, A; - - F_checkpoint = F; - U_checkpoint = U = gsVector::Zero(assembler.numDofs(),1); - V_checkpoint = V = gsVector::Zero(assembler.numDofs(),1); - A_checkpoint = A = gsVector::Zero(assembler.numDofs(),1); - - - real_t time = 0; - // Plot initial solution - if (plot) - { - gsMultiPatch<> solution; - gsVector<> displacements = U; - solution = assembler.constructDisplacement(displacements); - solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here - gsField<> solField(patches,solution); - std::string fileName = dirname + "/solution" + util::to_string(timestep); - gsWriteParaview<>(solField, fileName, 500); - fileName = "solution" + util::to_string(timestep) + "0"; - collection.addTimestep(fileName,time,".vts"); - } - - gsMatrix<> points(2,1); - points.col(0)<<0.5,1; - - gsStructuralAnalysisOutput writer("./output/pointData.csv",points); - writer.init({"x","y","z"},{"time"}); // point1 - x, point1 - y, point1 - z, time - - gsMatrix<> pointDataMatrix; - gsMatrix<> otherDataMatrix(1,1); - - // gsDebugVar(dt); - - // Time integration loop - while (participant.isCouplingOngoing()) - { - if (participant.requiresWritingCheckpoint()) - { - U_checkpoint = U; - V_checkpoint = V; - A_checkpoint = A; - - gsInfo<<"Checkpoint written:\n"; - gsInfo<<"\t ||U|| = "<step(time,dt,U,V,A); - solVector = U; - - gsInfo<<"Finished\n"; - - // potentially adjust non-matching timestep sizes - dt = std::min(dt,precice_dt); - - gsMultiPatch<> solution; - gsVector<> displacements = U; - solution = assembler.constructDisplacement(displacements); - // write heat fluxes to interface - geometryControlPoints = solution.patch(0).coefs().transpose(); - participant.writeData(GeometryControlPointMesh,GeometryControlPointData,geometryControlPointIDs,geometryControlPoints); - - if (get_writeTime) - t_write +=participant.writeTime(); - - - // do the coupling - precice_dt =participant.advance(dt); - - if (participant.requiresReadingCheckpoint()) - { - U = U_checkpoint; - V = V_checkpoint; - A = A_checkpoint; - timestep = timestep_checkpoint; - } - else - { - // gsTimeIntegrator advances the time step - // advance variables - time += dt; - timestep++; - - gsField<> solField(patches,solution); - if (timestep % plotmod==0 && plot) - { - // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here - std::string fileName = dirname + "/solution" + util::to_string(timestep); - gsWriteParaview<>(solField, fileName, 500); - fileName = "solution" + util::to_string(timestep) + "0"; - collection.addTimestep(fileName,time,".vts"); - } - solution.patch(0).eval_into(points,pointDataMatrix); - otherDataMatrix< +#include +#include +#include +// #include + +#include +#include + +#ifdef gsStructuralAnalysis_ENABLED +#include +#include +#include +#include +#include +#include +#include +#endif + +using namespace gismo; + +int main(int argc, char *argv[]) +{ + //! [Parse command line] + bool plot = false; + index_t plotmod = 1; + index_t numRefine = 0; + index_t numElevate = 0; + std::string precice_config; + int method = 3; // 1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson, 6 RK4 + + gsCmdLine cmd("Flow over heated plate for PreCICE."); + cmd.addInt( "e", "degreeElevation", + "Number of degree elevation steps to perform before solving (0: equalize degree in all directions)", numElevate ); + cmd.addInt( "r", "uniformRefine", "Number of Uniform h-refinement loops", numRefine ); + cmd.addString( "c", "config", "PreCICE config file", precice_config ); + cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot); + cmd.addInt("m","plotmod", "Modulo for plotting, i.e. if plotmod==1, plots every timestep", plotmod); + cmd.addInt("M", "method","1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson, 6: RK4",method); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + //! [Read input file] + GISMO_ASSERT(gsFileManager::fileExists(precice_config),"No precice config file has been defined"); + + // Generate domain + gsMultiPatch<> patches; + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(-0.05,0.0,0.05,1.0)); + + // Create bases + gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) + + // Set degree + bases.setDegree( bases.maxCwiseDegree() + numElevate); + + // h-refine each basis + for (int r =0; r < numRefine; ++r) + bases.uniformRefine(); + numRefine = 0; + + gsInfo << "Patches: "<< patches.nPatches() <<", degree: "<< bases.minCwiseDegree() <<"\n"; + + real_t rho = 3000; + real_t E = 4e6; + real_t nu = 0.3; + + // Set the interface for the precice coupling + std::vector couplingInterfaces(3); + couplingInterfaces[0] = patchSide(0,boundary::east); + couplingInterfaces[1] = patchSide(0,boundary::north); + couplingInterfaces[2] = patchSide(0,boundary::west); + + /* + * Initialize the preCICE participant + * + * + */ + std::string participantName = "Solid"; + gsPreCICE participant(participantName, precice_config); + + /* + * Data initialization + * + * This participant manages the geometry. The follow meshes and data are made available: + * + * - Meshes: + * + KnotMesh This mesh contains the knots as mesh vertices + * + ControlPointMesh: This mesh contains the control points as mesh vertices + * + ForceMesh: This mesh contains the integration points as mesh vertices + * + * - Data: + * + ControlPointData: This data is defined on the ControlPointMesh and stores the displacement of the control points + * + ForceData: This data is defined on the ForceMesh and stores pressure/forces + */ + + std::string SolidMesh = "Solid-Mesh"; + std::string StressData = "Stress"; + std::string DisplacementData = "Displacement"; + + // Get the quadrature nodes on the coupling interface + gsOptionList quadOptions = gsExprAssembler<>::defaultOptions(); + + // Get the quadrature points + gsMatrix<> quad_uv = gsQuadrature::getAllNodes(bases.basis(0),quadOptions,couplingInterfaces); // Quadrature points in the parametric domain + gsMatrix<> quad_xy = patches.patch(0).eval(quad_uv); // Quadrature points in the physical domain + gsVector quad_xyIDs; // needed for writing + participant.addMesh(SolidMesh,quad_xy,quad_xyIDs); + + // Define precice interface + real_t precice_dt = participant.initialize(); + +// ---------------------------------------------------------------------------------------------- + + // Define boundary conditions + gsBoundaryConditions<> bcInfo; + // Dirichlet side + gsConstantFunction<> g_D(0,patches.geoDim()); + // Coupling side + // gsFunctionExpr<> g_C("1","0",patches.geoDim()); + gsPreCICEFunction g_C(&participant,SolidMesh,StressData,patches,patches.geoDim(),false); + // Add all BCs + // Coupling interface + bcInfo.addCondition(0, boundary::north, condition_type::neumann , &g_C,-1,true); + bcInfo.addCondition(0, boundary::east, condition_type::neumann , &g_C,-1,true); + bcInfo.addCondition(0, boundary::west, condition_type::neumann , &g_C,-1,true); + // Bottom side (prescribed temp) + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &g_D, 0); + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &g_D, 1); + // Assign geometry map + bcInfo.setGeoMap(patches); + +// ---------------------------------------------------------------------------------------------- + // source function, rhs + gsConstantFunction<> g(0.,0.,2); + + // source function, rhs + gsConstantFunction<> gravity(0.,0.,2); + + // creating mass assembler + gsMassAssembler massAssembler(patches,bases,bcInfo,gravity); + massAssembler.options().setReal("Density",rho); + massAssembler.assemble(); + + // creating stiffness assembler. + gsElasticityAssembler assembler(patches,bases,bcInfo,g); + assembler.options().setReal("YoungsModulus",E); + assembler.options().setReal("PoissonsRatio",nu); + assembler.assemble(); + + gsMatrix Minv; + gsSparseMatrix<> M = massAssembler.matrix(); + gsSparseMatrix<> K = assembler.matrix(); + gsSparseMatrix<> K_T; + + // Time step + real_t dt = precice_dt; + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + // RHS of the projection + gsMatrix<> solVector; + solVector.setZero(assembler.numDofs(),1); + + std::vector > fixedDofs = assembler.allFixedDofs(); + // Assemble the RHS + gsVector<> F = assembler.rhs(); + gsVector<> F_checkpoint, U_checkpoint, V_checkpoint, A_checkpoint, U, V, A; + + F_checkpoint = F; + U_checkpoint = U = gsVector::Zero(assembler.numDofs(),1); + V_checkpoint = V = gsVector::Zero(assembler.numDofs(),1); + A_checkpoint = A = gsVector::Zero(assembler.numDofs(),1); + + // Define the solution collection for Paraview + gsParaviewCollection collection("./output/solution"); + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&fixedDofs](gsMatrix const &x, gsSparseMatrix & m) { + // to do: add time dependency of forcing + assembler.assemble(x, fixedDofs); + m = assembler.matrix(); + return true; + }; + + // Function for the Residual + gsStructuralAnalysisOps::TResidual_t Residual = [&assembler,&fixedDofs](gsMatrix const &x, real_t t, gsVector & result) + { + assembler.assemble(x,fixedDofs); + result = assembler.rhs(); + return true; + }; + + + gsSparseMatrix<> C = gsSparseMatrix<>(assembler.numDofs(),assembler.numDofs()); + gsStructuralAnalysisOps::Damping_t Damping = [&C](const gsVector &, gsSparseMatrix & m) { m = C; return true; }; + gsStructuralAnalysisOps::Mass_t Mass = [&M]( gsSparseMatrix & m) { m = M; return true; }; + + gsDynamicBase * timeIntegrator; + if (method==1) + timeIntegrator = new gsDynamicExplicitEuler(Mass,Damping,Jacobian,Residual); + else if (method==2) + timeIntegrator = new gsDynamicImplicitEuler(Mass,Damping,Jacobian,Residual); + else if (method==3) + timeIntegrator = new gsDynamicNewmark(Mass,Damping,Jacobian,Residual); + else if (method==4) + timeIntegrator = new gsDynamicBathe(Mass,Damping,Jacobian,Residual); + else if (method==5) + { + timeIntegrator = new gsDynamicWilson(Mass,Damping,Jacobian,Residual); + timeIntegrator->options().setReal("gamma",1.4); + } + else if (method==6) + timeIntegrator = new gsDynamicRK4(Mass,Damping,Jacobian,Residual); + else + GISMO_ERROR("Method "<options().setReal("DT",dt); + timeIntegrator->options().setReal("TolU",1e-3); + timeIntegrator->options().setSwitch("Verbose",true); + + real_t time = 0; + + // Plot initial solution + if (plot) + { + gsMultiPatch<> solution; + assembler.constructSolution(solVector,fixedDofs,solution); + + gsField<> solField(patches,solution); + std::string fileName = "./output/solution" + util::to_string(timestep); + gsWriteParaview<>(solField, fileName, 500); + fileName = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileName,time,".vts"); + } + + gsMatrix<> points(2,1); + points.col(0)<<0.5,1; + + gsStructuralAnalysisOutput writer("pointData.csv",points); + writer.init({"x","y"},{"time"}); // point1 - x, point1 - y, time + + gsMatrix<> pointDataMatrix; + gsMatrix<> otherDataMatrix(1,1); + + // Time integration loop + while (participant.isCouplingOngoing()) + { + if (participant.requiresWritingCheckpoint()) + { + U_checkpoint = U; + V_checkpoint = V; + A_checkpoint = A; + + gsInfo<<"Checkpoint written:\n"; + gsInfo<<"\t ||U|| = "<step(time,dt,U,V,A); + solVector = U; + gsInfo<<"Finished\n"; + + // potentially adjust non-matching timestep sizes + dt = std::min(dt,precice_dt); + + + gsMultiPatch<> solution; + assembler.constructSolution(solVector,fixedDofs,solution); + // write heat fluxes to interface + gsMatrix<> result(patches.geoDim(),quad_uv.cols()); + solution.patch(0).eval_into(quad_uv,result); + participant.writeData(SolidMesh,DisplacementData,quad_xyIDs,result); + + // do the coupling + precice_dt =participant.advance(dt); + + if (participant.requiresReadingCheckpoint()) + { + U = U_checkpoint; + V = V_checkpoint; + A = A_checkpoint; + timestep = timestep_checkpoint; + } + else + { + // gsTimeIntegrator advances + // advance variables + time += dt; + timestep++; + + gsField<> solField(patches,solution); + if (timestep % plotmod==0 && plot) // Generate Paraview output for visualization + { + // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here + std::string fileName = "./output/solution" + util::to_string(timestep); + gsWriteParaview<>(solField, fileName, 500); + fileName = "solution" + util::to_string(timestep) + "0"; + collection.addTimestep(fileName,time,".vts"); + } + + solution.patch(0).eval_into(points,pointDataMatrix); + otherDataMatrix< -j <#threads> +make install +``` + + +You can find the corresponding `run.sh` script for running the case in the folders corresponding to the participant you want to use: + +for the fluid participant with OpenFOAM, +```bash +cd fluid-openfoam +OpenFOAM // Enter the OpenFOAM environment +./run.sh +``` + +and + +```bash +cd solid-gismo +./run.sh +``` + +## Post-processing +The results of this tutorial are comparable to the simulation results communicated with force under the perpendicular-flap tutorials. +![G+Smo stress](https://github.com/Crazy-Rich-Meghan/tutorials/blob/perpendicular-flap-gismo-elasticity-stress/perpendicular-flap-stress/images/tutorials-perpendicular-flap-displacement-openfoam-gismo-elasticity.png) + +Additionally, the mesh convergence study data is available under the `images/data` directory. A sample plot illustrating the convergence of x-displacement over time is shown below: +![G+Smo converfence](https://github.com/Crazy-Rich-Meghan/tutorials/blob/perpendicular-flap-gismo-elasticity-stress/perpendicular-flap-stress/images/x_displacement_vs_time.png) diff --git a/examples/perpendicular-flap/fluid-openfoam/0/U b/examples/perpendicular-flap/fluid-openfoam/0/U new file mode 100644 index 0000000..5ad1913 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/0/U @@ -0,0 +1,41 @@ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + object U; +} + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (10 0 0); + +boundaryField +{ + inlet + { + type fixedValue; + value $internalField; + } + outlet + { + type zeroGradient; + } + flap + { + type movingWallVelocity; + value uniform (0 0 0); + } + upperWall + { + type noSlip; + } + lowerWall + { + type noSlip; + } + frontAndBack + { + type empty; + } +} diff --git a/examples/perpendicular-flap/fluid-openfoam/0/p b/examples/perpendicular-flap/fluid-openfoam/0/p new file mode 100644 index 0000000..9ab4557 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/0/p @@ -0,0 +1,45 @@ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p; +} + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + inlet + { + type zeroGradient; + } + + outlet + { + type fixedValue; + value uniform 0; + } + + flap + { + type zeroGradient; + } + + upperWall + { + type zeroGradient; + } + + lowerWall + { + type zeroGradient; + } + + frontAndBack + { + type empty; + } +} diff --git a/examples/perpendicular-flap/fluid-openfoam/0/phi b/examples/perpendicular-flap/fluid-openfoam/0/phi new file mode 100644 index 0000000..9221e03 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/0/phi @@ -0,0 +1,44 @@ +FoamFile +{ + version 2.0; + format ascii; + class surfaceScalarField; + object phi; +} + +dimensions [0 3 -1 0 0 0 0]; + +internalField uniform 0; +boundaryField +{ + inlet + { + type calculated; + value $internalField; + } + outlet + { + type calculated; + value $internalField; + } + flap + { + type calculated; + value uniform 0; + } + upperWall + { + type calculated; + value uniform 0; + } + lowerWall + { + type calculated; + value uniform 0; + } + frontAndBack + { + type empty; + value nonuniform 0; + } +} diff --git a/examples/perpendicular-flap/fluid-openfoam/0/pointDisplacement b/examples/perpendicular-flap/fluid-openfoam/0/pointDisplacement new file mode 100644 index 0000000..4ac8684 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/0/pointDisplacement @@ -0,0 +1,47 @@ +FoamFile +{ + version 2.0; + format ascii; + class pointVectorField; + object pointDisplacement; +} + +dimensions [0 1 0 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + inlet + { + type fixedValue; + value uniform (0 0 0); + } + + outlet + { + type fixedValue; + value uniform (0 0 0); + } + + flap + { + type fixedValue; + value $internalField; + } + + upperWall + { + type slip; + } + + lowerWall + { + type slip; + } + + frontAndBack + { + type empty; + } +} diff --git a/examples/perpendicular-flap/fluid-openfoam/clean.sh b/examples/perpendicular-flap/fluid-openfoam/clean.sh new file mode 100755 index 0000000..b64fc51 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/clean.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_openfoam . diff --git a/examples/perpendicular-flap/fluid-openfoam/constant/dynamicMeshDict b/examples/perpendicular-flap/fluid-openfoam/constant/dynamicMeshDict new file mode 100644 index 0000000..2a3c1a2 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/constant/dynamicMeshDict @@ -0,0 +1,18 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object dynamicMeshDict; +} + +dynamicFvMesh dynamicMotionSolverFvMesh; + +motionSolverLibs ("libfvMotionSolvers.so"); + +solver displacementLaplacian; +// OpenFOAM9 or newer: rename "solver" to "motionSolver" + +displacementLaplacianCoeffs { + diffusivity quadratic inverseDistance (flap); +} diff --git a/examples/perpendicular-flap/fluid-openfoam/constant/transportProperties b/examples/perpendicular-flap/fluid-openfoam/constant/transportProperties new file mode 100644 index 0000000..5383ada --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/constant/transportProperties @@ -0,0 +1,11 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object transportProperties; +} + +transportModel Newtonian; + +nu nu [ 0 2 -1 0 0 0 0 ] 1; diff --git a/examples/perpendicular-flap/fluid-openfoam/constant/turbulenceProperties b/examples/perpendicular-flap/fluid-openfoam/constant/turbulenceProperties new file mode 100644 index 0000000..592f6d5 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/constant/turbulenceProperties @@ -0,0 +1,9 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object turbulenceProperties; +} + +simulationType laminar; diff --git a/examples/perpendicular-flap/fluid-openfoam/run.sh b/examples/perpendicular-flap/fluid-openfoam/run.sh new file mode 100755 index 0000000..8f55fbf --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/run.sh @@ -0,0 +1,12 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + +blockMesh + +../../tools/run-openfoam.sh "$@" +. ../../tools/openfoam-remove-empty-dirs.sh && openfoam_remove_empty_dirs + +close_log diff --git a/examples/perpendicular-flap/fluid-openfoam/system/blockMeshDict b/examples/perpendicular-flap/fluid-openfoam/system/blockMeshDict new file mode 100644 index 0000000..875f422 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/system/blockMeshDict @@ -0,0 +1,145 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} + +x0 -3.; +x1 -0.05; +x2 0.05; +x3 3.; + +y0 0.; +y1 1.; +y2 4.; + +z0 0; +z1 1; + +vertices +( + ($x0 $y0 $z0 ) // 0 + ($x1 $y0 $z0 ) // 1 + ($x2 $y0 $z0 ) // 2 + ($x3 $y0 $z0 ) // 3 + ($x0 $y1 $z0 ) // 4 + ($x1 $y1 $z0 ) // 5 + ($x2 $y1 $z0 ) // 6 + ($x3 $y1 $z0 ) // 7 + ($x0 $y2 $z0 ) // 8 + ($x1 $y2 $z0 ) // 9 + ($x2 $y2 $z0 ) // 10 + ($x3 $y2 $z0 ) // 11 + + ($x0 $y0 $z1 ) // 12 + ($x1 $y0 $z1 ) // 13 + ($x2 $y0 $z1 ) // 14 + ($x3 $y0 $z1 ) // 15 + ($x0 $y1 $z1 ) // 16 + ($x1 $y1 $z1 ) // 17 + ($x2 $y1 $z1 ) // 18 + ($x3 $y1 $z1 ) // 19 + ($x0 $y2 $z1 ) // 20 + ($x1 $y2 $z1 ) // 21 + ($x2 $y2 $z1 ) // 22 + ($x3 $y2 $z1 ) // 23 +); + +// Grading +h1 30; +h2 3; +v1 15; +v2 30; + +blocks +( + hex ( 0 1 5 4 12 13 17 16 ) + ($h1 $v1 1 ) + simpleGrading (0.5 1 1) + + hex ( 2 3 7 6 14 15 19 18 ) + ($h1 $v1 1) + simpleGrading (2 1 1) + + hex ( 4 5 9 8 16 17 21 20 ) + ($h1 $v2 1) + simpleGrading (0.5 2 1) + + hex ( 5 6 10 9 17 18 22 21 ) + ($h2 $v2 1) + simpleGrading (1 2 1) + + hex ( 6 7 11 10 18 19 23 22 ) + ($h1 $v2 1 ) + simpleGrading (2 2 1) +); + +boundary +( + inlet + { + type patch; + faces + ( + ( 0 4 16 12 ) + ( 4 8 20 16 ) + ); + } + outlet + { + type patch; + faces + ( + ( 3 7 19 15 ) + ( 7 11 23 19 ) + ); + } + flap + { + type wall; + faces + ( + ( 1 5 17 13 ) + ( 5 6 18 17 ) + ( 6 2 14 18 ) + ); + } + upperWall + { + type wall; + faces + ( + ( 8 9 21 20 ) + ( 9 10 22 21 ) + ( 10 11 23 22 ) + ); + } + lowerWall + { + type wall; + faces + ( + ( 0 1 13 12 ) + ( 2 3 15 14 ) + ); + } + frontAndBack + { + type empty; + faces + ( + ( 0 1 5 4 ) + ( 2 3 7 6 ) + ( 4 5 9 8 ) + ( 5 6 10 9 ) + ( 6 7 11 10 ) + ( 12 13 17 16 ) + ( 14 15 19 18 ) + ( 16 17 21 20 ) + ( 17 18 22 21 ) + ( 18 19 23 22 ) + ); + } +); diff --git a/examples/perpendicular-flap/fluid-openfoam/system/controlDict b/examples/perpendicular-flap/fluid-openfoam/system/controlDict new file mode 100644 index 0000000..9d35646 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/system/controlDict @@ -0,0 +1,45 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object controlDict; +} + +application pimpleFoam; // latest OpenFOAM +// application pimpleDyMFoam; // OpenFOAM v1712, OpenFOAM 5.x, or older + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 5; + +deltaT 0.01; + +writeControl adjustableRunTime; + +writeInterval 0.1; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression off; + +timeFormat general; + +timePrecision 6; + +functions +{ + preCICE_Adapter + { + type preciceAdapterFunctionObject; + libs ("libpreciceAdapterFunctionObject.so"); + } +} diff --git a/examples/perpendicular-flap/fluid-openfoam/system/decomposeParDict b/examples/perpendicular-flap/fluid-openfoam/system/decomposeParDict new file mode 100644 index 0000000..31d721d --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/system/decomposeParDict @@ -0,0 +1,17 @@ +FoamFile { + version 2.0; + class dictionary; + object decomposeParDict; + format ascii; +} + +numberOfSubdomains 4; + +method simple; + +simpleCoeffs +{ + n (2 2 1); + delta 0.001; +} + diff --git a/examples/perpendicular-flap/fluid-openfoam/system/fvSchemes b/examples/perpendicular-flap/fluid-openfoam/system/fvSchemes new file mode 100644 index 0000000..80c0961 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/system/fvSchemes @@ -0,0 +1,39 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object fvSchemes; +} + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss linear; +} + +divSchemes +{ + default none; + div(phi,U) bounded Gauss upwind; + div((nuEff*dev2(T(grad(U))))) Gauss linear; +} + +laplacianSchemes +{ + default Gauss linear corrected; +} + +interpolationSchemes +{ + default linear; +} + +snGradSchemes +{ + default corrected; +} diff --git a/examples/perpendicular-flap/fluid-openfoam/system/fvSolution b/examples/perpendicular-flap/fluid-openfoam/system/fvSolution new file mode 100644 index 0000000..064d7f3 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/system/fvSolution @@ -0,0 +1,75 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object fvSolution; +} + +solvers +{ + + p + { + + solver PCG; + preconditioner DIC; + tolerance 1e-8; + relTol 1.0e-3; + } + + pFinal + { + $p; + relTol 0; + } + + pcorr + { + $p; + } + + pcorrFinal + { + $pcorr; + relTol 0; + } + + Phi + { + $p; + } + + "(U|cellDisplacement)" + { + solver smoothSolver; + smoother symGaussSeidel; + tolerance 1e-6; + relTol 1e-4; + minIter 2; + } + + "(U|cellDisplacement)Final" + { + $U; + relTol 0; + } +} + +PIMPLE +{ + nCorrectors 4; + nNonOrthogonalCorrectors 1; + // tolerance 1.0e-14; + // relTol 5e-3; + consistent true; + correctPhi true; + momentumPredictor true; + nOuterCorrectors 1; +} + + +potentialFlow +{ + nNonOrthogonalCorrectors 10; +} diff --git a/examples/perpendicular-flap/fluid-openfoam/system/preciceDict b/examples/perpendicular-flap/fluid-openfoam/system/preciceDict new file mode 100644 index 0000000..1196109 --- /dev/null +++ b/examples/perpendicular-flap/fluid-openfoam/system/preciceDict @@ -0,0 +1,38 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object preciceDict; +} + +preciceConfig "../precice-config.xml"; + +participant Fluid; + +modules (FSI); + +interfaces +{ + Interface1 + { + mesh Fluid-Mesh; + patches (flap); + locations faceCenters; + + readData + ( + Displacement + ); + + writeData + ( + Force + ); + }; +}; + +FSI +{ + rho rho [1 -3 0 0 0 0 0] 1; +} diff --git a/examples/perpendicular-flap/precice-config.xml b/examples/perpendicular-flap/precice-config.xml new file mode 100644 index 0000000..1274f05 --- /dev/null +++ b/examples/perpendicular-flap/precice-config.xml @@ -0,0 +1,64 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/perpendicular-flap/solid-gismo/clean.sh b/examples/perpendicular-flap/solid-gismo/clean.sh new file mode 100755 index 0000000..826662a --- /dev/null +++ b/examples/perpendicular-flap/solid-gismo/clean.sh @@ -0,0 +1,25 @@ +#!/bin/bash + +# Cleaning script for solid-gismo-elasticity directory + +echo "Cleaning unnecessary files..." + +# Remove precice-profiling directory +if [ -d "precice-profiling" ]; then + rm -rf precice-profiling + echo "Deleted 'precice-profiling' folder." +fi + +# Remove precice-run directory +if [ -d "../precice-run" ]; then + rm -rf ../precice-run + echo "Deleted 'precice-run' folder." +fi + +# Remove files ending with .pvd, .vts, .vtp, .log, and .txt +for ext in pvd vts vtp log txt; do + find . -type f -name "*.$ext" -exec rm -f {} \; + echo "Deleted all *.$ext files." +done + +echo "Cleaning completed!" diff --git a/examples/perpendicular-flap/solid-gismo/run.sh b/examples/perpendicular-flap/solid-gismo/run.sh new file mode 100644 index 0000000..cb8fa38 --- /dev/null +++ b/examples/perpendicular-flap/solid-gismo/run.sh @@ -0,0 +1,4 @@ +#!/bin/bash +set -e -u +perpendicular-flap-vertex-gismo -c ../precice-config.xml --plot -r 2 +