Perpendicular flap example with OpenFOAM

Author: Crazy-Rich-Meghan

examples/perpendicular-flap-vertex-gismo.cpp

@@ -0,0 +1,338 @@
+/** @file solid-gismo-elasticity.cpp
+
+    @brief Elasticity participant for the PreCICE example "perpendicular-flap".
+
+
+    This file is part of the 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 (2023 - ..., TU Delft), H.M. Verhelst (2019 - 2024, TU Delft)
+*/
+
+#include <gismo.h>
+#include <gsPreCICE/gsPreCICE.h>
+#include <gsPreCICE/gsPreCICEUtils.h>
+#include <gsPreCICE/gsPreCICEFunction.h>
+// #include <gsPreCICE/gsPreCICEVectorFunction.h>
+
+#include <gsElasticity/gsMassAssembler.h>
+#include <gsElasticity/gsElasticityAssembler.h>
+
+#ifdef gsStructuralAnalysis_ENABLED
+#include <gsStructuralAnalysis/src/gsDynamicSolvers/gsDynamicExplicitEuler.h>
+#include <gsStructuralAnalysis/src/gsDynamicSolvers/gsDynamicImplicitEuler.h>
+#include <gsStructuralAnalysis/src/gsDynamicSolvers/gsDynamicNewmark.h>
+#include <gsStructuralAnalysis/src/gsDynamicSolvers/gsDynamicBathe.h>
+#include <gsStructuralAnalysis/src/gsDynamicSolvers/gsDynamicWilson.h>
+#include <gsStructuralAnalysis/src/gsDynamicSolvers/gsDynamicRK4.h>
+#include <gsStructuralAnalysis/src/gsStructuralAnalysisTools/gsStructuralAnalysisUtils.h>
+#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<patchSide> 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<real_t> 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<index_t> 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<real_t> 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<real_t> massAssembler(patches,bases,bcInfo,gravity);
+    massAssembler.options().setReal("Density",rho);
+    massAssembler.assemble();
+
+    // creating stiffness assembler.
+    gsElasticityAssembler<real_t> assembler(patches,bases,bcInfo,g);
+    assembler.options().setReal("YoungsModulus",E);
+    assembler.options().setReal("PoissonsRatio",nu);
+    assembler.assemble();
+
+    gsMatrix<real_t> 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<gsMatrix<> > 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<real_t>::Zero(assembler.numDofs(),1);
+    V_checkpoint = V = gsVector<real_t>::Zero(assembler.numDofs(),1);
+    A_checkpoint = A = gsVector<real_t>::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<real_t>::Jacobian_t Jacobian = [&assembler,&fixedDofs](gsMatrix<real_t> const &x, gsSparseMatrix<real_t> & m) {
+        // to do: add time dependency of forcing
+        assembler.assemble(x, fixedDofs);
+        m = assembler.matrix();
+        return true;
+    };
+
+    // Function for the Residual
+    gsStructuralAnalysisOps<real_t>::TResidual_t Residual = [&assembler,&fixedDofs](gsMatrix<real_t> const &x, real_t t, gsVector<real_t> & result)
+    {
+        assembler.assemble(x,fixedDofs);
+        result = assembler.rhs();
+        return true;
+    };
+
+
+    gsSparseMatrix<> C = gsSparseMatrix<>(assembler.numDofs(),assembler.numDofs());
+    gsStructuralAnalysisOps<real_t>::Damping_t Damping = [&C](const gsVector<real_t> &, gsSparseMatrix<real_t> & m) { m = C; return true; };
+    gsStructuralAnalysisOps<real_t>::Mass_t    Mass    = [&M](                          gsSparseMatrix<real_t> & m) { m = M; return true; };
+
+    gsDynamicBase<real_t> * timeIntegrator;
+    if (method==1)
+        timeIntegrator = new gsDynamicExplicitEuler<real_t,true>(Mass,Damping,Jacobian,Residual);
+    else if (method==2)
+        timeIntegrator = new gsDynamicImplicitEuler<real_t,true>(Mass,Damping,Jacobian,Residual);
+    else if (method==3)
+        timeIntegrator = new gsDynamicNewmark<real_t,true>(Mass,Damping,Jacobian,Residual);
+    else if (method==4)
+        timeIntegrator = new gsDynamicBathe<real_t,true>(Mass,Damping,Jacobian,Residual);
+    else if (method==5)
+    {
+        timeIntegrator = new gsDynamicWilson<real_t,true>(Mass,Damping,Jacobian,Residual);
+        timeIntegrator->options().setReal("gamma",1.4);
+    }
+    else if (method==6)
+        timeIntegrator = new gsDynamicRK4<real_t,true>(Mass,Damping,Jacobian,Residual);
+    else
+        GISMO_ERROR("Method "<<method<<" not known");
+
+    timeIntegrator->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<real_t> 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|| = "<<U.norm()<<"\n";
+            gsInfo<<"\t ||V|| = "<<V.norm()<<"\n";
+            gsInfo<<"\t ||A|| = "<<A.norm()<<"\n";
+
+
+            timestep_checkpoint = timestep;
+        }
+
+        assembler.assemble();
+        F = assembler.rhs();
+
+        // solve gismo timestep
+        gsInfo << "Solving timestep " << time << "...\n";
+        timeIntegrator->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<<time;
+            writer.add(pointDataMatrix,otherDataMatrix);
+        }
+    }
+
+    if (plot)
+    {
+        collection.save();
+    }
+
+
+    return  EXIT_SUCCESS;
+}

examples/perpendicular-flap/README.md

@@ -0,0 +1,59 @@
+---
+title: Perpendicular Flap (with IGA Solid Participant Communicating Stress Data)
+keywords: G+Smo, perpendicular flap
+summary: This tutorial demonstrates the isogeometric analysis (IGA) solid solver version of the “Perpendicular Flap” tutorial. It focuses on using G+Smo to handle solid-structure interactions by exchanging stress data during simulations.
+---
+
+
+## Overview
+
+This example demostrates how the Geometry + Simulation Modules (G+Smo) can be utilised through a module-type adapter to couple with other codes using preCICE. G+Smo offers a robust framework for isogeometric analysis (IGA), seamlessly integrating geometric representations with numerical solvers to enable advanced simulations and efficient code coupling.
+
+## Setup
+This tutorial uses the same setup as the **original** perpendicular flap example, with the key difference being the switch in the communicated information to stress. 
+
+## Requirements
+
+To run the tutorial you need to install the following components:
+- [preCICE](https://precice.org/quickstart.html)
+- [G+Smo and gsPreCICE](https://github.com/gismo/gismo)
+- [OpenFOAM](https://openfoam.org/download/)
+
+## Run the tutorial
+
+In the G+Smo build folder, build the tutorial file.
+
+```
+make <participant file> -j <#threads>
+```
+
+Go to the gismo-executable folder and link the compiled executable to the gismo_executable.
+
+```bash
+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:
+
+for the fluid participant with OpenFOAM, 
+```bash
+cd fluid-openfoam
+OpenFOAM<version> // 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)