From 1e37ac42bee1288c6686e3ff5ca892e27e19ca1d Mon Sep 17 00:00:00 2001 From: Crazy-Rich-Meghan Date: Thu, 5 Dec 2024 17:12:12 +0100 Subject: [PATCH 01/14] Reorder directories --- ...titioned-heat-conduction-IGA-Dirichlet.cpp | 0 ...artitioned-heat-conduction-IGA-Neumann.cpp | 0 .../gismo-left/clean.sh | 6 - .../gismo-left/run.sh | 4 - .../gismo-right/clean.sh | 6 - .../gismo-right/run.sh | 4 - .../partitioned-heat-conduction-direct.cpp | 407 ------------- .../precice-config.xml | 77 --- .../gismo-left/clean.sh | 6 - .../gismo-left/run.sh | 4 - .../gismo-right/clean.sh | 6 - .../gismo-right/run.sh | 4 - .../partitioned-heat-conduction.cpp | 36 +- .../partitioned-heat-conduction/.DS_Store | Bin 0 -> 8196 bytes .../partitioned-heat-conduction/README.md | 68 +++ .../.DS_Store | Bin 0 -> 8196 bytes .../dirichlet-gismo}/.DS_Store | Bin 6148 -> 6148 bytes .../dirichlet-gismo/clean.sh | 27 + .../dirichlet-gismo/run.sh | 4 + .../gismo-executable}/.DS_Store | Bin 6148 -> 6148 bytes .../gismo-executable/creat_symlink.sh | 48 ++ .../neumann-gismo}/.DS_Store | Bin 6148 -> 6148 bytes .../neumann-gismo/clean.sh | 27 + .../neumann-gismo/run.sh | 4 + .../precice-config.xml | 0 .../.DS_Store | Bin 6148 -> 6148 bytes .../dirichlet-gismo/clean.sh | 27 + .../dirichlet-gismo/run.sh | 4 + .../gismo-executable/.DS_Store | Bin 0 -> 6148 bytes .../gismo-executable/creat_symlink.sh | 26 + .../neumann-gismo/clean.sh | 27 + .../neumann-gismo/run.sh | 4 + .../precice-config.xml | 44 +- .../perpendicular-flap-IGA-IGA-fluid.cpp | 0 .../perpendicular-flap-IGA-IGA-solid.cpp | 0 .../precice-config.xml | 84 --- .../precice-config.xml | 72 --- .../solid-gismo-IGA.cpp | 463 -------------- .../vertical-beam-vertex-vertex-fluid.cpp | 211 ------- .../precice-config.xml | 74 --- .../perpendicular-flap-solid-gismo.cpp | 564 ------------------ .../vertical-beam-vertex-vertex-solid.cpp | 430 ------------- src/gsLookupFunction.h => gsLookupFunction.h | 40 +- src/gsPreCICE.h => gsPreCICE.h | 42 +- ...gsPreCICEFunction.h => gsPreCICEFunction.h | 4 +- src/gsPreCICEUtils.h => gsPreCICEUtils.h | 44 +- 46 files changed, 365 insertions(+), 2533 deletions(-) rename examples/{partitioned-heat-conduction-IGA-IGA/dirichlet => }/partitioned-heat-conduction-IGA-Dirichlet.cpp (100%) rename examples/{partitioned-heat-conduction-IGA-IGA/neumann => }/partitioned-heat-conduction-IGA-Neumann.cpp (100%) delete mode 100755 examples/partitioned-heat-conduction-direct/gismo-left/clean.sh delete mode 100755 examples/partitioned-heat-conduction-direct/gismo-left/run.sh delete mode 100755 examples/partitioned-heat-conduction-direct/gismo-right/clean.sh delete mode 100755 examples/partitioned-heat-conduction-direct/gismo-right/run.sh delete mode 100644 examples/partitioned-heat-conduction-direct/partitioned-heat-conduction-direct.cpp delete mode 100644 examples/partitioned-heat-conduction-direct/precice-config.xml delete mode 100755 examples/partitioned-heat-conduction-vertex-vertex/gismo-left/clean.sh delete mode 100755 examples/partitioned-heat-conduction-vertex-vertex/gismo-left/run.sh delete mode 100755 examples/partitioned-heat-conduction-vertex-vertex/gismo-right/clean.sh delete mode 100755 examples/partitioned-heat-conduction-vertex-vertex/gismo-right/run.sh rename examples/{partitioned-heat-conduction-vertex-vertex => }/partitioned-heat-conduction.cpp (92%) create mode 100644 examples/partitioned-heat-conduction/.DS_Store create mode 100644 examples/partitioned-heat-conduction/README.md create mode 100644 examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/.DS_Store rename examples/{partitioned-heat-conduction-direct => partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo}/.DS_Store (84%) create mode 100755 examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/clean.sh create mode 100755 examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/run.sh rename examples/{partitioned-heat-conduction-vertex-vertex => partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/gismo-executable}/.DS_Store (84%) create mode 100755 examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/gismo-executable/creat_symlink.sh rename examples/{partitioned-heat-conduction-IGA-IGA => partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo}/.DS_Store (77%) create mode 100755 examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/clean.sh create mode 100755 examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/run.sh rename examples/{ => partitioned-heat-conduction}/partitioned-heat-conduction-IGA-IGA/precice-config.xml (100%) rename examples/{perpendicular-flap-vertex-vertex => partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex}/.DS_Store (76%) create mode 100755 examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/dirichlet-gismo/clean.sh create mode 100755 examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/dirichlet-gismo/run.sh create mode 100644 examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/gismo-executable/.DS_Store create mode 100755 examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/gismo-executable/creat_symlink.sh create mode 100755 examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/neumann-gismo/clean.sh create mode 100755 examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/neumann-gismo/run.sh rename examples/{ => partitioned-heat-conduction}/partitioned-heat-conduction-vertex-vertex/precice-config.xml (71%) rename examples/{perpendicular-flap-IGA-IGA => }/perpendicular-flap-IGA-IGA-fluid.cpp (100%) rename examples/{perpendicular-flap-IGA-IGA => }/perpendicular-flap-IGA-IGA-solid.cpp (100%) delete mode 100644 examples/perpendicular-flap-IGA-IGA/precice-config.xml delete mode 100644 examples/perpendicular-flap-IGA-Spine/precice-config.xml delete mode 100644 examples/perpendicular-flap-IGA-Spine/solid-gismo-IGA.cpp delete mode 100644 examples/perpendicular-flap-vertex-vertex/dummy-fluid-gismo/vertical-beam-vertex-vertex-fluid.cpp delete mode 100644 examples/perpendicular-flap-vertex-vertex/precice-config.xml delete mode 100644 examples/perpendicular-flap-vertex-vertex/solid-gismo/perpendicular-flap-solid-gismo.cpp delete mode 100644 examples/perpendicular-flap-vertex-vertex/solid-gismo/vertical-beam-vertex-vertex-solid.cpp rename src/gsLookupFunction.h => gsLookupFunction.h (80%) rename src/gsPreCICE.h => gsPreCICE.h (87%) rename src/gsPreCICEFunction.h => gsPreCICEFunction.h (97%) rename src/gsPreCICEUtils.h => gsPreCICEUtils.h (89%) diff --git a/examples/partitioned-heat-conduction-IGA-IGA/dirichlet/partitioned-heat-conduction-IGA-Dirichlet.cpp b/examples/partitioned-heat-conduction-IGA-Dirichlet.cpp similarity index 100% rename from examples/partitioned-heat-conduction-IGA-IGA/dirichlet/partitioned-heat-conduction-IGA-Dirichlet.cpp rename to examples/partitioned-heat-conduction-IGA-Dirichlet.cpp diff --git a/examples/partitioned-heat-conduction-IGA-IGA/neumann/partitioned-heat-conduction-IGA-Neumann.cpp b/examples/partitioned-heat-conduction-IGA-Neumann.cpp similarity index 100% rename from examples/partitioned-heat-conduction-IGA-IGA/neumann/partitioned-heat-conduction-IGA-Neumann.cpp rename to examples/partitioned-heat-conduction-IGA-Neumann.cpp diff --git a/examples/partitioned-heat-conduction-direct/gismo-left/clean.sh b/examples/partitioned-heat-conduction-direct/gismo-left/clean.sh deleted file mode 100755 index 708528e..0000000 --- a/examples/partitioned-heat-conduction-direct/gismo-left/clean.sh +++ /dev/null @@ -1,6 +0,0 @@ -#!/bin/sh -set -e -u - -. ../../tools/cleaning-tools.sh - -clean_gismo . diff --git a/examples/partitioned-heat-conduction-direct/gismo-left/run.sh b/examples/partitioned-heat-conduction-direct/gismo-left/run.sh deleted file mode 100755 index 7630e55..0000000 --- a/examples/partitioned-heat-conduction-direct/gismo-left/run.sh +++ /dev/null @@ -1,4 +0,0 @@ -#!/bin/bash -set -e -u - -../../../../../build/bin/partitioned-heat-conduction -c ../precice-config.xml --plot -s 0 diff --git a/examples/partitioned-heat-conduction-direct/gismo-right/clean.sh b/examples/partitioned-heat-conduction-direct/gismo-right/clean.sh deleted file mode 100755 index 708528e..0000000 --- a/examples/partitioned-heat-conduction-direct/gismo-right/clean.sh +++ /dev/null @@ -1,6 +0,0 @@ -#!/bin/sh -set -e -u - -. ../../tools/cleaning-tools.sh - -clean_gismo . diff --git a/examples/partitioned-heat-conduction-direct/gismo-right/run.sh b/examples/partitioned-heat-conduction-direct/gismo-right/run.sh deleted file mode 100755 index 026edd8..0000000 --- a/examples/partitioned-heat-conduction-direct/gismo-right/run.sh +++ /dev/null @@ -1,4 +0,0 @@ -#!/bin/bash -set -e -u - -../../../../../build/bin/partitioned-heat-conduction -c ../precice-config.xml --plot -s 1 diff --git a/examples/partitioned-heat-conduction-direct/partitioned-heat-conduction-direct.cpp b/examples/partitioned-heat-conduction-direct/partitioned-heat-conduction-direct.cpp deleted file mode 100644 index 776ae84..0000000 --- a/examples/partitioned-heat-conduction-direct/partitioned-heat-conduction-direct.cpp +++ /dev/null @@ -1,407 +0,0 @@ -/** @file heat-equation-coupling.cpp - - @brief Heat equation participant for a double coupled heat equation - - 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): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) -*/ - -#include -#include -#include - -using namespace gismo; - -int main(int argc, char *argv[]) -{ - //! [Parse command line] - bool plot = false; - index_t plotmod = 1; - index_t numRefine = 5; - index_t numElevate = 0; - short_t side = 0; - 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","plotmod", "Modulo for plotting, i.e. if plotmod==1, plots every timestep", plotmod); - cmd.addInt("s","side", "Patchside of interface", side); - try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } - - //! [Read input file] - - gsMultiPatch<> patches; - if (side==0) //left - patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,0.0,1.0,1.0)); - else if (side==1) //right - patches.addPatch(gsNurbsCreator<>::BSplineRectangle(1.0,0.0,2.0,1.0)); - else - GISMO_ERROR("Side unknown"); - - real_t alpha = 3; - real_t beta = 1.3; - real_t time = 0; - real_t k_temp = 1; - - // Set external heat-flux to zero - gsConstantFunction<> f(beta-2-2*alpha,2); - gsFunctionExpr<> u_ex("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); - - gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) - - 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"; - -// ---------------------------------------------------------------------------------------------- - typedef gsExprAssembler<>::geometryMap geometryMap; - typedef gsExprAssembler<>::space space; - typedef gsExprAssembler<>::solution solution; - - gsExprAssembler<> A(1,1); - - gsInfo<<"Active options:\n"<< A.options() <<"\n"; - - A.setIntegrationElements(bases); - - gsExprEvaluator<> ev(A); - -// ---------------------------------------------------------------------------------------------- - boxSide couplingSide; - if (side==0) //left - couplingSide = boxSide(2); - else if (side==1) //right - couplingSide = boxSide(1); - else - GISMO_ERROR("Side unknown"); - - patchSide couplingInterface(0,couplingSide); - typename gsBasis::domainIter domIt = bases.basis(couplingInterface.patch).makeDomainIterator(couplingInterface.side()); - index_t rows = patches.targetDim(); - gsMatrix<> nodes; - // Start iteration over elements - gsVector<> tmp; - index_t k=0; - - gsOptionList quadOptions = A.options(); - // NEED A DIFFERENT RULE FOR dirichlet::interpolation --> see: gsGaussRule bdQuRule(basis, 1.0, 1, iter->side().direction()); - /* - quadOptions.addInt("quRule","Quadrature rule [1:GaussLegendre,2:GaussLobatto]",1); - quadOptions.addReal("quA","Number of quadrature points: quA*deg + quB",1.0); - quadOptions.addInt("quB","Number of quadrature points: quA*deg + quB",1); - */ - - index_t quadSize = 0; - typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT - for (; domIt->good(); domIt->next(), k++ ) - { - QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); - quadSize+=QuRule->numNodes(); - } - gsMatrix<> uv(rows,quadSize); - gsMatrix<> xy(rows,quadSize); - - index_t offset = 0; - - for (domIt->reset(); domIt->good(); domIt->next(), k++ ) - { - QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); - // Map the Quadrature rule to the element - QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(), - nodes, tmp); - uv.block(0,offset,rows,QuRule->numNodes()) = nodes; - - gsMatrix<> tmp2; - patches.patch(couplingInterface.patch).eval_into(nodes,tmp2); - xy.block(0,offset,rows,QuRule->numNodes()) = patches.patch(couplingInterface.patch).eval(nodes); - offset += QuRule->numNodes(); - } - - GISMO_ASSERT(side==0 || side==1,"Side must be west or east"); - std::string readName = (side == 0) ? "Dirichlet" : "Neumann"; - std::string writeName = (side == 0) ? "Neumann" : "Dirichlet"; - - gsPreCICE interface(membername, precice_config); - std::string readMeshName = readName + "-Mesh"; - std::string writeMeshName = writeName + "-Mesh"; - interface.addMesh(readMeshName,xy); - - interface.requiresInitialData(); - - real_t precice_dt = interface.initialize(); - - std::string tempName = "Temperature"; - std::string fluxName = "Heat-Flux"; - - gsMatrix<> bbox(2,2); - bbox<<0.9,1.1,-0.1,1.1; - interface.setMeshAccessRegion(writeMeshName,bbox); - - std::vector writeIDs; - gsMatrix<> writePoints; - interface.getMeshVertexIDsAndCoordinates(writeMeshName,writeIDs,writePoints); - -/* - TO DO: - * Write the evaluation on the other coordinates on the other IDs -*/ -// ---------------------------------------------------------------------------------------------- - - gsBoundaryConditions<> bcInfo; - gsPreCICEFunction g_CD(&interface,readMeshName,(side==0 ? tempName : fluxName),patches,1); - gsPreCICEFunction g_CN(&interface,readMeshName,(side==0 ? tempName : fluxName),patches,1); - gsFunction<> * g_C = &u_ex; - - bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &u_ex, 0, false, 0); - bcInfo.addCondition(0, boundary::north, condition_type::dirichlet, &u_ex, 0, false, 0); - if (side==0) - { - bcInfo.addCondition(0, couplingSide,condition_type::dirichlet , g_C, 0, false, 0); - bcInfo.addCondition(0, couplingSide.opposite(),condition_type::dirichlet, &u_ex, 0, false, 0); - } - else - { - bcInfo.addCondition(0, couplingSide,condition_type::neumann , &g_CN); - bcInfo.addCondition(0, couplingSide.opposite(),condition_type::dirichlet, &u_ex, 0, false, 0); - } - - bcInfo.setGeoMap(patches); - // gsDebugVar(bcInfo); - -// ---------------------------------------------------------------------------------------------- - - // Generate system matrix and load vector - // gsInfo << "Assembling mass and stiffness...\n"; - - // Set the geometry map - geometryMap G = A.getMap(patches); - - // Set the discretization space - space u = A.getSpace(bases); - - // Set the source term - auto ff = A.getCoeff(f, G); - - // Set the solution - gsMatrix<> solVector, solVector_ini; - solution u_sol = A.getSolution(u, solVector); - solution u_ini = A.getSolution(u, solVector); - - u.setup(bcInfo, dirichlet::homogeneous, 0); - A.initSystem(); - gsDebugVar(A.numDofs()); - A.assemble( u * u.tr() * meas(G)); - gsSparseMatrix<> M = A.matrix(); - - - // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner - gsSparseSolver<>::CGDiagonal solver; - - real_t dt = 0.01; - - // Project u_wall as initial condition (violates Dirichlet side on precice interface) - auto uex = A.getCoeff(u_ex, G); - // RHS of the projection - u.setup(bcInfo, dirichlet::l2Projection, 0); - A.initSystem(); - A.assemble( u * u.tr() * meas(G), u * uex * meas(G) ); - solver.compute(A.matrix()); - solVector_ini = solVector = solver.solve(A.rhs()); - - gsMatrix<> result(1,uv.cols()), tmp2; - // Write initial data - if (interface.requiresWritingCheckpoint()) - { - for (index_t k=0; k!=uv.cols(); k++) - { - if (side==0) - { - gsWarn<<"Write the flux here!!!\n"; - tmp2 = ev.eval( - jac(u_sol) * nv(G).normalized(),writePoints.col(k)); - result(0,k) = tmp2.at(0); - } - else - { - tmp2 = ev.eval(u_sol,writePoints.col(k)); - result(0,k) = tmp2.at(0); - } - } - // gsDebugVar(result); - interface.writeData(writeMeshName,writeName=="Dirichlet" ? fluxName : tempName,writeIDs,result); - } - // interface.initialize_data(); - - // interface.readBlockScalarData(readMeshName,side==0 ? tempName : fluxName,xy,result); - // gsDebugVar(result); - - // Initialize the RHS for assembly - if (side==0) - g_C = &g_CD; - gsDebugVar(bcInfo); - A.initSystem(); - A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * uex * meas(G) ); - gsSparseMatrix<> K = A.matrix(); - - // Assemble the RHS - gsVector<> F = dt*A.rhs() + M*solVector; - gsVector<> F0 = F; - gsVector<> F_checkpoint = F; - gsMatrix<> solVector_checkpoint = solVector; - - gsParaviewCollection collection("solution"); - gsParaviewCollection exact_collection("exact_solution"); - - index_t timestep = 0; - index_t timestep_checkpoint = 0; - real_t t_checkpoint = 0; - if (plot) - { - std::string fileName = "solution_" + util::to_string(timestep); - ev.options().setSwitch("plot.elements", true); - ev.options().setInt("plot.npts", 1000); - ev.writeParaview( u_sol , G, fileName); - for (size_t p=0; p!=patches.nPatches(); p++) - { - fileName = "solution_" + util::to_string(timestep) + std::to_string(p); - collection.addTimestep(fileName,time,".vts"); - } - - // fileName = "exact_solution_" + util::to_string(timestep); - // ev.writeParaview( uex , G, fileName); - // for (size_t p=0; p!=patches.nPatches(); p++) - // { - // fileName = "exact_solution_" + util::to_string(timestep) + std::to_string(p); - // exact_collection.addTimestep(fileName,time,".vts"); - // } - - ev.writeParaview( u_sol , G, "initial_solution"); - - } - - // time += precice_dt; - while (interface.isCouplingOngoing()) - { - // read temperature from interface - // if (interface.isReadDataAvailable()) - // { - u_ex = gsFunctionExpr<>("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); - u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: - - A.initSystem(); - A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * ff * meas(G) ); - K = A.matrix(); - auto g_Neumann = A.getBdrFunction(G); - A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); - F = dt*A.rhs() + M*solVector; - - // gsMatrix<> result; - // interface.readBlockScalarData(readMeshName,side==0 ? tempName : fluxName,xy,result); - // gsDebugVar(result); - // } - - // save checkpoint - if (interface.requiresWritingCheckpoint()) - { - gsDebugVar("Write checkpoint"); - F_checkpoint = F0; - t_checkpoint = time; - timestep_checkpoint = timestep; - solVector_checkpoint = solVector; - } - - // potentially adjust non-matching timestep sizes - dt = std::min(dt,precice_dt); - - // solve gismo timestep - gsInfo << "Solving timestep " << timestep*dt << "..."; - solVector = solver.compute(M + dt*K).solve(F); - gsInfo<<"Finished\n"; - // write heat fluxes to interface - gsMatrix<> result(1,uv.cols()), tmp; - for (index_t k=0; k!=uv.cols(); k++) - { - // gsDebugVar(ev.eval(nv(G),uv.col(k))); - // tmp = ev.eval(k_temp * igrad(u_sol,G),uv.col(k)); - // Only exchange y component - // result(0,k) = -tmp.at(1); - // - // - if (side==0) - { - gsWarn<<"Write the flux here!!!\n"; - tmp = ev.eval(jac(u_sol) * nv(G).normalized(),uv.col(k)); - result(0,k) = tmp.at(0); - } - else - { - tmp = ev.eval(u_sol,uv.col(k)); - result(0,k) = tmp.at(0); - } - } - // TODO - interface.writeData(readMeshName,side==0 ? fluxName : tempName,xy,result); - - // do the coupling - precice_dt = interface.advance(dt); - - // advance variables - time += dt; - timestep += 1; - F0 = F; - - if (interface.requiresReadingCheckpoint()) - { - gsDebugVar("Read checkpoint"); - F0 = F_checkpoint; - time = t_checkpoint; - timestep = timestep_checkpoint; - solVector = solVector_checkpoint; - } - else - { - if (timestep % plotmod==0 && plot) - { - std::string fileName = "solution_" + util::to_string(timestep); - ev.options().setSwitch("plot.elements", true); - ev.options().setInt("plot.npts", 1000); - ev.writeParaview( u_sol , G, fileName); - for (size_t p=0; p!=patches.nPatches(); p++) - { - fileName = "solution_" + util::to_string(timestep) + std::to_string(p); - collection.addTimestep(fileName,time,".vts"); - } - - // fileName = "exact_solution_" + util::to_string(timestep); - // ev.writeParaview( uex , G, fileName); - // for (size_t p=0; p!=patches.nPatches(); p++) - // { - // fileName = "exact_solution_" + util::to_string(timestep) + std::to_string(p); - // exact_collection.addTimestep(fileName,time,".vts"); - // } - } - } - } - - if (plot) - { - collection.save(); - exact_collection.save(); - } - - - return EXIT_SUCCESS; -} diff --git a/examples/partitioned-heat-conduction-direct/precice-config.xml b/examples/partitioned-heat-conduction-direct/precice-config.xml deleted file mode 100644 index 46d6bdf..0000000 --- a/examples/partitioned-heat-conduction-direct/precice-config.xml +++ /dev/null @@ -1,77 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/clean.sh b/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/clean.sh deleted file mode 100755 index 708528e..0000000 --- a/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/clean.sh +++ /dev/null @@ -1,6 +0,0 @@ -#!/bin/sh -set -e -u - -. ../../tools/cleaning-tools.sh - -clean_gismo . diff --git a/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/run.sh b/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/run.sh deleted file mode 100755 index 7630e55..0000000 --- a/examples/partitioned-heat-conduction-vertex-vertex/gismo-left/run.sh +++ /dev/null @@ -1,4 +0,0 @@ -#!/bin/bash -set -e -u - -../../../../../build/bin/partitioned-heat-conduction -c ../precice-config.xml --plot -s 0 diff --git a/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/clean.sh b/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/clean.sh deleted file mode 100755 index 708528e..0000000 --- a/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/clean.sh +++ /dev/null @@ -1,6 +0,0 @@ -#!/bin/sh -set -e -u - -. ../../tools/cleaning-tools.sh - -clean_gismo . diff --git a/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/run.sh b/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/run.sh deleted file mode 100755 index 026edd8..0000000 --- a/examples/partitioned-heat-conduction-vertex-vertex/gismo-right/run.sh +++ /dev/null @@ -1,4 +0,0 @@ -#!/bin/bash -set -e -u - -../../../../../build/bin/partitioned-heat-conduction -c ../precice-config.xml --plot -s 1 diff --git a/examples/partitioned-heat-conduction-vertex-vertex/partitioned-heat-conduction.cpp b/examples/partitioned-heat-conduction.cpp similarity index 92% rename from examples/partitioned-heat-conduction-vertex-vertex/partitioned-heat-conduction.cpp rename to examples/partitioned-heat-conduction.cpp index 80c2551..54a222c 100644 --- a/examples/partitioned-heat-conduction-vertex-vertex/partitioned-heat-conduction.cpp +++ b/examples/partitioned-heat-conduction.cpp @@ -8,7 +8,7 @@ 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): H.M. Verhelst + Author(s): H.M. Verhelst (2019-2024, TUDelft), J. Li (2023-..., TUDelft) */ #include @@ -22,7 +22,7 @@ int main(int argc, char *argv[]) //! [Parse command line] bool plot = false; index_t plotmod = 1; - index_t numRefine = 5; + index_t numRefine = 2; // Number of uniform refinement (increase the number also requires to modify the rbf radius in the precice_config.xml) index_t numElevate = 0; short_t side = 0; std::string precice_config("../precice_config.xml"); @@ -34,7 +34,7 @@ int main(int argc, char *argv[]) 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("s","side", "Patchside of interface", side); + cmd.addInt("s","side", "Patchside of interface: 0 (Dirichlet), 1 (Neumann)", side); try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } //! [Read input file] @@ -112,8 +112,8 @@ int main(int argc, char *argv[]) QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); quadSize+=QuRule->numNodes(); } - gsMatrix<> uv(rows,quadSize); - gsMatrix<> xy(rows,quadSize); + gsMatrix<> uv(rows,quadSize); // Coordinates of the quadrature points in parameter space + gsMatrix<> xy(rows,quadSize); // Coordinates of the quadrature points in physical space index_t offset = 0; @@ -141,8 +141,6 @@ int main(int argc, char *argv[]) gsPreCICE interface(membername, precice_config); std::string meshName = membername + "-Mesh"; - - gsDebugVar(xy.dim()); interface.addMesh(meshName,xy); interface.requiresInitialData(); @@ -201,7 +199,7 @@ int main(int argc, char *argv[]) gsSparseMatrix<> M = A.matrix(); - // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner + // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner gsSparseSolver<>::CGDiagonal solver; real_t dt = 0.01; @@ -215,7 +213,7 @@ int main(int argc, char *argv[]) solver.compute(A.matrix()); solVector_ini = solVector = solver.solve(A.rhs()); - gsMatrix<> result(1,uv.cols()), tmp2; + gsMatrix<> result(1,uv.cols()), tmp2; // Write initial data if (interface.requiresWritingCheckpoint()) { @@ -234,8 +232,6 @@ int main(int argc, char *argv[]) } } gsDebugVar(result); - gsDebugVar(xy); - // gsDebugVar(result); interface.writeData(meshName,side==0 ? fluxName : tempName,xy,result); } // interface.initialize_data(); @@ -291,10 +287,8 @@ int main(int argc, char *argv[]) while (interface.isCouplingOngoing()) { // read temperature from interface - // if (interface.isReadDataAvailable()) - // { u_ex = gsFunctionExpr<>("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); - u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: + u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: dirichlet::l2Projection is used to project the Dirichlet BCs A.initSystem(); A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * ff * meas(G) ); @@ -303,11 +297,6 @@ int main(int argc, char *argv[]) A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); F = dt*A.rhs() + M*solVector; - // gsMatrix<> result; - // interface.readBlockScalarData(meshName,side==0 ? tempName : fluxName,xy,result); - // gsDebugVar(result); - // } - // save checkpoint if (interface.requiresWritingCheckpoint()) { @@ -347,7 +336,6 @@ int main(int argc, char *argv[]) result(0,k) = tmp.at(0); } } - // TODO interface.writeData(meshName,side==0 ? fluxName : tempName,xy,result); // do the coupling @@ -379,14 +367,6 @@ int main(int argc, char *argv[]) fileName = "solution_" + util::to_string(timestep) + std::to_string(p); collection.addTimestep(fileName,time,".vts"); } - - // fileName = "exact_solution_" + util::to_string(timestep); - // ev.writeParaview( uex , G, fileName); - // for (size_t p=0; p!=patches.nPatches(); p++) - // { - // fileName = "exact_solution_" + util::to_string(timestep) + std::to_string(p); - // exact_collection.addTimestep(fileName,time,".vts"); - // } } } } diff --git a/examples/partitioned-heat-conduction/.DS_Store b/examples/partitioned-heat-conduction/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..aeb464e46745c2c1ec31f573000d8177c362c00c GIT binary patch literal 8196 zcmeHMU2hUW6ur};UDBBFphmGm zUE^4xygE?G5CE7%w@j!b4v-vMV_oA|p`>C?6+H+|6*|NaijHxc&4G1|V}*)NLeWX+ zm4%K_gj^jwTbh%oD>S81z$h@Q0LSk4s822>GIswHzx!08L&Whe`b%^mqk5kPs7-82 zpd?R#PDVaA_cKMc}rYwMR-S()o-FdjOYn+?*2o4n`N8*IUw)Lm~s=nUJ1+s!!2+#u>Ck`VYEn0$B_M1DMM#DmCB zCD%~{R^H0D3wKXWwyTeecC}U;744HsrS_y~Z&x3UMtSSbX8Gw}`%Bo1;%`b)p|2&n zrP|Z-J8G9C$3$VC`;w+Sa*`;Fqb_FT4{(2Bza8Na@B!-3F}!!GbzOY81ZhtDa|zR* z^KG8?Xa~<$i*|tzs7HQ3^Ii_pEP>~LkV9W7l^*O-nfCFxhj`d~@Mwy5OI`ymYH!hI zUvuwg*aY;DDyZ**-vrCoj4nn}@MCJ@RH;iMz(vtPn|t{x3DGK}-k)A6F8a#U1@pD! zG;tF&vEMQ*kFkdyY?;NYFLr4y$~$wa&z^DYg=oIMue82%M!_-3C~$=etg5v)IR4+c z`1k)S+@V?8C}0#=wgMv8a#{_nGWimFenF14ZRAhLoY-!xP*PCHa2zPZap3e1L)2}U dGN-O_tPnj&fB%Po#lAO~_rG+PDXxhEzX6h%RJ;HH literal 0 HcmV?d00001 diff --git a/examples/partitioned-heat-conduction/README.md b/examples/partitioned-heat-conduction/README.md new file mode 100644 index 0000000..1009251 --- /dev/null +++ b/examples/partitioned-heat-conduction/README.md @@ -0,0 +1,68 @@ + +## Setup + +We solve a partitioned heat equation. For information on the non-partitioned case, please refer to [1, p.37ff]. In this tutorial the computational domain is partitioned and coupled via preCICE. The coupling roughly follows the approach described in [2]. + +![Case setup of partitioned-heat-conduction case](https://github.com/precice/tutorials/blob/master/partitioned-heat-conduction/images/tutorials-partitioned-heat-conduction-setup.png) + +Case setup from [2]. `D` denotes the Dirichlet participant and `N` denotes the Neumann participant. + +The heat equation is solved on a rectangular domain `Omega = [0,2] x [0,1]` with given Dirichlet boundary conditions. We split the domain at `x_c = 1` using a straight vertical line, the coupling interface. The left part of the domain will be referred to as the Dirichlet partition and the right part as the Neumann partition. To couple the two participants we use Dirichlet-Neumann coupling. Here, the Dirichlet participant receives Dirichlet boundary conditions (`Temperature`) at the coupling interface and solves the heat equation using these boundary conditions on the left part of the domain. Then the Dirichlet participant computes the resulting heat flux (`Flux`) from the solution and sends it to the Neumann participant. The Neumann participant uses the flux as a Neumann boundary condition to solve the heat equation on the right part of the domain. We then extract the temperature from the solution and send it back to the Dirichlet participant. This establishes the coupling between the two participants. + +This simple case allows us to compare the solution for the partitioned case to a known analytical solution (method of manufactures solutions, see [1, p.37ff]). For more usage examples and details, please refer to [3, sect. 4.1]. + +## Configuration + +preCICE configuration (image generated using the [precice-config-visualizer](https://precice.org/tooling-config-visualization.html)): + +![preCICE configuration visualization](https://github.com/precice/tutorials/blob/master/partitioned-heat-conduction/images/tutorials-partitioned-heat-conduction-precice-config.png) + + +## Running the simulation + +In the G+Smo build folder, build the tutorial file. + +``` +make partitioned-heat-conduction -j <#threads> +make partitioned-heat-conduction-IGA-dirichlet -j <#threads> +make partitioned-heat-conduction-IGA-neumann -j <#threads> +``` + +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: + +```bash +cd dirichlet-gismo +./run.sh +``` + +and + +```bash +cd neumann-gismo +./run.sh +``` + + + +## Visualization + + +For G+Smo, please use the file `solution.pvd` in both dirichlet-gismo and neumann-gismo directories. + +![Animation of the partitioned heat equation](https://github.com/Crazy-Rich-Meghan/tutorials/blob/partitioned-heat-conduction/partitioned-heat-conduction/images/tutorial-partitioned-heat-conduction-gismo.gif) + + + +## References + +[1] Azahar Monge and Philipp Birken. "Convergence Analysis of the Dirichlet-Neumann Iteration for Finite Element Discretizations." (2016). Proceedings in Applied Mathematics and Mechanics. [doi](https://doi.org/10.1002/pamm.201610355) +[2] Benjamin Rüth, Benjamin Uekermann, Miriam Mehl, Philipp Birken, Azahar Monge, and Hans Joachim Bungartz. "Quasi-Newton waveform iteration for partitioned surface-coupled multiphysics applications." (2020). International Journal for Numerical Methods in Engineering. [doi](https://doi.org/10.1002/nme.6443) +[3] Tobias Eppacher. "Parallel-in-Time Integration with preCICE" (2024). Bachlor's thesis at Technical University of Munich. [pdf](https://mediatum.ub.tum.de/1755012) diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/.DS_Store b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..fcf242967c60dc3550512971a6603daeac0efcef GIT binary patch literal 8196 zcmeHMO^?z*7=C9#w`8*slQmgAHSxMeST?d76VpY!day5cO{uoa_^Rd_f+v;juqBF_NGwr<3d_0{G-hqf%V;oe8_KCOi0-Udn-@j_qkvJsC}0#Y3T%P`II~5t&N=s0sY#6jMuGoQ0kJ>W zhyoprbA@v2Kq6fLU=gonKpuX8C4p{#oAhpM{-CtGxG+$0IL{J3~nb zgAPnyJ`dx7Pa1q22Z`jms=z8)rB?aobXu)cDhGD$@OXA$Pfre!KdRlE%}Umd+js9j zX}yX1F@LW{#`Nw$yGtYekUz{IX-7V~EtolqMBS>}=KVpcksL;!H=9~RtFU+L^Mlr> z)BJY&Jh%XB(EHWiy>*9i#N#e{Bm^we;@jgMr$g@$+cd_ljz_+@-?{7=aa23Tgs9iW zgs4{QCl{R&xdTnliCp6)6?Z9yZey$<54jjhE>yj2{r=R7${1o=xeT*4D=EV|y^3a! zzkSiGD3kGH1^paj&Faahf9ad2HUyJ33T(In1=YWO;r~B6|Nei&M`5~W6fg?>jsmRM zbej#-aQ5qYxiaBv+la3bMG$s#g_43qy5m6Ujsq8e7$R@ORB#-PbA@<=@YgQ_PSI$_ PiEhq+`B_Zzk1OyUr^#ER literal 0 HcmV?d00001 diff --git a/examples/partitioned-heat-conduction-direct/.DS_Store b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/.DS_Store similarity index 84% rename from examples/partitioned-heat-conduction-direct/.DS_Store rename to examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/.DS_Store index 03d268db68702b5b0c43d094f4be83e377cb9a85..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6 100644 GIT binary patch delta 68 zcmZoMXfc=|&Zs)EP;8=}A_oHyFfuR*Y`kd3KJh`-W_At%4o20D7r!%4<`+>E1WGX^ QfYbm1h~2Q+QRFZ)066Lo1poj5 delta 155 zcmZoMXfc=|&e%RNQEZ}~q9_vs0|O%ig8&0VK7%ep4nrzK8bis%M0Hn?1Q$a(LncEp zLoTvpN^x>dQht68P!|I?ih?4b;&g@#h7!1vjfv*$lMO^zHgj|Ea4@!SO#IF~nO{Uv T5M;Ol5GMdJ7;JVFIm`?I=P4jO diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/clean.sh b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/clean.sh new file mode 100755 index 0000000..ad9a90c --- /dev/null +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/clean.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +# Cleaning script for partitioned-heat-conduction 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/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 new file mode 100755 index 0000000..d201a7f --- /dev/null +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/run.sh @@ -0,0 +1,4 @@ +#!/bin/bash +set -e -u + +../gismo-executable/gismo-executable-dirichlet -c ../precice-config.xml diff --git a/examples/partitioned-heat-conduction-vertex-vertex/.DS_Store b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/gismo-executable/.DS_Store similarity index 84% rename from examples/partitioned-heat-conduction-vertex-vertex/.DS_Store rename to examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/gismo-executable/.DS_Store index d84888521f894e3f1a2e8f2936336e4ff7f87f34..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6 100644 GIT binary patch delta 67 zcmZoMXfc=|&Zs)EP;8=}A_oHyFfuR*Y}^>eKJh@*W_At%4o20D8^1G8<`+@q1WGX^ TfYeMj;Zfe4AhLvcVgm~RE=&+7 delta 124 zcmZoMXfc=|&e%3FQEZ}~q9`K+0|O8XFfimZ=rZImq%x#2luS%ipKKt)BFDv$&XCDa z%#e#LnNpmbla!yI!?-aqiXEhEGdBkh2V>jD#P7_L`9%~(nSe$r0C55kgN#1_#0;Aq IMfNiT0F-hX2LJ#7 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 new file mode 100755 index 0000000..8c1c480 --- /dev/null +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/gismo-executable/creat_symlink.sh @@ -0,0 +1,48 @@ +#!/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-IGA-IGA/.DS_Store b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/.DS_Store similarity index 77% rename from examples/partitioned-heat-conduction-IGA-IGA/.DS_Store rename to examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/.DS_Store index adfa1e4353cddddfb3beed7498a6cc3f165735f3..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6 100644 GIT binary patch delta 70 zcmZoMXfc=|#>AjHu~2NHo+1YW5HK<@2y8B7p2o7dfmw=qGdl-A2T%b}*mB z03loQy^GJj^SxAYO+;?GDJMi@A{s&&dm{`H;dR!I^we?+bYYKz7F5wAy_D;Yw;BE- z1N`nrRMHbYLtF2!sB-TsPC2q%%_p_s364n=Q^iIQD-b84KnZo) zVmJwhJyE|@u@RJTa@u@2J=tl8;=<|JKMC&SQb8My0b^jFfxbMBx&GfA@BjCM?8z7~ z2L2TTZjjBgDXyfowR1VHwF&ed%EEq);5G!4P>SKprT79G1ong{z*Mmjgau-M1Og2< J7z01bz!&qqThagk diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/clean.sh b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/clean.sh new file mode 100755 index 0000000..ad9a90c --- /dev/null +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/clean.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +# Cleaning script for partitioned-heat-conduction 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/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 new file mode 100755 index 0000000..2a424ec --- /dev/null +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/run.sh @@ -0,0 +1,4 @@ +#!/bin/bash +set -e -u + +../gismo-executable/gismo-executable-neumann -c ../precice-config.xml diff --git a/examples/partitioned-heat-conduction-IGA-IGA/precice-config.xml b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/precice-config.xml similarity index 100% rename from examples/partitioned-heat-conduction-IGA-IGA/precice-config.xml rename to examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/precice-config.xml diff --git a/examples/perpendicular-flap-vertex-vertex/.DS_Store b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/.DS_Store similarity index 76% rename from examples/perpendicular-flap-vertex-vertex/.DS_Store rename to examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/.DS_Store index 5278b6dbf8c65eb335908e7dde50b0a5e78725ae..71bf82071030689a72d67b16677410902ef7aa10 100644 GIT binary patch delta 256 zcmZoMXfc=|#>B`mu~2NHo}wrV0|Nsi1A_nqLo!1KLk>eKLkWZK#KPtEAPIhk6oyQO zA|OVRWJm|{iy3km@>7bFbCUA&bAb8@DJ=^w%FD^mONZD7v<_|vgD%j{3Ls8qC*}M#SK>1Rjm5D%B9#98dH{1!EKQfE4Y+|us+RV`PM@6IwFut*)UEU!Y2sN`jG`*>e| zJ`G#14#f{6ckit?og(FckGKRa9Tad*JuV$-_jsPC`g)zsHPv44Q+nTfccjGcY+?); z1IECK0p83Z$yCrrW55_N296By{@|gEsbV7-w+;-k1ptm=j)FP&5**_dQ^iIQcOXte zffDMp#c&c1yH~$du@RJTa@u@2UD;`e;^OLfzK`zYQb8My0b}5jfh)P3@%jJU-~V3* z*_AP14E!qw+&EihOQfV{Yb!ZEYd!Q0%EEq);4uV~cof5zkK#LM6xcmafT?052n)pi N2t*oeFa~~=fiGr$Z4>|i diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/dirichlet-gismo/clean.sh b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/dirichlet-gismo/clean.sh new file mode 100755 index 0000000..ad9a90c --- /dev/null +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/dirichlet-gismo/clean.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +# Cleaning script for partitioned-heat-conduction 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/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 new file mode 100755 index 0000000..2612a8c --- /dev/null +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/dirichlet-gismo/run.sh @@ -0,0 +1,4 @@ +#!/bin/bash +set -e -u + +../gismo-executable/gismo-executable -c ../precice-config.xml --plot -s 0 diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/gismo-executable/.DS_Store b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/gismo-executable/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6 GIT binary patch literal 6148 zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0 $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/clean.sh b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/neumann-gismo/clean.sh new file mode 100755 index 0000000..ad9a90c --- /dev/null +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/neumann-gismo/clean.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +# Cleaning script for partitioned-heat-conduction 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/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 new file mode 100755 index 0000000..6e6bddd --- /dev/null +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/neumann-gismo/run.sh @@ -0,0 +1,4 @@ +#!/bin/bash +set -e -u + +../gismo-executable/gismo-executable -c ../precice-config.xml --plot -s 1 diff --git a/examples/partitioned-heat-conduction-vertex-vertex/precice-config.xml b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/precice-config.xml similarity index 71% rename from examples/partitioned-heat-conduction-vertex-vertex/precice-config.xml rename to examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/precice-config.xml index fcb1124..baab5d5 100644 --- a/examples/partitioned-heat-conduction-vertex-vertex/precice-config.xml +++ b/examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/precice-config.xml @@ -7,8 +7,8 @@ enabled="true" /> - - + + @@ -21,33 +21,27 @@ - + - - - + constraint="consistent" /> - + - - - + constraint="consistent" /> @@ -57,21 +51,31 @@ - + - - - + initialize="true" + substeps="true" /> + + + + + + diff --git a/examples/perpendicular-flap-IGA-IGA/perpendicular-flap-IGA-IGA-fluid.cpp b/examples/perpendicular-flap-IGA-IGA-fluid.cpp similarity index 100% rename from examples/perpendicular-flap-IGA-IGA/perpendicular-flap-IGA-IGA-fluid.cpp rename to examples/perpendicular-flap-IGA-IGA-fluid.cpp diff --git a/examples/perpendicular-flap-IGA-IGA/perpendicular-flap-IGA-IGA-solid.cpp b/examples/perpendicular-flap-IGA-IGA-solid.cpp similarity index 100% rename from examples/perpendicular-flap-IGA-IGA/perpendicular-flap-IGA-IGA-solid.cpp rename to examples/perpendicular-flap-IGA-IGA-solid.cpp diff --git a/examples/perpendicular-flap-IGA-IGA/precice-config.xml b/examples/perpendicular-flap-IGA-IGA/precice-config.xml deleted file mode 100644 index a3c2785..0000000 --- a/examples/perpendicular-flap-IGA-IGA/precice-config.xml +++ /dev/null @@ -1,84 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/perpendicular-flap-IGA-Spine/precice-config.xml b/examples/perpendicular-flap-IGA-Spine/precice-config.xml deleted file mode 100644 index 491fae2..0000000 --- a/examples/perpendicular-flap-IGA-Spine/precice-config.xml +++ /dev/null @@ -1,72 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/perpendicular-flap-IGA-Spine/solid-gismo-IGA.cpp b/examples/perpendicular-flap-IGA-Spine/solid-gismo-IGA.cpp deleted file mode 100644 index 1901261..0000000 --- a/examples/perpendicular-flap-IGA-Spine/solid-gismo-IGA.cpp +++ /dev/null @@ -1,463 +0,0 @@ -/** @file flow-over-heated-plate.cpp - - @brief Heat equation participant for the PreCICE example "flow over heated plate" - - 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 (TU Delft, 2023-...) - - Give knots/control points for the boundary curve -*/ - -#include -#include -#include -#include -// #include - -#include -#include - -#ifdef gsStructuralAnalysis_ENABLED -#include -#include -#include -#include -#include -#include - -#include -#endif - -using namespace gismo; - -template class basis_type> -inline void getKnots(const gsBasis & source, std::vector> & tensorKnots) -{ - if ( const basis_type * basis = dynamic_cast*>(&source) ) - for (index_t d=0; d!=DIM; d++) - tensorKnots[d] = basis->knots(d); -} - -int main(int argc, char *argv[]) -{ - //! [Parse command line] - bool plot = false; - index_t plotmod = 1; - index_t numRefine = 1; - 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)); - - // p-refine - if (numElevate!=0) - patches.degreeElevate(numElevate); - - // h-refine - for (int r =0; r < numRefine; ++r) - patches.uniformRefine(); - - // patches.embed(3); - // Create bases - gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) - - 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 KnotMesh = "KnotMesh"; - std::string ControlPointMesh= "ControlPointMesh"; - std::string ForceMesh = "ForceMesh"; - - std::string ControlPointData= "ControlPointData"; - std::string ForceData = "ForceData"; - - - - // Step 3: initialize the participant - // 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,ForceMesh,ForceData,patches,patches.geoDim(),false); - // Add all BCs - // Coupling interface - // bcInfo.addCondition(0, boundary::north, condition_type::neumann , &g_C); - bcInfo.addCondition(0, couplingInterfaces[0], condition_type::neumann , &g_C, -1, true); - bcInfo.addCondition(0, couplingInterfaces[1], condition_type::neumann , &g_C, -1, true); - bcInfo.addCondition(0, couplingInterfaces[2], condition_type::neumann , &g_C, -1, true); - // bcInfo.addCondition(0, boundary::west, condition_type::neumann , &g_C); - // 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); - - std::vector::uPtr> boundaries(couplingInterfaces.size()); - gsMultiPatch<> boundaryGeometry; - - // Derive the basis for boundaries - - // Finish this for both knot vectors and control points - gsMultiBasis<> boundaryBases; - - size_t vector_size; - size_t N = 0; - - // Pre-calculate total size N to resize derived_knots just once. - for(index_t i= 0; i < couplingInterfaces.size();++i) - { - boundaries[i] = patches.patch(0).boundary(couplingInterfaces[i].side()); // Add boundary coefficients to a vector - auto& basis_temp = boundaries[i]->basis(); // Assuming this returns a pointer to gsBasis - boundaryBases.addBasis(&basis_temp); - N += knotsToVector(boundaryBases.basis(i)).size(); - } - - gsVector<> derived_knots; - derived_knots.resize(N); - - // boundaryBases.basis(i).knots().asMatrix() - - // static_cast>(boundaryBases.basis(i)).knots().asMatrix() - - size_t offset = 0; // Tracks the current fill position in derived_knots. - for (index_t i = 0; i < couplingInterfaces.size(); ++i) - { - // No need to call knotsToVector twice for the same basis. - const gsVector<> temp_knots = knotsToVector(boundaryBases.basis(i)); - derived_knots.segment(offset, temp_knots.size()) = temp_knots; - offset += temp_knots.size(); // Increment offset by the size of the current knot vector. - } - gsDebugVar(derived_knots); - - // Control Points - // In one matrix - // gsMatrix<> derived_control_points; - - size_t totalControlPoints = 0; - size_t controlpoint_offset = 0; - for (index_t i = 0; i < couplingInterfaces.size(); ++i) - { - totalControlPoints += boundaries[i]->coefs().rows(); - } - gsMatrix<> derived_control_points(totalControlPoints, boundaries[0]->coefs().cols()); - derived_control_points.setZero(); - - for (index_t i = 0; i < couplingInterfaces.size(); ++i) - { - const gsMatrix<> temp_control_points = boundaries[i]->coefs(); - for (index_t j = 0; j < temp_control_points.rows(); ++j) - { - derived_control_points.row(controlpoint_offset + j) << temp_control_points.row(j); - } - controlpoint_offset += temp_control_points.rows(); - } - - gsDebugVar(derived_control_points); - - // Step 1: write the meshes to PreCICE - // Step 1a: KnotMesh - // get the knots in a matrix, using the utility function knotsToMatrix - gsMatrix<> knots = knotsToMatrix(bases.basis(0)); // - participant.addMesh(KnotMesh,derived_knots.transpose()); - - // Step 1b: ControlPointMesh - // get the control points, in the format where every column is a control point - gsVector controlPointIDs; // needed for writing - gsMatrix<> controlPoints = derived_control_points; - // gsDebugVar(controlPoints); - - participant.addMesh(ControlPointMesh,derived_control_points.transpose(),controlPointIDs); - - // Step 1c: ForceMesh - // Get the quadrature nodes on the coupling interface - gsOptionList quadOptions = gsAssembler<>::defaultOptions(); - - // Get the quadrature points - gsMatrix<> quadPoints = gsQuadrature::getAllNodes(bases.basis(0),quadOptions,couplingInterfaces); - gsDebugVar(quadPoints.dim()); - participant.addMesh(ForceMesh,quadPoints); - - - // Step 2 (not needed???) - // Needed for direct mesh coupling - gsMatrix<> bbox(2,2); - bbox.col(0).setConstant(-1e300); - bbox.col(1).setConstant(1e300); - bbox.transposeInPlace(); - participant.setMeshAccessRegion(ForceMesh,bbox); - - - real_t precice_dt = participant.initialize(); - - - - // for (index_t i = 0; i < couplingInterfaces.size(); ++i) { - // derived_control_points.push_back(boundaries[i]->coefs().transpose()); - // } - // std::cout << "Control Points: " << derived_control_points[0]<< std::endl; - - - // gsDebugVar(boundaryBases.basis(0)); - // derived_knots.push_back(knotsToVector(boundaryBases.basis(0))); - - // gsDebugVar(derived_knots); - - - // Maybe also make the fluid part as well (for testing). - -//---------------------------------------------------------------------------------------------- - // // 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.options().setInt("MaterialLaw",material_law::hooke); - 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); - - // solution.patch(0).coefs() -= patches.patch(0).coefs();// assuming 1 patch here - 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); - - 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; - assembler.constructSolution(solVector,fixedDofs,solution); - // write heat fluxes to interface - // solution.embed(3); - controlPoints = solution.patch(0).coefs().transpose(); - - participant.writeData(ControlPointMesh,ControlPointData,controlPointIDs,controlPoints); - - // 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) - { - // 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< -#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 loadCase = 0; - bool get_readTime = false; - bool get_writeTime = false; - std::string precice_config; - - gsCmdLine cmd("Flow over heated plate for PreCICE."); - 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("l","loadCase", "Load case: 0=constant load, 1='spring' load", loadCase); - 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; } - - //! [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.0,0.0,0.5,1.0)); - - patches.addAutoBoundaries(); - patches.embed(3); - - // Generate discrete fluid mesh - gsMatrix<> bbox = patches.patch(0).support(); - - int numPoints = 100; - gsMatrix<> pointGrid = gsPointGrid<>(bbox, numPoints); - gsMatrix<> mesh = patches.patch(0).eval(pointGrid); - - // Embed the 2D geometry to 3D - - /* - * Initialize the preCICE participant - * - * - */ - std::string participantName = "Fluid"; - gsPreCICE participant(participantName, precice_config); - - /* - * Data initialization - * - * This participant creates its own mesh, and it writes and reads displacements and stresses on its mesh. - * The follow meshes and data are made available: - * - * - Meshes: - * + FluidMesh: This mesh contains the integration points as mesh vertices - * - * - Data: - * + DisplacementData: This data is defined on the FluidMesh and stores the displacement at the integration points - * + StressData: This data is defined on the FluidMesh and stores pressure/forces at the integration points - */ - std::string FluidMesh = "FluidMesh"; - std::string StressData = "StressData"; - std::string DisplacementData = "DisplacementData"; - - - // Step 1: FluidMesh - // Get the quadrature points - gsVector FluidMeshIDs; - gsDebugVar(mesh.dim()); - participant.addMesh(FluidMesh,mesh,FluidMeshIDs); - gsDebugVar(mesh); - - - // Step 2: initialize the participant - real_t precice_dt = participant.initialize(); - gsDebugVar("Initialize Fluid"); - - /* - * Collect the geometry - * - */ - - real_t t = 0, 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<> StressPointData(3,mesh.cols()); - - // Time integration loop - while (participant.isCouplingOngoing()) - { - if (participant.requiresWritingCheckpoint()) - gsInfo<<"Writing Checkpoint\n"; - - // Read control point displacements - gsMatrix<> meshPointDisplacements; - participant.readData(FluidMesh,DisplacementData,FluidMeshIDs,meshPointDisplacements); - - if (get_readTime) - t_read += participant.readTime(); - - - if (loadCase==0) - { - // Write data at the quadrature points - StressPointData.setZero(); - StressPointData.row(2).setConstant(-5e3); - } - else if (loadCase==1) - { - StressPointData.setZero(); - for (index_t k = 0; k!=meshPointDisplacements.cols(); k++) - { - StressPointData(2,k) = -meshPointDisplacements(2,k); - } - // impact loading at t==0 - if (t==0) - StressPointData.row(2).array() += -5e3; - } - else - GISMO_ERROR("Load case "< 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"); - } - - // solution.patch(0).eval_into(points,pointDataMatrix); - // otherDataMatrix< - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/perpendicular-flap-vertex-vertex/solid-gismo/perpendicular-flap-solid-gismo.cpp b/examples/perpendicular-flap-vertex-vertex/solid-gismo/perpendicular-flap-solid-gismo.cpp deleted file mode 100644 index 6d12c4b..0000000 --- a/examples/perpendicular-flap-vertex-vertex/solid-gismo/perpendicular-flap-solid-gismo.cpp +++ /dev/null @@ -1,564 +0,0 @@ -/** @file flow-over-heated-plate.cpp - - @brief Heat equation participant for the PreCICE example "flow over heated plate" - - 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/. - - Solid solver: quadrature point - Fluid solver: displacements on its control points - - Author(s): J.Li (TU Delft, 2023-...) -*/ - -#include -#include -#include -#include -#include -// #include - -#include -#include - - -#ifdef gsKLShell_ENABLED -#include -#include -#include -#include -#include -#endif - -#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 numRefine = 2; - index_t numElevate = 1; - std::string precice_config; - int method = 3; // 1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe, 5: Wilson, 6 RK4 - - std::string dirname = "./output"; - - - 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",method); - 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; } - - //! [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)); - - // Embed the 2D geometry to 3D - gsMultiPatch<> solutions; - patches.addAutoBoundaries(); - patches.embed(3); - // patches.embed(3); - // source function, rhs - gsConstantFunction<> g(0.,0.,3); - - // source function, rhs - gsConstantFunction<> gravity(0.,0.,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"; - - real_t rho = 3000; - real_t E = 4e6; - real_t nu = 0.3; - // real_t nu = 0.3; - real_t mu = E / (2.0 * (1.0 + nu)); - real_t lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); - - /* - * Initialize the preCICE participant - * - * - */ - std::string participantName = "Solid"; - gsPreCICE participant(participantName, precice_config); - - // ---------------------------------------------------------------------------------------------- - // typedef gsExprAssembler<>::geometryMap geometryMap; - // typedef gsExprAssembler<>::space space; - // typedef gsExprAssembler<>::solution solution; - - // gsExprAssembler<> A_quad(1,1); - - // gsInfo<<"Active options:\n"<< A.options() <<"\n"; - - // A.setIntegrationElements(bases); - - // gsExprEvaluator<> ev(A); - - - - - - /* - * Data initialization - * - * This participant creates its own mesh, and it writes and reads displacements and stresses on its mesh. - * The follow meshes and data are made available: - * - * - Meshes: - * + SolidMesh: This mesh contains the integration points as mesh vertices - * - * - Data: - * + DisplacementData: This data is defined on the SolidMesh and stores the displacement at the integration points - * + StressData: This data is defined on the SolidMesh and stores pressure/forces at the integration points - */ - std::string SolidMesh = "Solid-Mesh"; - std::string StressData = "Force"; - std::string DisplacementData = "Displacement"; - // std::string ForceMesh = "Fluid-Mesh"; - - std::vector couplingInterfaces(3); - couplingInterfaces[0] = patchSide(0,boundary::west); - couplingInterfaces[1] = patchSide(0,boundary::north); - couplingInterfaces[2] = patchSide(0,boundary::east); - - // gsOptionList quadOptions = A_quad.options(); - - - // Step 1: SolidMesh - // Get the quadrature nodes on the coupling interface - gsOptionList quadOptions = gsExprAssembler<>::defaultOptions(); - - // Get the quadrature points - gsVector quadPointIDs; - // gsMatrix<> quadPoints = gsQuadrature::getAllNodes(bases.basis(0),quadOptions,couplingInterfaces); - // gsMatrix<> datapoints = patches.patch(0).eval(quadPoints); - // gsMatrix<> quadPointsUpdate(2,13); - - // Set the first point to (0, 0) - // quadPointsUpdate(0, 0) = 0; - // quadPointsUpdate(1, 0) = 0; - - // index_t quadSize = 0; - - // for (std::vector::const_iterator it = couplingInterfaces.begin(); it!=couplingInterfaces.end(); it++) - // { - // // Get a domain iterator on the coupling interface - // typename gsBasis::domainIter domIt = bases.basis(it->patch).makeDomainIterator(it->side()); - - // // First obtain the size of all quadrature points - // typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT - // for (; domIt->good(); domIt->next() ) - // { - // QuRule = gsQuadrature::getPtr(bases.basis(it->patch), quadOptions,it->side().direction()); - // quadSize+=QuRule->numNodes(); - // } - // } - // quadPointsUpdate.block(0, 1, 2, 12) = datapoints.block(0, 0, 2, 12); - - // gsDebugVar(quadPointsUpdate); - - - - - // Initialize parametric coordinates - gsMatrix<> quadPoints = gsQuadrature::getAllNodes(bases.basis(0),quadOptions); - - - // gsMatrix<> quadPoints(patches.domainDim(),quadSize); - // quadPoints.setZero(); - // // Initialize physical coordinates - // gsMatrix<> datapoints(patches.targetDim(),quadSize); - // datapoints.setZero(); - - // Grab all quadrature points - index_t offset = 0; - // Set the dimension of the points - gsMatrix<> nodes; - // Start iteration over elements - gsVector<> tmp; - - - // for (std::vector::const_iterator it = couplingInterfaces.begin(); it!=couplingInterfaces.end(); it++) - // { - // // Get a domain iterator on the coupling interface - // typename gsBasis::domainIter domIt = bases.basis(it->patch).makeDomainIterator(it->side()); - // typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT - // for (domIt->reset(); domIt->good(); domIt->next()) - // { - // QuRule = gsQuadrature::getPtr(bases.basis(it->patch), quadOptions,it->side().direction()); - // // Map the Quadrature rule to the element - // QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(), - // nodes, tmp); - // quadPoints.block(0,offset,patches.domainDim(),QuRule->numNodes()) = nodes; - - // gsMatrix<> tmp2; - // patches.patch(it->patch).eval_into(nodes,tmp2); - // datapoints.block(0,offset,patches.targetDim(),QuRule->numNodes()) = patches.patch(it->patch).eval(nodes); - // offset += QuRule->numNodes(); - // } - // } - - - gsDebugVar(quadPoints); - - - - participant.addMesh(SolidMesh,quadPoints,quadPointIDs); //Set the vertices to be datapoints (quadpoints in physical domain) - - - - gsMatrix<> quadPointsData(3, quadPoints.cols()); - - - // gsDebugVar(quadPoints); - // gsDebugVar(quadPoints.dim()); - quadPointsData.setZero(); - - - - // Step 2: initialize the participant - real_t precice_dt = participant.initialize(); - - gsMultiPatch<> forceMesh; - - gsMultiBasis<> forceBases; - std::vector::uPtr> couplingBoundaries(couplingInterfaces.size()); - couplingBoundaries[0] = patches.patch(0).boundary(couplingInterfaces[0].side()); - couplingBoundaries[1] = patches.patch(0).boundary(couplingInterfaces[1].side()); - couplingBoundaries[2] = patches.patch(0).boundary(couplingInterfaces[2].side()); - - - //--------------------------Option: using spline to represent the stress field------------------- - forceMesh.addPatch(patches.patch(0).boundary(couplingInterfaces[0].side())); - forceMesh.addPatch(patches.patch(0).boundary(couplingInterfaces[1].side())); - forceMesh.addPatch(patches.patch(0).boundary(couplingInterfaces[2].side())); - - gsDebugVar(forceMesh.basis(0).size()); - gsDebugVar(forceMesh.basis(1).size()); - - index_t N_cps = forceMesh.basis(0).size(); //Number of control points in y direction - index_t M_cps = forceMesh.basis(1).size(); //Number of control points in y direction - - - - - // gsDebugVar(forceMesh.basis(0).eval(quadPoints)); - - - //Collect the geometry - - // participant.getMeshVertexIDsAndCoordinates(ForceMesh,quadPointIDs,quadPoints); - - - - -//---------------------------------------------------------------------------------------------- - - // Define boundary conditions - gsBoundaryConditions<> bcInfo; - - // Bottom side - bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, nullptr, -1); - bcInfo.addCondition(0, boundary::south, condition_type::clamped, nullptr, 2); - - // // Top side - // bcInfo.addCondition(0, boundary::west, condition_type::neumann, &forceMesh.patch(0)); - // bcInfo.addCondition(0, boundary::north, condition_type::neumann, &forceMesh.patch(1)); - // bcInfo.addCondition(0, boundary::east, condition_type::neumann, &forceMesh.patch(2)); - - // Assign geometry map - bcInfo.setGeoMap(patches); - gsLookupFunction surfForce(quadPoints, quadPointsData); - - gsDebugVar(quadPoints); - gsDebugVar(quadPointsData); - - // Set up the material matrices - gsFunctionExpr<> E_modulus(std::to_string(E),3); - gsFunctionExpr<> PoissonRatio(std::to_string(nu),3); - gsFunctionExpr<> Density(std::to_string(rho),3); - - //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; - - gsOptionList options; - - 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); - - - gsThinShellAssembler<3, real_t, false> assembler(patches, bases, bcInfo, surfForce, materialMatrix); - gsOptionList assemblerOptions = quadOptions.wrapIntoGroup("ExprAssembler"); - 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(); - gsInfo << "Got here \n"; - assembler.assemble(); - gsSparseMatrix<> K = assembler.matrix(); - - - // Define the solution collection for Paraview - gsFileManager::mkdir(dirname); - gsParaviewCollection collection(dirname + "/solution"); - - // Time step - real_t dt = precice_dt; - - 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; - real_t t_read = 0; - real_t t_write = 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("Got here"); - - // 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|| = "< tempData(2,quadPointsData.cols()); - dt = std::min(dt,precice_dt); - participant.readData(SolidMesh,StressData, quadPointIDs, tempData); - - quadPointsData.row(0) = tempData.row(0).array() ; - quadPointsData.row(1) = tempData.row(1).array() ; - - gsDebugVar(quadPointsData); - 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 - - gsMultiPatch<> solution; - gsVector<> displacements = U; - solution = assembler.constructDisplacement(displacements); - // Information to fluid needs to be 2D - gsMatrix<> tempDisplacement(2, quadPoints.cols()); - gsMatrix<> quadPointsDisplacements = solution.patch(0).eval(quadPoints); - tempDisplacement.row(0) = quadPointsDisplacements.row(0); - tempDisplacement.row(1) = quadPointsDisplacements.row(1); - gsDebugVar(quadPointsDisplacements); - participant.writeData(SolidMesh,DisplacementData,quadPointIDs,tempDisplacement); - - 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 - // 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); - if (get_readTime) - t_read += participant.readTime(); - otherDataMatrix< -#include -#include -//#include -#include -// #include - -#include -#include - - -#ifdef gsKLShell_ENABLED -#include -#include -#include -#include -#include -#endif - -#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 numRefine = 1; - 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 - - std::string dirname = "./output"; - - - 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",method); - 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; } - - //! [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.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"; - - real_t rho = 3000; - real_t E = 4e6; - real_t nu = 0.3; - // real_t nu = 0.3; - // real_t mu = E / (2.0 * (1.0 + nu)); - // real_t lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); - - /* - * Initialize the preCICE participant - * - * - */ - std::string participantName = "Solid"; - gsPreCICE participant(participantName, precice_config); - - - - - /* - * Data initialization - * - * This participant creates its own mesh, and it writes and reads displacements and stresses on its mesh. - * The follow meshes and data are made available: - * - * - Meshes: - * + SolidMesh: This mesh contains the integration points as mesh vertices - * - * - Data: - * + DisplacementData: This data is defined on the SolidMesh and stores the displacement at the integration points - * + StressData: This data is defined on the SolidMesh and stores pressure/forces at the integration points - */ - std::string SolidMesh = "SolidMesh"; - std::string StressData = "StressData"; - std::string DisplacementData = "DisplacementData"; - - - // Step 1: SolidMesh - // Get the quadrature nodes on the coupling interface - gsOptionList quadOptions = gsExprAssembler<>::defaultOptions(); - - // Get the quadrature points - gsVector quadPointIDs; - gsMatrix<> quadPoints = gsQuadrature::getAllNodes(bases.basis(0),quadOptions); - participant.addMesh(SolidMesh,quadPoints,quadPointIDs); - gsMatrix<> quadPointsData(3, quadPoints.cols()); - quadPointsData.setZero(); - - - // Step 2: initialize the participant - real_t precice_dt = participant.initialize(); - -//---------------------------------------------------------------------------------------------- - - // Define boundary conditions - 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 - // gsPreCICEFunction surfForce(&participant,SolidMesh,StressData,patches,patches.parDim(),patches.geoDim(),false); - gsLookupFunction surfForce(quadPoints, quadPointsData); - - - - -//------------------------------------------IF gsElasticity is Used------------------------------ - // // 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.options().setInt("MaterialLaw",material_law::hooke); - // assembler.assemble(); - - - // gsMatrix Minv; - // gsSparseMatrix<> M = massAssembler.matrix(); - // gsSparseMatrix<> K = assembler.matrix(); - // gsSparseMatrix<> K_T; -//------------------------------------------IF gsElasticity is Used------------------------------ - - - // Set up the material matrices - gsFunctionExpr<> E_modulus(std::to_string(E),3); - gsFunctionExpr<> PoissonRatio(std::to_string(nu),3); - gsFunctionExpr<> Density(std::to_string(rho),3); - - //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; - - gsOptionList options; - - 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); - - gsThinShellAssembler<3, real_t, true> assembler(patches, bases, bcInfo, surfForce, materialMatrix); - gsOptionList assemblerOptions = quadOptions.wrapIntoGroup("ExprAssembler"); - 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 - gsFileManager::mkdir(dirname); - gsParaviewCollection collection(dirname + "/solution"); - - // Time step - real_t dt = precice_dt; - - 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; - real_t t_read = 0; - real_t t_write = 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); - - // 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); - gsMatrix<> quadPointsDisplacements = solution.patch(0).eval(quadPoints); - participant.writeData(SolidMesh,DisplacementData,quadPointIDs,quadPointsDisplacements); - - 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 - // 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); - if (get_readTime) - t_read += participant.readTime(); - otherDataMatrix< class gsLookupFunction : public gsFunction @@ -95,11 +92,26 @@ class gsLookupFunction : public gsFunction virtual void eval_into(const gsMatrix& u, gsMatrix& result) const { index_t col; + gsDebugVar(u); result.resize(this->targetDim(),u.cols()); result.setZero(); + + // // 打印 m_map 的内容 + // gsInfo << "Content of m_map in eval_into:\n"; + // // Print m_map content without structured bindings + // for (auto it = m_map.begin(); it != m_map.end(); ++it) + // { + // gsInfo << "Key: "; + // for (int i = 0; i < it->first.size(); ++i) + // { + // gsInfo << it->first[i] << " "; + // } + // gsInfo << "=> Value: " << it->second << "\n"; + // } + for (index_t k = 0; k!= u.cols(); k++) { - + // gsDebugVar(u.col(k)); GISMO_ASSERT(m_map.find(u.col(k))!=m_map.end(),"Coordinate " + std::to_string(k) + " not registered in the table"); col = m_map.at(u.col(k)); result.col(k) = m_data.col(col); @@ -137,6 +149,14 @@ class gsLookupFunction : public gsFunction GISMO_ASSERT(m_points.cols()==m_data.cols(),"Points and data must have the same number of columns"); for (index_t k = 0; k != m_points.cols(); k++) m_map.insert({m_points.col(k),k}); // m_map.at(vector) returns the column index of vector + gsInfo << "Content of m_map:\n"; + for (const auto& [key, value] : m_map) { + gsInfo << "Key: "; + for (int i = 0; i < key.size(); ++i) { + gsInfo << key[i] << " "; + } + gsInfo << "=> Value: " << value << "\n"; + } } /// See \a gsFunction diff --git a/src/gsPreCICE.h b/gsPreCICE.h similarity index 87% rename from src/gsPreCICE.h rename to gsPreCICE.h index 127c0a5..a932158 100644 --- a/src/gsPreCICE.h +++ b/gsPreCICE.h @@ -8,46 +8,7 @@ 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): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) - - @details Summary of Functions: - 1. Constructors: - - gsPreCICE() DEFAULT - - gsPreCICE(std::string participantName, std::string configurationFileName): - Initializes the solver interface with custom participant and configuration names. - - 2. Coupling Status: - - isCouplingOngoing(): Checks if coupling between solvers is still ongoing. - - isTimeWindowComplete(): Checks if the time window for communication is complete. - - requiresInitialData(): Checks if initial data is required (preCICE v3.0.0). - - requiresReadingCheckpoint(): Checks if checkpoint reading is required. - - requiresWritingCheckpoint(): Checks if checkpoint writing is required. - - 3. Initialization: - - initialize(): Initializes the solver interface, sets time variables, and returns the maximum timestep size. - - addMesh(): Adds a mesh to the preCICE interface by providing a mesh name and a matrix of points - (optionally returns vertex IDs). - - 4. Data Handling - - readData(): Reads scalar data from a mesh by coordinates or vertex IDs. - - writeData(): Writes scalar data to a mesh by coordinates or vertex IDs. - - 5. Mesh Handling - - getMeshVertexIDsFromPositions(): Retrieves mesh vertex IDs from positions, - ensuring mesh point compatibility. - - setMeshAccessRegion(): Sets an access region for a mesh using bounding boxes. - - getMeshVertexIDsAndCoordinates(): Retrieves the vertex IDs and coordinates for a given mesh name. - - getMeshID(): Returns the mesh ID based on a mesh name. - - 6. Timing Information - - readTime(): Returns the total time spent reading data. - - writeTime(): Returns the total time spent writing data. - - initializeTime(): Returns the total time spent during initialization. - - 7. Advance & Finalize: - - advance(T dt): Advances the solver interface by a timestep. - - finalize(): Finalizes the preCICE solver interface when the coupling is complete. - + Author(s): H.M. Verhelst (TU Delft, 2019-2024), J.Li (TU Delft, 2023 - ...) */ #pragma once @@ -165,7 +126,6 @@ struct mapCompare * @param[in] meshName The mesh name * @param[in] points The points stored in columns */ - void addMesh(const std::string & meshName, const gsMatrix & points, gsVector & vertexIDs) { const index_t dMesh = m_interface.getMeshDimensions(meshName); diff --git a/src/gsPreCICEFunction.h b/gsPreCICEFunction.h similarity index 97% rename from src/gsPreCICEFunction.h rename to gsPreCICEFunction.h index bdfcb94..1ea5f8c 100644 --- a/src/gsPreCICEFunction.h +++ b/gsPreCICEFunction.h @@ -8,7 +8,7 @@ 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): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) + Author(s): H.M. Verhelst (TU Delft, 2019-2024), J.Li (TU Delft, 2023 - ...) */ #pragma once @@ -22,7 +22,7 @@ namespace gismo /** * @brief Class defining a gsFunction that reads from the precice::SolverInterface - * @usage Combine with gsElasticity module to update boundary conditions. + * * @tparam T Number format */ template diff --git a/src/gsPreCICEUtils.h b/gsPreCICEUtils.h similarity index 89% rename from src/gsPreCICEUtils.h rename to gsPreCICEUtils.h index 391265f..9f9c2b5 100644 --- a/src/gsPreCICEUtils.h +++ b/gsPreCICEUtils.h @@ -8,7 +8,7 @@ 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): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) + Author(s): H.M. Verhelst (TU Delft, 2019-2024), J.Li (TU Delft, 2023 - ...) */ #pragma once @@ -40,6 +40,7 @@ inline void getKnots(const gsBasis & source, std::vector> & t * * @return { description_of_the_return_value } */ + template gsVector knotsToVector(const gsBasis & basis) { @@ -170,6 +171,16 @@ gsMatrix knotVectorUnpack(const gsMatrix & knots, index_t numBoundaries) // return std::make_pair(knotMatrix, coefMatrix); // } +/** + * @brief Pack the knot and ratio matrices of a gsMultiPatch object into a single matrix. + * + * @param[in] mp The gsMultiPatch object + * + * @tparam T { description } + * + * @return A tuple containing the packed knot matrix, the packed ratio matrix, the number of columns in each patch's knot matrix, and the number of columns in each patch's ratio matrix. + */ + template std::tuple, gsMatrix, std::vector, std::vector> packMultiPatch(const gsMultiPatch &mp) { std::vector> knotMatrices; @@ -233,7 +244,7 @@ std::tuple, gsMatrix, std::vector, std::vector> * * @tparam T { description } * - * @return { description_of_the_return_value } + * @return A gsMultiPatch object. */ @@ -263,7 +274,17 @@ gsMultiPatch unpackMultiPatch(const gsMatrix &knotMatrix, const gsMatrix gsMatrix unPackControlPoints(const gsMatrix & controlPoints, const gsMatrix & kv_unpacked, index_t knot_index, index_t cp_index) @@ -296,8 +317,13 @@ gsMatrix unPackControlPoints(const gsMatrix & controlPoints, const gsMatri return unpackedCps; } - - +/** + * @brief Convert a matrix of knot vectors into a gsBasis object. + * + * @param[in] knots The matrix of knot vectors + * + * @return A shared pointer to the gsBasis object. + */ template typename gsBasis::Ptr knotMatrixToBasis(const gsMatrix & knots) @@ -313,10 +339,10 @@ typename gsBasis::Ptr knotMatrixToBasis(const gsMatrix & knots) std::copy_if(knots.row(d).begin(),knots.row(d).end(), std::back_inserter(tmp), [](T a){return !math::isnan(a);}); - KVs[d] = gsKnotVector(tmp);/////// MEMLEAK - // KVs[d] = std::move(gsKnotVector(tmp)); + // KVs[d] = gsKnotVector(tmp);/////// MEMLEAK + KVs[d] = std::move(gsKnotVector(tmp)); gsDebugVar(d); - gsDebugVar(KVs[d].size); + gsDebugVar(KVs[d]); gsDebug<<"(gsPreCICEUtils.h: There is a memory leak in the line above)\n"; } @@ -336,5 +362,3 @@ typename gsBasis::Ptr knotMatrixToBasis(const gsMatrix & knots) } } //namespace gismo - - From 51eaa1fea734120c10bcee4554f3cbaf70441026 Mon Sep 17 00:00:00 2001 From: Crazy-Rich-Meghan Date: Thu, 5 Dec 2024 17:19:47 +0100 Subject: [PATCH 02/14] Update documentation for the adapter function --- gsLookupFunction.h | 24 ++++++++---------------- gsPreCICE.h | 33 +++++++++++++++++++++++++++++++++ gsPreCICEFunction.h | 2 +- 3 files changed, 42 insertions(+), 17 deletions(-) diff --git a/gsLookupFunction.h b/gsLookupFunction.h index c37b8e3..98e3593 100644 --- a/gsLookupFunction.h +++ b/gsLookupFunction.h @@ -21,10 +21,15 @@ namespace gismo { /** - * @brief Class defining a function that looks up registered data on points - * - * @tparam T Number format + * @brief Class defining a function that looks up registered data on points. + * @usage Combine with gsThinShellAssembler to update the stress on the solid. + * @details The gsLookupFunction enables + * efficient data lookups based on spatial coordinates. When given a set of points and corresponding + * data, it creates a mapping that allows for quick retrieval of the data based on the point + * coordinates. + * @param T Number format */ + template class gsLookupFunction : public gsFunction { @@ -96,19 +101,6 @@ class gsLookupFunction : public gsFunction result.resize(this->targetDim(),u.cols()); result.setZero(); - // // 打印 m_map 的内容 - // gsInfo << "Content of m_map in eval_into:\n"; - // // Print m_map content without structured bindings - // for (auto it = m_map.begin(); it != m_map.end(); ++it) - // { - // gsInfo << "Key: "; - // for (int i = 0; i < it->first.size(); ++i) - // { - // gsInfo << it->first[i] << " "; - // } - // gsInfo << "=> Value: " << it->second << "\n"; - // } - for (index_t k = 0; k!= u.cols(); k++) { // gsDebugVar(u.col(k)); diff --git a/gsPreCICE.h b/gsPreCICE.h index a932158..cdb984d 100644 --- a/gsPreCICE.h +++ b/gsPreCICE.h @@ -9,6 +9,38 @@ file, You can obtain one at http://mozilla.org/MPL/2.0/. Author(s): H.M. Verhelst (TU Delft, 2019-2024), J.Li (TU Delft, 2023 - ...) + + @details Summary of Functions: + 1. Constructors: + - gsPreCICE() DEFAULT + - gsPreCICE(std::string participantName, std::string configurationFileName): + Initializes the solver interface with custom participant and configuration names. + 2. Coupling Status: + - isCouplingOngoing(): Checks if coupling between solvers is still ongoing. + - isTimeWindowComplete(): Checks if the time window for communication is complete. + - requiresInitialData(): Checks if initial data is required (preCICE v3.0.0). + - requiresReadingCheckpoint(): Checks if checkpoint reading is required. + - requiresWritingCheckpoint(): Checks if checkpoint writing is required. + 3. Initialization: + - initialize(): Initializes the solver interface, sets time variables, and returns the maximum timestep size. + - addMesh(): Adds a mesh to the preCICE interface by providing a mesh name and a matrix of points + (optionally returns vertex IDs). + 4. Data Handling + - readData(): Reads scalar data from a mesh by coordinates or vertex IDs. + - writeData(): Writes scalar data to a mesh by coordinates or vertex IDs. + 5. Mesh Handling + - getMeshVertexIDsFromPositions(): Retrieves mesh vertex IDs from positions, + ensuring mesh point compatibility. + - setMeshAccessRegion(): Sets an access region for a mesh using bounding boxes. + - getMeshVertexIDsAndCoordinates(): Retrieves the vertex IDs and coordinates for a given mesh name. + - getMeshID(): Returns the mesh ID based on a mesh name. + 6. Timing Information + - readTime(): Returns the total time spent reading data. + - writeTime(): Returns the total time spent writing data. + - initializeTime(): Returns the total time spent during initialization. + 7. Advance & Finalize: + - advance(T dt): Advances the solver interface by a timestep. + - finalize(): Finalizes the preCICE solver interface when the coupling is complete. */ #pragma once @@ -319,6 +351,7 @@ struct mapCompare * @see getMeshVertexSize() to get the amount of vertices in the mesh * @see getMeshDimensions() to get the spacial dimensionality of the mesh */ + void getMeshVertexIDsAndCoordinates(const std::string & meshName, gsVector & IDs, gsMatrix & coords) const { const index_t nPoints = m_interface.getMeshVertexSize(meshName); diff --git a/gsPreCICEFunction.h b/gsPreCICEFunction.h index 1ea5f8c..f203f4a 100644 --- a/gsPreCICEFunction.h +++ b/gsPreCICEFunction.h @@ -22,7 +22,7 @@ namespace gismo /** * @brief Class defining a gsFunction that reads from the precice::SolverInterface - * + * @details The gsPreCICEFunction enables efficient data lookups based on spatial coordinates. When given a set of points and corresponding * @tparam T Number format */ template From c22c6bb3433ef8b553af130a81b5a07a9034b9fd Mon Sep 17 00:00:00 2001 From: Jingya Li <70752306+Crazy-Rich-Meghan@users.noreply.github.com> Date: Thu, 5 Dec 2024 17:45:04 +0100 Subject: [PATCH 03/14] Update .gitignore to exclude .DS_Store file --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index d2dd0fe..f7fe7bd 100644 --- a/.gitignore +++ b/.gitignore @@ -38,6 +38,8 @@ *.pvd *.vtk +# Files from macOS +.DS_Store # Log files *.log From c254482ff56d35b75f0e106cb4ae62eb4fe0b5d6 Mon Sep 17 00:00:00 2001 From: Jingya Li <70752306+Crazy-Rich-Meghan@users.noreply.github.com> Date: Thu, 5 Dec 2024 17:45:34 +0100 Subject: [PATCH 04/14] Delete examples/partitioned-heat-conduction/.DS_Store --- examples/partitioned-heat-conduction/.DS_Store | Bin 8196 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 examples/partitioned-heat-conduction/.DS_Store diff --git a/examples/partitioned-heat-conduction/.DS_Store b/examples/partitioned-heat-conduction/.DS_Store deleted file mode 100644 index aeb464e46745c2c1ec31f573000d8177c362c00c..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8196 zcmeHMU2hUW6ur};UDBBFphmGm zUE^4xygE?G5CE7%w@j!b4v-vMV_oA|p`>C?6+H+|6*|NaijHxc&4G1|V}*)NLeWX+ zm4%K_gj^jwTbh%oD>S81z$h@Q0LSk4s822>GIswHzx!08L&Whe`b%^mqk5kPs7-82 zpd?R#PDVaA_cKMc}rYwMR-S()o-FdjOYn+?*2o4n`N8*IUw)Lm~s=nUJ1+s!!2+#u>Ck`VYEn0$B_M1DMM#DmCB zCD%~{R^H0D3wKXWwyTeecC}U;744HsrS_y~Z&x3UMtSSbX8Gw}`%Bo1;%`b)p|2&n zrP|Z-J8G9C$3$VC`;w+Sa*`;Fqb_FT4{(2Bza8Na@B!-3F}!!GbzOY81ZhtDa|zR* z^KG8?Xa~<$i*|tzs7HQ3^Ii_pEP>~LkV9W7l^*O-nfCFxhj`d~@Mwy5OI`ymYH!hI zUvuwg*aY;DDyZ**-vrCoj4nn}@MCJ@RH;iMz(vtPn|t{x3DGK}-k)A6F8a#U1@pD! zG;tF&vEMQ*kFkdyY?;NYFLr4y$~$wa&z^DYg=oIMue82%M!_-3C~$=etg5v)IR4+c z`1k)S+@V?8C}0#=wgMv8a#{_nGWimFenF14ZRAhLoY-!xP*PCHa2zPZap3e1L)2}U dGN-O_tPnj&fB%Po#lAO~_rG+PDXxhEzX6h%RJ;HH From 3b5432ad73165b56dce240d2926fa1067a7c6a83 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst <38376601+hverhelst@users.noreply.github.com> Date: Thu, 5 Dec 2024 18:11:24 +0100 Subject: [PATCH 05/14] Update ci.yml --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ea95106..3ee120d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -72,4 +72,4 @@ jobs: -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" \ -D CTEST_SITE="${{ matrix.os }}_[actions]" \ -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON;-Dprecice_DIR=$precice_DIR" \ - -D GISMO_OPTIONAL="${{ github.event.repository.name }}" -Q + -D GISMO_OPTIONAL="gsElasticity\\;${{ github.event.repository.name }}'" -Q From c049332433e9ea9dec464c67668e4996e90e9e45 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst <38376601+hverhelst@users.noreply.github.com> Date: Thu, 5 Dec 2024 18:23:16 +0100 Subject: [PATCH 06/14] Update ci.yml --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3ee120d..744a289 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -72,4 +72,4 @@ jobs: -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" \ -D CTEST_SITE="${{ matrix.os }}_[actions]" \ -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON;-Dprecice_DIR=$precice_DIR" \ - -D GISMO_OPTIONAL="gsElasticity\\;${{ github.event.repository.name }}'" -Q + -D GISMO_OPTIONAL="gsElasticity\\;${{ github.event.repository.name }}" -Q From d0f6d6f5fe180fc3d45036fe993fb7163025b6b5 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst <38376601+hverhelst@users.noreply.github.com> Date: Fri, 6 Dec 2024 09:29:02 +0100 Subject: [PATCH 07/14] Update README.md --- README.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 36cffe9..1d97c6f 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,6 @@ |Dependency|[preCICE v.3](https://github.com/gismo/gsPreCICE)| |Maintainer|[j.li-9@tudelft.nl](mailto:j.li-9@tudelft.nl)| - ## Start here 1. Install [preCICE](https://precice.org/quickstart.html). @@ -50,6 +49,8 @@ cd fluid- **Note**: You need to perform steps 2 and 5 if you want to run the simulation with other libraries. +## Examples +- [Partitioned Heat Conduction](examples/partitioned-heat-conduction/README.md) ## Versions @@ -57,4 +58,4 @@ The latest supported G+Smo version is v24.08.0, and the latest supported preCICE ## References -[1] Benjamin Uekermann, Hans-Joachim Bungartz, Lucia Cheung Yau, Gerasimos Chourdakis and Alexander Rusch. Official preCICE Adapters for Standard Open-Source Solvers. In Proceedings of the _7th GACM Colloquium on Computational Mechanics for Young Scientists from Academia_, 2017. \ No newline at end of file +[1] Benjamin Uekermann, Hans-Joachim Bungartz, Lucia Cheung Yau, Gerasimos Chourdakis and Alexander Rusch. Official preCICE Adapters for Standard Open-Source Solvers. In Proceedings of the _7th GACM Colloquium on Computational Mechanics for Young Scientists from Academia_, 2017. From a377c6e3384bd17769b8ddd205f1ce916dba49f5 Mon Sep 17 00:00:00 2001 From: Jingya Li <70752306+Crazy-Rich-Meghan@users.noreply.github.com> Date: Fri, 6 Dec 2024 09:34:35 +0100 Subject: [PATCH 08/14] Update maintainer in README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 1d97c6f..ec61c31 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ |Repository|[gismo/gismo/gsPreCICE](https://github.com/gismo/gsPreCICE)| |Developer|[Hugo Verhelst](https://github.com/hverhelst), [Jingya Li](https://github.com/Crazy-Rich-Meghan) | |Dependency|[preCICE v.3](https://github.com/gismo/gsPreCICE)| -|Maintainer|[j.li-9@tudelft.nl](mailto:j.li-9@tudelft.nl)| +|Maintainer|[Jingya Li](https://github.com/Crazy-Rich-Meghan), [Hugo Verhelst](https://github.com/hverhelst)| ## Start here From 53d335da0053226d24f7c4ad8348251fac1755f1 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Sun, 8 Dec 2024 19:45:41 +0100 Subject: [PATCH 09/14] add `gsPreCICE-examples` target in CMake add perform CI only on `gsPreCICE-examples` --- .github/workflows/ci.yml | 5 +++-- CMakeLists.txt | 9 ++++++++- examples/CMakeLists.txt | 17 +++++++++++++---- 3 files changed, 24 insertions(+), 7 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 744a289..da870a1 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -33,7 +33,7 @@ jobs: with: path: ./gismo/optional/${{ github.event.repository.name }} #token: ${{ secrets.GH_PAT }} - + # Step to install preCICE - name: Install preCICE if: runner.os == 'Linux' @@ -62,7 +62,7 @@ jobs: make -j 10 sudo make install precice_DIR=$(pwd) - + - name: "Run for ${{ matrix.os }}" shell: bash working-directory: ${{ runner.workspace }} @@ -71,5 +71,6 @@ jobs: ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake \ -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" \ -D CTEST_SITE="${{ matrix.os }}_[actions]" \ + -D LABELS_FOR_SUBPROJECTS="gsPreCICE-examples" \ -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON;-Dprecice_DIR=$precice_DIR" \ -D GISMO_OPTIONAL="gsElasticity\\;${{ github.event.repository.name }}" -Q diff --git a/CMakeLists.txt b/CMakeLists.txt index 4345bb7..422b18c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -53,4 +53,11 @@ set(gismo_UNITTESTS ${gismo_UNITTESTS} ${unittests_SRCS} CACHE INTERNAL "gismo list of unittests") # Look for CMakeLists.txt in examples -add_subdirectory(examples) +# add example files +if(GISMO_BUILD_EXAMPLES) + add_custom_target(${PROJECT_NAME}-examples) + add_subdirectory(examples) +else() + add_subdirectory(examples EXCLUDE_FROM_ALL) +endif(GISMO_BUILD_EXAMPLES) + diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index e55ccfc..213bfe1 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -11,8 +11,11 @@ aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR} FILES) foreach(file ${FILES}) add_gismo_executable(${file}) get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") - set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") + if(GISMO_TEST_PRECICE) + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}-examples") + endif() + set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") + add_dependencies(${PROJECT_NAME}-examples ${tarname}) # install the example executables (optionally) install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) endforeach(file ${FILES}) @@ -24,8 +27,11 @@ FOREACH(subdir ${SUBDIRS}) foreach(file ${FILES}) add_gismo_executable(${file}) get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + if(GISMO_TEST_PRECICE) + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + endif() set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") + add_dependencies(${PROJECT_NAME}-examples ${tarname}) # install the example executables (optionally) install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) endforeach(file ${FILES}) @@ -40,8 +46,11 @@ FOREACH(subdir ${SUBDIRS}) foreach(file ${FILES}) add_gismo_executable(${file}) get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + if(GISMO_TEST_PRECICE) + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + endif() set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") + add_dependencies(${PROJECT_NAME}-examples ${tarname}) # install the example executables (optionally) install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) endforeach(file ${FILES}) From 1149f37da8a767d3ea91911791de70a600e5b254 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Sun, 8 Dec 2024 19:58:17 +0100 Subject: [PATCH 10/14] trigger new CI --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index da870a1..c8ad725 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -71,6 +71,6 @@ jobs: ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake \ -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" \ -D CTEST_SITE="${{ matrix.os }}_[actions]" \ - -D LABELS_FOR_SUBPROJECTS="gsPreCICE-examples" \ + -D LABELS_FOR_SUBPROJECTS="gsPreCICE-examples" \ -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON;-Dprecice_DIR=$precice_DIR" \ -D GISMO_OPTIONAL="gsElasticity\\;${{ github.event.repository.name }}" -Q From 5c3e2502d50697d17e879d789367ec91a86f1a44 Mon Sep 17 00:00:00 2001 From: Crazy-Rich-Meghan Date: Mon, 9 Dec 2024 20:39:13 +0100 Subject: [PATCH 11/14] Resolve reviews --- #CMakeLists.txt# | 86 ++++ .gitignore | 2 + examples/flow-over-heated-plate.cpp | 316 ++++++++++++++ .../partitioned-heat-conduction-direct.cpp | 407 ++++++++++++++++++ .../partitioned-heat-conduction/README.md | 4 +- gsLookupFunction.h | 25 +- gsPreCICEUtils.h | 20 +- 7 files changed, 835 insertions(+), 25 deletions(-) create mode 100644 #CMakeLists.txt# create mode 100644 examples/flow-over-heated-plate.cpp create mode 100644 examples/partitioned-heat-conduction-direct.cpp diff --git a/#CMakeLists.txt# b/#CMakeLists.txt# new file mode 100644 index 0000000..0280e55 --- /dev/null +++ b/#CMakeLists.txt# @@ -0,0 +1,86 @@ +###################################################################### +## CMakeLists.txt --- gsPreCICE +## This file is part of the G+Smo library. +## +###################################################################### + +if(POLICY CMP0076) +cmake_policy(SET CMP0076 NEW) +endif() + + +## gsPreCICE module +project(gsPreCICE) + +# Apply G+Smo config +include(gsConfig) + + +## Collect files +aux_header_directory (${CMAKE_CURRENT_SOURCE_DIR}/src ${PROJECT_NAME}_H ) + + + + +if( (NOT GISMO_BUILD_LIB) ) + aux_instance_directory (${CMAKE_CURRENT_SOURCE_DIR} ${PROJECT_NAME}_INS) + if(${PROJECT_NAME}_INS) + LIST( REMOVE_ITEM ${PROJECT_NAME}_CPP ${${PROJECT_NAME}_INS}) + endif() +endif() + +# Add object library +add_library(${PROJECT_NAME} OBJECT + ${${PROJECT_NAME}_H} + ${${PROJECT_NAME}_HPP} + ${${PROJECT_NAME}_CPP} ) + +set_target_properties(${PROJECT_NAME} PROPERTIES + COMPILE_DEFINITIONS gismo_EXPORTS + POSITION_INDEPENDENT_CODE ON + LINKER_LANGUAGE CXX + FOLDER "G+Smo modules" ) + +set(gismo_MODULES ${gismo_MODULES} $ + CACHE INTERNAL "G+Smo modules" ) + +#Symlink include dir (in case your headers are in /src) +execute_process( COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_CURRENT_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/${PROJECT_NAME}) + + +# Search for preCICE +find_package(precice REQUIRED CONFIG) + +target_link_libraries(${PROJECT_NAME} PRIVATE precice::precice) + +#set_property(TARGET ${PROJECT_NAME} APPEND PROPERTY INTERFACE_LINK_LIBRARIES precice::precice) + +target_include_directories(${PROJECT_NAME} INTERFACE +$ +) + +add_dependencies(${PROJECT_NAME} precice) + + +install(DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/${PROJECT_NAME}" + DESTINATION include/gismo + FILES_MATCHING PATTERN "*.h" ) + +# add filedata folder +add_definitions(-DELAST_DATA_DIR="${CMAKE_CURRENT_SOURCE_DIR}/filedata/") + +# add example files +aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR}/examples FILES) +foreach(file ${FILES}) + add_gismo_executable(${file}) + get_filename_component(tarname ${file} NAME_WE) # name without extension + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") + # install the example executables (optionally) + install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) +endforeach(file ${FILES}) + +# add unittests +aux_gs_cpp_directory(${PROJECT_SOURCE_DIR}/unittests unittests_SRCS) +set(gismo_UNITTESTS ${gismo_UNITTESTS} ${unittests_SRCS} + CACHE INTERNAL "gismo list of unittests") diff --git a/.gitignore b/.gitignore index f7fe7bd..06d74d9 100644 --- a/.gitignore +++ b/.gitignore @@ -39,7 +39,9 @@ *.vtk # Files from macOS +*.DS_Store .DS_Store +.DS_Store? # Log files *.log diff --git a/examples/flow-over-heated-plate.cpp b/examples/flow-over-heated-plate.cpp new file mode 100644 index 0000000..d4cca4e --- /dev/null +++ b/examples/flow-over-heated-plate.cpp @@ -0,0 +1,316 @@ +/** @file flow-over-heated-plate.cpp + + @brief Heat equation participant for the PreCICE example "flow over heated plate" + + 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): H.M. Verhelst +*/ + +#include +#include +#include + +using namespace gismo; + +int main(int argc, char *argv[]) +{ + //! [Parse command line] + bool plot = false; + index_t plotmod = 1; + index_t numRefine = 5; + index_t numElevate = 0; + std::string precice_config("../precice_config.xml"); + + 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); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + //! [Read input file] + + // Generate domain + gsMultiPatch<> patches; + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,-0.25,1.0,0.0)); + + // Set external heat-flux to zero + gsFunctionExpr<> f("0",2); + + // 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"; + + // Set heat conduction coefficient + real_t k_temp = 100; + // Set bottom wall temp + real_t u_wall = 310; + +// ---------------------------------------------------------------------------------------------- + typedef gsExprAssembler<>::geometryMap geometryMap; + typedef gsExprAssembler<>::space space; + typedef gsExprAssembler<>::solution solution; + + gsExprAssembler<> A(1,1); + + gsInfo<<"Active options:\n"<< A.options() <<"\n"; + + A.setIntegrationElements(bases); + + gsExprEvaluator<> ev(A); + +// ---------------------------------------------------------------------------------------------- + // Set the interface for the precice coupling + patchSide couplingInterface(0,boundary::north); + + // Get a domain iterator on the coupling interface + typename gsBasis::domainIter domIt = bases.basis(couplingInterface.patch).makeDomainIterator(couplingInterface.side()); + + // Set the dimension of the points + gsMatrix<> nodes; + // Start iteration over elements + gsVector<> tmp; + + gsOptionList quadOptions = A.options(); + + // First obtain the size of all quadrature points + index_t quadSize = 0; + typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT + for (; domIt->good(); domIt->next() ) + { + QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); + quadSize+=QuRule->numNodes(); + } + + // Initialize parametric coordinates + gsMatrix<> uv(patches.domainDim(),quadSize); + // Initialize physical coordinates + gsMatrix<> xy(patches.targetDim(),quadSize); + + // Grab all quadrature points + index_t offset = 0; + for (domIt->reset(); domIt->good(); domIt->next()) + { + QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); + // Map the Quadrature rule to the element + QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(), + nodes, tmp); + uv.block(0,offset,patches.domainDim(),QuRule->numNodes()) = nodes; + + gsMatrix<> tmp2; + patches.patch(couplingInterface.patch).eval_into(nodes,tmp2); + xy.block(0,offset,patches.targetDim(),QuRule->numNodes()) = patches.patch(couplingInterface.patch).eval(nodes); + offset += QuRule->numNodes(); + } + + // Define precice interface + gsPreCICE interface("Solid", precice_config); + interface.addMesh("Solid-Mesh",xy); + real_t precice_dt = interface.initialize(); + + index_t meshID = interface.getMeshID("Solid-Mesh"); + index_t tempID = interface.getDataID("Temperature",meshID); + index_t fluxID = interface.getDataID("Heat-Flux",meshID); + +// ---------------------------------------------------------------------------------------------- + + // Define boundary conditions + gsBoundaryConditions<> bcInfo; + // Dirichlet side + gsConstantFunction<> g_D(u_wall,2); + // Homogeneous Neumann + gsConstantFunction<> g_N(0.0,2); + // Coupling side + gsPreCICEFunction g_C(&interface,meshID,tempID,patches); + // Add all BCs + // Isolated Neumann sides + bcInfo.addCondition(0, boundary::east, condition_type::neumann , &g_N, 0, false, 0); + bcInfo.addCondition(0, boundary::west, condition_type::neumann , &g_N, 0, false, 0); + // Bottom side (prescribed temp) + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &g_D, 0, false, 0); + // Coupling interface + bcInfo.addCondition(0, boundary::north, condition_type::dirichlet, &g_C, 0, false, 0); + // Assign geometry map + bcInfo.setGeoMap(patches); + +// ---------------------------------------------------------------------------------------------- + + // Time integration coefficient (0.0 = explicit, 1.0 = implicit) + real_t theta = 1.0; + + // Set the geometry map + geometryMap G = A.getMap(patches); + + // Set the discretization space + space u = A.getSpace(bases); + u.setup(bcInfo, dirichlet::homogeneous, 0); + + // Set the source term + auto ff = A.getCoeff(f, G); + + // Set the solution + gsMatrix<> solVector; + solution u_sol = A.getSolution(u, solVector); + + // Assemble mass matrix + A.initSystem(); + A.assemble( u * u.tr() * meas(G)); + gsSparseMatrix<> M = A.matrix(); + + // Assemble stiffness matrix (NOTE: also adds the dirichlet BCs inside the matrixs) + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G) ); + gsSparseMatrix<> K = A.matrix(); + + // Enforce Neumann conditions to right-hand side + auto g_Neumann = A.getBdrFunction(G); + A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); + + // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner + gsSparseSolver<>::CGDiagonal solver; + + // Time step + real_t dt = 0.01; + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + gsConstantFunction<> uwall_fun(u_wall,2); + auto uwall = A.getCoeff(uwall_fun, G); + // RHS of the projection + A.assemble( u * uwall * meas(G) ); + solver.compute(M); + solVector = solver.solve(A.rhs()); + + // Initialize the RHS for assembly + A.initVector(); + u.setup(bcInfo, dirichlet::l2Projection, 0); + + // Assemble the RHS + gsVector<> F = dt*A.rhs() + (M-dt*(1.0-theta)*K)*solVector; + gsVector<> F0 = F; + gsVector<> F_checkpoint = F; + + // Define the solution collection for Paraview + gsParaviewCollection collection("solution"); + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + real_t time = 0; + + // Plot initial solution + if (plot) + { + std::string fileName = "solution_" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + } + + // Time integration loop + while (interface.isCouplingOngoing()) + { + // read temperature from interface + if (interface.isReadDataAvailable()) + { + u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G) ); + K = A.matrix(); + + // Then assemble + auto g_Neumann = A.getBdrFunction(G); + A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); + A.assemble( u * ff * meas(G) ); + F = theta*dt*A.rhs() + (1.0-theta)*dt*F0 + (M-dt*(1.0-theta)*K)*solVector; + } + + // save checkpoint + if (interface.isActionRequired(interface.actionWriteIterationCheckpoint())) + { + F_checkpoint = F0; + timestep_checkpoint = timestep; + interface.markActionFulfilled(interface.actionWriteIterationCheckpoint()); + } + + // potentially adjust non-matching timestep sizes + dt = std::min(dt,precice_dt); + + // solve gismo timestep + gsInfo << "Solving timestep " << timestep*dt << "..."; + solVector = solver.compute(M + dt*theta*K).solve(F); + gsInfo<<"Finished\n"; + // write heat fluxes to interface + if (interface.isWriteDataRequired(dt)) + { + gsMatrix<> result(1,uv.cols()), tmp; + for (index_t k=0; k!=uv.cols(); k++) + { + // gsDebugVar(ev.eval(nv(G),uv.col(k))); + tmp = ev.eval(k_temp * igrad(u_sol,G),uv.col(k)); + // Only exchange y component + result(0,k) = -tmp.at(1); + } + interface.writeBlockScalarData(meshID,fluxID,xy,result); + } + + // do the coupling + precice_dt = interface.advance(dt); + + // advance variables + timestep += 1; + F0 = F; + + if (interface.isActionRequired(interface.actionReadIterationCheckpoint())) + { + F0 = F_checkpoint; + timestep = timestep_checkpoint; + interface.markActionFulfilled(interface.actionReadIterationCheckpoint()); + } + else + { + time += dt; + if (timestep % plotmod==0 && plot) + { + std::string fileName = "solution_" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + } + } + + + } + + if (plot) + { + collection.save(); + } + + + return EXIT_SUCCESS; +} diff --git a/examples/partitioned-heat-conduction-direct.cpp b/examples/partitioned-heat-conduction-direct.cpp new file mode 100644 index 0000000..776ae84 --- /dev/null +++ b/examples/partitioned-heat-conduction-direct.cpp @@ -0,0 +1,407 @@ +/** @file heat-equation-coupling.cpp + + @brief Heat equation participant for a double coupled heat equation + + 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): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) +*/ + +#include +#include +#include + +using namespace gismo; + +int main(int argc, char *argv[]) +{ + //! [Parse command line] + bool plot = false; + index_t plotmod = 1; + index_t numRefine = 5; + index_t numElevate = 0; + short_t side = 0; + 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","plotmod", "Modulo for plotting, i.e. if plotmod==1, plots every timestep", plotmod); + cmd.addInt("s","side", "Patchside of interface", side); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + //! [Read input file] + + gsMultiPatch<> patches; + if (side==0) //left + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,0.0,1.0,1.0)); + else if (side==1) //right + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(1.0,0.0,2.0,1.0)); + else + GISMO_ERROR("Side unknown"); + + real_t alpha = 3; + real_t beta = 1.3; + real_t time = 0; + real_t k_temp = 1; + + // Set external heat-flux to zero + gsConstantFunction<> f(beta-2-2*alpha,2); + gsFunctionExpr<> u_ex("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); + + gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) + + 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"; + +// ---------------------------------------------------------------------------------------------- + typedef gsExprAssembler<>::geometryMap geometryMap; + typedef gsExprAssembler<>::space space; + typedef gsExprAssembler<>::solution solution; + + gsExprAssembler<> A(1,1); + + gsInfo<<"Active options:\n"<< A.options() <<"\n"; + + A.setIntegrationElements(bases); + + gsExprEvaluator<> ev(A); + +// ---------------------------------------------------------------------------------------------- + boxSide couplingSide; + if (side==0) //left + couplingSide = boxSide(2); + else if (side==1) //right + couplingSide = boxSide(1); + else + GISMO_ERROR("Side unknown"); + + patchSide couplingInterface(0,couplingSide); + typename gsBasis::domainIter domIt = bases.basis(couplingInterface.patch).makeDomainIterator(couplingInterface.side()); + index_t rows = patches.targetDim(); + gsMatrix<> nodes; + // Start iteration over elements + gsVector<> tmp; + index_t k=0; + + gsOptionList quadOptions = A.options(); + // NEED A DIFFERENT RULE FOR dirichlet::interpolation --> see: gsGaussRule bdQuRule(basis, 1.0, 1, iter->side().direction()); + /* + quadOptions.addInt("quRule","Quadrature rule [1:GaussLegendre,2:GaussLobatto]",1); + quadOptions.addReal("quA","Number of quadrature points: quA*deg + quB",1.0); + quadOptions.addInt("quB","Number of quadrature points: quA*deg + quB",1); + */ + + index_t quadSize = 0; + typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT + for (; domIt->good(); domIt->next(), k++ ) + { + QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); + quadSize+=QuRule->numNodes(); + } + gsMatrix<> uv(rows,quadSize); + gsMatrix<> xy(rows,quadSize); + + index_t offset = 0; + + for (domIt->reset(); domIt->good(); domIt->next(), k++ ) + { + QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); + // Map the Quadrature rule to the element + QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(), + nodes, tmp); + uv.block(0,offset,rows,QuRule->numNodes()) = nodes; + + gsMatrix<> tmp2; + patches.patch(couplingInterface.patch).eval_into(nodes,tmp2); + xy.block(0,offset,rows,QuRule->numNodes()) = patches.patch(couplingInterface.patch).eval(nodes); + offset += QuRule->numNodes(); + } + + GISMO_ASSERT(side==0 || side==1,"Side must be west or east"); + std::string readName = (side == 0) ? "Dirichlet" : "Neumann"; + std::string writeName = (side == 0) ? "Neumann" : "Dirichlet"; + + gsPreCICE interface(membername, precice_config); + std::string readMeshName = readName + "-Mesh"; + std::string writeMeshName = writeName + "-Mesh"; + interface.addMesh(readMeshName,xy); + + interface.requiresInitialData(); + + real_t precice_dt = interface.initialize(); + + std::string tempName = "Temperature"; + std::string fluxName = "Heat-Flux"; + + gsMatrix<> bbox(2,2); + bbox<<0.9,1.1,-0.1,1.1; + interface.setMeshAccessRegion(writeMeshName,bbox); + + std::vector writeIDs; + gsMatrix<> writePoints; + interface.getMeshVertexIDsAndCoordinates(writeMeshName,writeIDs,writePoints); + +/* + TO DO: + * Write the evaluation on the other coordinates on the other IDs +*/ +// ---------------------------------------------------------------------------------------------- + + gsBoundaryConditions<> bcInfo; + gsPreCICEFunction g_CD(&interface,readMeshName,(side==0 ? tempName : fluxName),patches,1); + gsPreCICEFunction g_CN(&interface,readMeshName,(side==0 ? tempName : fluxName),patches,1); + gsFunction<> * g_C = &u_ex; + + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &u_ex, 0, false, 0); + bcInfo.addCondition(0, boundary::north, condition_type::dirichlet, &u_ex, 0, false, 0); + if (side==0) + { + bcInfo.addCondition(0, couplingSide,condition_type::dirichlet , g_C, 0, false, 0); + bcInfo.addCondition(0, couplingSide.opposite(),condition_type::dirichlet, &u_ex, 0, false, 0); + } + else + { + bcInfo.addCondition(0, couplingSide,condition_type::neumann , &g_CN); + bcInfo.addCondition(0, couplingSide.opposite(),condition_type::dirichlet, &u_ex, 0, false, 0); + } + + bcInfo.setGeoMap(patches); + // gsDebugVar(bcInfo); + +// ---------------------------------------------------------------------------------------------- + + // Generate system matrix and load vector + // gsInfo << "Assembling mass and stiffness...\n"; + + // Set the geometry map + geometryMap G = A.getMap(patches); + + // Set the discretization space + space u = A.getSpace(bases); + + // Set the source term + auto ff = A.getCoeff(f, G); + + // Set the solution + gsMatrix<> solVector, solVector_ini; + solution u_sol = A.getSolution(u, solVector); + solution u_ini = A.getSolution(u, solVector); + + u.setup(bcInfo, dirichlet::homogeneous, 0); + A.initSystem(); + gsDebugVar(A.numDofs()); + A.assemble( u * u.tr() * meas(G)); + gsSparseMatrix<> M = A.matrix(); + + + // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner + gsSparseSolver<>::CGDiagonal solver; + + real_t dt = 0.01; + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + auto uex = A.getCoeff(u_ex, G); + // RHS of the projection + u.setup(bcInfo, dirichlet::l2Projection, 0); + A.initSystem(); + A.assemble( u * u.tr() * meas(G), u * uex * meas(G) ); + solver.compute(A.matrix()); + solVector_ini = solVector = solver.solve(A.rhs()); + + gsMatrix<> result(1,uv.cols()), tmp2; + // Write initial data + if (interface.requiresWritingCheckpoint()) + { + for (index_t k=0; k!=uv.cols(); k++) + { + if (side==0) + { + gsWarn<<"Write the flux here!!!\n"; + tmp2 = ev.eval( - jac(u_sol) * nv(G).normalized(),writePoints.col(k)); + result(0,k) = tmp2.at(0); + } + else + { + tmp2 = ev.eval(u_sol,writePoints.col(k)); + result(0,k) = tmp2.at(0); + } + } + // gsDebugVar(result); + interface.writeData(writeMeshName,writeName=="Dirichlet" ? fluxName : tempName,writeIDs,result); + } + // interface.initialize_data(); + + // interface.readBlockScalarData(readMeshName,side==0 ? tempName : fluxName,xy,result); + // gsDebugVar(result); + + // Initialize the RHS for assembly + if (side==0) + g_C = &g_CD; + gsDebugVar(bcInfo); + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * uex * meas(G) ); + gsSparseMatrix<> K = A.matrix(); + + // Assemble the RHS + gsVector<> F = dt*A.rhs() + M*solVector; + gsVector<> F0 = F; + gsVector<> F_checkpoint = F; + gsMatrix<> solVector_checkpoint = solVector; + + gsParaviewCollection collection("solution"); + gsParaviewCollection exact_collection("exact_solution"); + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + real_t t_checkpoint = 0; + if (plot) + { + std::string fileName = "solution_" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + + // fileName = "exact_solution_" + util::to_string(timestep); + // ev.writeParaview( uex , G, fileName); + // for (size_t p=0; p!=patches.nPatches(); p++) + // { + // fileName = "exact_solution_" + util::to_string(timestep) + std::to_string(p); + // exact_collection.addTimestep(fileName,time,".vts"); + // } + + ev.writeParaview( u_sol , G, "initial_solution"); + + } + + // time += precice_dt; + while (interface.isCouplingOngoing()) + { + // read temperature from interface + // if (interface.isReadDataAvailable()) + // { + u_ex = gsFunctionExpr<>("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); + u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: + + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * ff * meas(G) ); + K = A.matrix(); + auto g_Neumann = A.getBdrFunction(G); + A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); + F = dt*A.rhs() + M*solVector; + + // gsMatrix<> result; + // interface.readBlockScalarData(readMeshName,side==0 ? tempName : fluxName,xy,result); + // gsDebugVar(result); + // } + + // save checkpoint + if (interface.requiresWritingCheckpoint()) + { + gsDebugVar("Write checkpoint"); + F_checkpoint = F0; + t_checkpoint = time; + timestep_checkpoint = timestep; + solVector_checkpoint = solVector; + } + + // potentially adjust non-matching timestep sizes + dt = std::min(dt,precice_dt); + + // solve gismo timestep + gsInfo << "Solving timestep " << timestep*dt << "..."; + solVector = solver.compute(M + dt*K).solve(F); + gsInfo<<"Finished\n"; + // write heat fluxes to interface + gsMatrix<> result(1,uv.cols()), tmp; + for (index_t k=0; k!=uv.cols(); k++) + { + // gsDebugVar(ev.eval(nv(G),uv.col(k))); + // tmp = ev.eval(k_temp * igrad(u_sol,G),uv.col(k)); + // Only exchange y component + // result(0,k) = -tmp.at(1); + // + // + if (side==0) + { + gsWarn<<"Write the flux here!!!\n"; + tmp = ev.eval(jac(u_sol) * nv(G).normalized(),uv.col(k)); + result(0,k) = tmp.at(0); + } + else + { + tmp = ev.eval(u_sol,uv.col(k)); + result(0,k) = tmp.at(0); + } + } + // TODO + interface.writeData(readMeshName,side==0 ? fluxName : tempName,xy,result); + + // do the coupling + precice_dt = interface.advance(dt); + + // advance variables + time += dt; + timestep += 1; + F0 = F; + + if (interface.requiresReadingCheckpoint()) + { + gsDebugVar("Read checkpoint"); + F0 = F_checkpoint; + time = t_checkpoint; + timestep = timestep_checkpoint; + solVector = solVector_checkpoint; + } + else + { + if (timestep % plotmod==0 && plot) + { + std::string fileName = "solution_" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + + // fileName = "exact_solution_" + util::to_string(timestep); + // ev.writeParaview( uex , G, fileName); + // for (size_t p=0; p!=patches.nPatches(); p++) + // { + // fileName = "exact_solution_" + util::to_string(timestep) + std::to_string(p); + // exact_collection.addTimestep(fileName,time,".vts"); + // } + } + } + } + + if (plot) + { + collection.save(); + exact_collection.save(); + } + + + return EXIT_SUCCESS; +} diff --git a/examples/partitioned-heat-conduction/README.md b/examples/partitioned-heat-conduction/README.md index 1009251..56378ac 100644 --- a/examples/partitioned-heat-conduction/README.md +++ b/examples/partitioned-heat-conduction/README.md @@ -23,9 +23,7 @@ preCICE configuration (image generated using the [precice-config-visualizer](htt In the G+Smo build folder, build the tutorial file. ``` -make partitioned-heat-conduction -j <#threads> -make partitioned-heat-conduction-IGA-dirichlet -j <#threads> -make partitioned-heat-conduction-IGA-neumann -j <#threads> +make -j <#threads> ``` Go to the gismo-executable folder and link the compiled executable to the gismo_executable. diff --git a/gsLookupFunction.h b/gsLookupFunction.h index 98e3593..09245a5 100644 --- a/gsLookupFunction.h +++ b/gsLookupFunction.h @@ -93,17 +93,26 @@ class gsLookupFunction : public gsFunction virtual short_t targetDim() const { return m_data.rows(); } - /// See \a gsFunction + /** \brief Evaluate the function at points \a u into \a result. + * This function evaluates a target function at a given set of input points and + * stores the results in the provided output matrix. The mapping between input + * points and corresponding results is defined by the internal mapping table \c m_map. + * The function ensures that all input points are valid and registered in the table + * before performing the evaluation. + * + * \param[in] u A \c gsMatrix , where each column represents a point in the parameter domain. + * \param[out] result A \c gsMatrix , where each column will contain the result of + * evaluating the function at the corresponding input point. + * + */ virtual void eval_into(const gsMatrix& u, gsMatrix& result) const { index_t col; - gsDebugVar(u); result.resize(this->targetDim(),u.cols()); result.setZero(); - for (index_t k = 0; k!= u.cols(); k++) { - // gsDebugVar(u.col(k)); + GISMO_ASSERT(m_map.find(u.col(k))!=m_map.end(),"Coordinate " + std::to_string(k) + " not registered in the table"); col = m_map.at(u.col(k)); result.col(k) = m_data.col(col); @@ -141,14 +150,6 @@ class gsLookupFunction : public gsFunction GISMO_ASSERT(m_points.cols()==m_data.cols(),"Points and data must have the same number of columns"); for (index_t k = 0; k != m_points.cols(); k++) m_map.insert({m_points.col(k),k}); // m_map.at(vector) returns the column index of vector - gsInfo << "Content of m_map:\n"; - for (const auto& [key, value] : m_map) { - gsInfo << "Key: "; - for (int i = 0; i < key.size(); ++i) { - gsInfo << key[i] << " "; - } - gsInfo << "=> Value: " << value << "\n"; - } } /// See \a gsFunction diff --git a/gsPreCICEUtils.h b/gsPreCICEUtils.h index 9f9c2b5..65dfd58 100644 --- a/gsPreCICEUtils.h +++ b/gsPreCICEUtils.h @@ -56,7 +56,7 @@ gsVector knotsToVector(const gsBasis & basis) GISMO_ERROR("Basis type not understood"); break; } - gsDebugVar(tensorKnots[0]); + std::vector sizes(DIM); sizes[0] = tensorKnots[0].size(); @@ -172,17 +172,18 @@ gsMatrix knotVectorUnpack(const gsMatrix & knots, index_t numBoundaries) // } /** - * @brief Pack the knot and ratio matrices of a gsMultiPatch object into a single matrix. + * @brief Pack the knot and control points matrices of a gsMultiPatch object into a single matrix. * - * @param[in] mp The gsMultiPatch object + * @param[in] mp The gsMultiPatch object. MultiPatch geometry object. * - * @tparam T { description } + * @tparam T Real type. * * @return A tuple containing the packed knot matrix, the packed ratio matrix, the number of columns in each patch's knot matrix, and the number of columns in each patch's ratio matrix. */ template -std::tuple, gsMatrix, std::vector, std::vector> packMultiPatch(const gsMultiPatch &mp) { +std::tuple, gsMatrix, std::vector, std::vector> packMultiPatch(const gsMultiPatch &mp) +{ std::vector> knotMatrices; knotMatrices.reserve(mp.nPatches()); std::vector> coefMatrices; @@ -242,7 +243,7 @@ std::tuple, gsMatrix, std::vector, std::vector> * @param[in] knotMatrices The knot matrices * @param[in] coefMatrices The ratio matrices * - * @tparam T { description } + * @tparam T Reak type. * * @return A gsMultiPatch object. */ @@ -339,10 +340,9 @@ typename gsBasis::Ptr knotMatrixToBasis(const gsMatrix & knots) std::copy_if(knots.row(d).begin(),knots.row(d).end(), std::back_inserter(tmp), [](T a){return !math::isnan(a);}); - // KVs[d] = gsKnotVector(tmp);/////// MEMLEAK - KVs[d] = std::move(gsKnotVector(tmp)); - gsDebugVar(d); - gsDebugVar(KVs[d]); + KVs[d] = give(gsKnotVector(tmp)); + + gsDebug<<"(gsPreCICEUtils.h: There is a memory leak in the line above)\n"; } From 834373e16945f8361911f04fae81538323964e13 Mon Sep 17 00:00:00 2001 From: Crazy-Rich-Meghan Date: Mon, 9 Dec 2024 20:49:37 +0100 Subject: [PATCH 12/14] Clean up --- #CMakeLists.txt# | 86 ---- examples/flow-over-heated-plate.cpp | 316 -------------- .../partitioned-heat-conduction-direct.cpp | 407 ------------------ .../.DS_Store | Bin 8196 -> 0 bytes .../dirichlet-gismo/.DS_Store | Bin 6148 -> 0 bytes .../gismo-executable/.DS_Store | Bin 6148 -> 0 bytes .../neumann-gismo/.DS_Store | Bin 6148 -> 0 bytes .../.DS_Store | Bin 6148 -> 0 bytes .../gismo-executable/.DS_Store | Bin 6148 -> 0 bytes 9 files changed, 809 deletions(-) delete mode 100644 #CMakeLists.txt# delete mode 100644 examples/flow-over-heated-plate.cpp delete mode 100644 examples/partitioned-heat-conduction-direct.cpp delete mode 100644 examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/.DS_Store delete mode 100644 examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/.DS_Store delete mode 100644 examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/gismo-executable/.DS_Store delete mode 100644 examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/neumann-gismo/.DS_Store delete mode 100644 examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/.DS_Store delete mode 100644 examples/partitioned-heat-conduction/partitioned-heat-conduction-vertex-vertex/gismo-executable/.DS_Store diff --git a/#CMakeLists.txt# b/#CMakeLists.txt# deleted file mode 100644 index 0280e55..0000000 --- a/#CMakeLists.txt# +++ /dev/null @@ -1,86 +0,0 @@ -###################################################################### -## CMakeLists.txt --- gsPreCICE -## This file is part of the G+Smo library. -## -###################################################################### - -if(POLICY CMP0076) -cmake_policy(SET CMP0076 NEW) -endif() - - -## gsPreCICE module -project(gsPreCICE) - -# Apply G+Smo config -include(gsConfig) - - -## Collect files -aux_header_directory (${CMAKE_CURRENT_SOURCE_DIR}/src ${PROJECT_NAME}_H ) - - - - -if( (NOT GISMO_BUILD_LIB) ) - aux_instance_directory (${CMAKE_CURRENT_SOURCE_DIR} ${PROJECT_NAME}_INS) - if(${PROJECT_NAME}_INS) - LIST( REMOVE_ITEM ${PROJECT_NAME}_CPP ${${PROJECT_NAME}_INS}) - endif() -endif() - -# Add object library -add_library(${PROJECT_NAME} OBJECT - ${${PROJECT_NAME}_H} - ${${PROJECT_NAME}_HPP} - ${${PROJECT_NAME}_CPP} ) - -set_target_properties(${PROJECT_NAME} PROPERTIES - COMPILE_DEFINITIONS gismo_EXPORTS - POSITION_INDEPENDENT_CODE ON - LINKER_LANGUAGE CXX - FOLDER "G+Smo modules" ) - -set(gismo_MODULES ${gismo_MODULES} $ - CACHE INTERNAL "G+Smo modules" ) - -#Symlink include dir (in case your headers are in /src) -execute_process( COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_CURRENT_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/${PROJECT_NAME}) - - -# Search for preCICE -find_package(precice REQUIRED CONFIG) - -target_link_libraries(${PROJECT_NAME} PRIVATE precice::precice) - -#set_property(TARGET ${PROJECT_NAME} APPEND PROPERTY INTERFACE_LINK_LIBRARIES precice::precice) - -target_include_directories(${PROJECT_NAME} INTERFACE -$ -) - -add_dependencies(${PROJECT_NAME} precice) - - -install(DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/${PROJECT_NAME}" - DESTINATION include/gismo - FILES_MATCHING PATTERN "*.h" ) - -# add filedata folder -add_definitions(-DELAST_DATA_DIR="${CMAKE_CURRENT_SOURCE_DIR}/filedata/") - -# add example files -aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR}/examples FILES) -foreach(file ${FILES}) - add_gismo_executable(${file}) - get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") - set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") - # install the example executables (optionally) - install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) -endforeach(file ${FILES}) - -# add unittests -aux_gs_cpp_directory(${PROJECT_SOURCE_DIR}/unittests unittests_SRCS) -set(gismo_UNITTESTS ${gismo_UNITTESTS} ${unittests_SRCS} - CACHE INTERNAL "gismo list of unittests") diff --git a/examples/flow-over-heated-plate.cpp b/examples/flow-over-heated-plate.cpp deleted file mode 100644 index d4cca4e..0000000 --- a/examples/flow-over-heated-plate.cpp +++ /dev/null @@ -1,316 +0,0 @@ -/** @file flow-over-heated-plate.cpp - - @brief Heat equation participant for the PreCICE example "flow over heated plate" - - 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): H.M. Verhelst -*/ - -#include -#include -#include - -using namespace gismo; - -int main(int argc, char *argv[]) -{ - //! [Parse command line] - bool plot = false; - index_t plotmod = 1; - index_t numRefine = 5; - index_t numElevate = 0; - std::string precice_config("../precice_config.xml"); - - 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); - try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } - - //! [Read input file] - - // Generate domain - gsMultiPatch<> patches; - patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,-0.25,1.0,0.0)); - - // Set external heat-flux to zero - gsFunctionExpr<> f("0",2); - - // 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"; - - // Set heat conduction coefficient - real_t k_temp = 100; - // Set bottom wall temp - real_t u_wall = 310; - -// ---------------------------------------------------------------------------------------------- - typedef gsExprAssembler<>::geometryMap geometryMap; - typedef gsExprAssembler<>::space space; - typedef gsExprAssembler<>::solution solution; - - gsExprAssembler<> A(1,1); - - gsInfo<<"Active options:\n"<< A.options() <<"\n"; - - A.setIntegrationElements(bases); - - gsExprEvaluator<> ev(A); - -// ---------------------------------------------------------------------------------------------- - // Set the interface for the precice coupling - patchSide couplingInterface(0,boundary::north); - - // Get a domain iterator on the coupling interface - typename gsBasis::domainIter domIt = bases.basis(couplingInterface.patch).makeDomainIterator(couplingInterface.side()); - - // Set the dimension of the points - gsMatrix<> nodes; - // Start iteration over elements - gsVector<> tmp; - - gsOptionList quadOptions = A.options(); - - // First obtain the size of all quadrature points - index_t quadSize = 0; - typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT - for (; domIt->good(); domIt->next() ) - { - QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); - quadSize+=QuRule->numNodes(); - } - - // Initialize parametric coordinates - gsMatrix<> uv(patches.domainDim(),quadSize); - // Initialize physical coordinates - gsMatrix<> xy(patches.targetDim(),quadSize); - - // Grab all quadrature points - index_t offset = 0; - for (domIt->reset(); domIt->good(); domIt->next()) - { - QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); - // Map the Quadrature rule to the element - QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(), - nodes, tmp); - uv.block(0,offset,patches.domainDim(),QuRule->numNodes()) = nodes; - - gsMatrix<> tmp2; - patches.patch(couplingInterface.patch).eval_into(nodes,tmp2); - xy.block(0,offset,patches.targetDim(),QuRule->numNodes()) = patches.patch(couplingInterface.patch).eval(nodes); - offset += QuRule->numNodes(); - } - - // Define precice interface - gsPreCICE interface("Solid", precice_config); - interface.addMesh("Solid-Mesh",xy); - real_t precice_dt = interface.initialize(); - - index_t meshID = interface.getMeshID("Solid-Mesh"); - index_t tempID = interface.getDataID("Temperature",meshID); - index_t fluxID = interface.getDataID("Heat-Flux",meshID); - -// ---------------------------------------------------------------------------------------------- - - // Define boundary conditions - gsBoundaryConditions<> bcInfo; - // Dirichlet side - gsConstantFunction<> g_D(u_wall,2); - // Homogeneous Neumann - gsConstantFunction<> g_N(0.0,2); - // Coupling side - gsPreCICEFunction g_C(&interface,meshID,tempID,patches); - // Add all BCs - // Isolated Neumann sides - bcInfo.addCondition(0, boundary::east, condition_type::neumann , &g_N, 0, false, 0); - bcInfo.addCondition(0, boundary::west, condition_type::neumann , &g_N, 0, false, 0); - // Bottom side (prescribed temp) - bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &g_D, 0, false, 0); - // Coupling interface - bcInfo.addCondition(0, boundary::north, condition_type::dirichlet, &g_C, 0, false, 0); - // Assign geometry map - bcInfo.setGeoMap(patches); - -// ---------------------------------------------------------------------------------------------- - - // Time integration coefficient (0.0 = explicit, 1.0 = implicit) - real_t theta = 1.0; - - // Set the geometry map - geometryMap G = A.getMap(patches); - - // Set the discretization space - space u = A.getSpace(bases); - u.setup(bcInfo, dirichlet::homogeneous, 0); - - // Set the source term - auto ff = A.getCoeff(f, G); - - // Set the solution - gsMatrix<> solVector; - solution u_sol = A.getSolution(u, solVector); - - // Assemble mass matrix - A.initSystem(); - A.assemble( u * u.tr() * meas(G)); - gsSparseMatrix<> M = A.matrix(); - - // Assemble stiffness matrix (NOTE: also adds the dirichlet BCs inside the matrixs) - A.initSystem(); - A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G) ); - gsSparseMatrix<> K = A.matrix(); - - // Enforce Neumann conditions to right-hand side - auto g_Neumann = A.getBdrFunction(G); - A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); - - // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner - gsSparseSolver<>::CGDiagonal solver; - - // Time step - real_t dt = 0.01; - - // Project u_wall as initial condition (violates Dirichlet side on precice interface) - gsConstantFunction<> uwall_fun(u_wall,2); - auto uwall = A.getCoeff(uwall_fun, G); - // RHS of the projection - A.assemble( u * uwall * meas(G) ); - solver.compute(M); - solVector = solver.solve(A.rhs()); - - // Initialize the RHS for assembly - A.initVector(); - u.setup(bcInfo, dirichlet::l2Projection, 0); - - // Assemble the RHS - gsVector<> F = dt*A.rhs() + (M-dt*(1.0-theta)*K)*solVector; - gsVector<> F0 = F; - gsVector<> F_checkpoint = F; - - // Define the solution collection for Paraview - gsParaviewCollection collection("solution"); - - index_t timestep = 0; - index_t timestep_checkpoint = 0; - real_t time = 0; - - // Plot initial solution - if (plot) - { - std::string fileName = "solution_" + util::to_string(timestep); - ev.options().setSwitch("plot.elements", true); - ev.options().setInt("plot.npts", 1000); - ev.writeParaview( u_sol , G, fileName); - for (size_t p=0; p!=patches.nPatches(); p++) - { - fileName = "solution_" + util::to_string(timestep) + std::to_string(p); - collection.addTimestep(fileName,time,".vts"); - } - } - - // Time integration loop - while (interface.isCouplingOngoing()) - { - // read temperature from interface - if (interface.isReadDataAvailable()) - { - u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: - A.initSystem(); - A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G) ); - K = A.matrix(); - - // Then assemble - auto g_Neumann = A.getBdrFunction(G); - A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); - A.assemble( u * ff * meas(G) ); - F = theta*dt*A.rhs() + (1.0-theta)*dt*F0 + (M-dt*(1.0-theta)*K)*solVector; - } - - // save checkpoint - if (interface.isActionRequired(interface.actionWriteIterationCheckpoint())) - { - F_checkpoint = F0; - timestep_checkpoint = timestep; - interface.markActionFulfilled(interface.actionWriteIterationCheckpoint()); - } - - // potentially adjust non-matching timestep sizes - dt = std::min(dt,precice_dt); - - // solve gismo timestep - gsInfo << "Solving timestep " << timestep*dt << "..."; - solVector = solver.compute(M + dt*theta*K).solve(F); - gsInfo<<"Finished\n"; - // write heat fluxes to interface - if (interface.isWriteDataRequired(dt)) - { - gsMatrix<> result(1,uv.cols()), tmp; - for (index_t k=0; k!=uv.cols(); k++) - { - // gsDebugVar(ev.eval(nv(G),uv.col(k))); - tmp = ev.eval(k_temp * igrad(u_sol,G),uv.col(k)); - // Only exchange y component - result(0,k) = -tmp.at(1); - } - interface.writeBlockScalarData(meshID,fluxID,xy,result); - } - - // do the coupling - precice_dt = interface.advance(dt); - - // advance variables - timestep += 1; - F0 = F; - - if (interface.isActionRequired(interface.actionReadIterationCheckpoint())) - { - F0 = F_checkpoint; - timestep = timestep_checkpoint; - interface.markActionFulfilled(interface.actionReadIterationCheckpoint()); - } - else - { - time += dt; - if (timestep % plotmod==0 && plot) - { - std::string fileName = "solution_" + util::to_string(timestep); - ev.options().setSwitch("plot.elements", true); - ev.options().setInt("plot.npts", 1000); - ev.writeParaview( u_sol , G, fileName); - for (size_t p=0; p!=patches.nPatches(); p++) - { - fileName = "solution_" + util::to_string(timestep) + std::to_string(p); - collection.addTimestep(fileName,time,".vts"); - } - } - } - - - } - - if (plot) - { - collection.save(); - } - - - return EXIT_SUCCESS; -} diff --git a/examples/partitioned-heat-conduction-direct.cpp b/examples/partitioned-heat-conduction-direct.cpp deleted file mode 100644 index 776ae84..0000000 --- a/examples/partitioned-heat-conduction-direct.cpp +++ /dev/null @@ -1,407 +0,0 @@ -/** @file heat-equation-coupling.cpp - - @brief Heat equation participant for a double coupled heat equation - - 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): H.M. Verhelst (University of Pavia), J.Li (TU Delft, 2023-...) -*/ - -#include -#include -#include - -using namespace gismo; - -int main(int argc, char *argv[]) -{ - //! [Parse command line] - bool plot = false; - index_t plotmod = 1; - index_t numRefine = 5; - index_t numElevate = 0; - short_t side = 0; - 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","plotmod", "Modulo for plotting, i.e. if plotmod==1, plots every timestep", plotmod); - cmd.addInt("s","side", "Patchside of interface", side); - try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } - - //! [Read input file] - - gsMultiPatch<> patches; - if (side==0) //left - patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,0.0,1.0,1.0)); - else if (side==1) //right - patches.addPatch(gsNurbsCreator<>::BSplineRectangle(1.0,0.0,2.0,1.0)); - else - GISMO_ERROR("Side unknown"); - - real_t alpha = 3; - real_t beta = 1.3; - real_t time = 0; - real_t k_temp = 1; - - // Set external heat-flux to zero - gsConstantFunction<> f(beta-2-2*alpha,2); - gsFunctionExpr<> u_ex("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); - - gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) - - 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"; - -// ---------------------------------------------------------------------------------------------- - typedef gsExprAssembler<>::geometryMap geometryMap; - typedef gsExprAssembler<>::space space; - typedef gsExprAssembler<>::solution solution; - - gsExprAssembler<> A(1,1); - - gsInfo<<"Active options:\n"<< A.options() <<"\n"; - - A.setIntegrationElements(bases); - - gsExprEvaluator<> ev(A); - -// ---------------------------------------------------------------------------------------------- - boxSide couplingSide; - if (side==0) //left - couplingSide = boxSide(2); - else if (side==1) //right - couplingSide = boxSide(1); - else - GISMO_ERROR("Side unknown"); - - patchSide couplingInterface(0,couplingSide); - typename gsBasis::domainIter domIt = bases.basis(couplingInterface.patch).makeDomainIterator(couplingInterface.side()); - index_t rows = patches.targetDim(); - gsMatrix<> nodes; - // Start iteration over elements - gsVector<> tmp; - index_t k=0; - - gsOptionList quadOptions = A.options(); - // NEED A DIFFERENT RULE FOR dirichlet::interpolation --> see: gsGaussRule bdQuRule(basis, 1.0, 1, iter->side().direction()); - /* - quadOptions.addInt("quRule","Quadrature rule [1:GaussLegendre,2:GaussLobatto]",1); - quadOptions.addReal("quA","Number of quadrature points: quA*deg + quB",1.0); - quadOptions.addInt("quB","Number of quadrature points: quA*deg + quB",1); - */ - - index_t quadSize = 0; - typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT - for (; domIt->good(); domIt->next(), k++ ) - { - QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); - quadSize+=QuRule->numNodes(); - } - gsMatrix<> uv(rows,quadSize); - gsMatrix<> xy(rows,quadSize); - - index_t offset = 0; - - for (domIt->reset(); domIt->good(); domIt->next(), k++ ) - { - QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); - // Map the Quadrature rule to the element - QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(), - nodes, tmp); - uv.block(0,offset,rows,QuRule->numNodes()) = nodes; - - gsMatrix<> tmp2; - patches.patch(couplingInterface.patch).eval_into(nodes,tmp2); - xy.block(0,offset,rows,QuRule->numNodes()) = patches.patch(couplingInterface.patch).eval(nodes); - offset += QuRule->numNodes(); - } - - GISMO_ASSERT(side==0 || side==1,"Side must be west or east"); - std::string readName = (side == 0) ? "Dirichlet" : "Neumann"; - std::string writeName = (side == 0) ? "Neumann" : "Dirichlet"; - - gsPreCICE interface(membername, precice_config); - std::string readMeshName = readName + "-Mesh"; - std::string writeMeshName = writeName + "-Mesh"; - interface.addMesh(readMeshName,xy); - - interface.requiresInitialData(); - - real_t precice_dt = interface.initialize(); - - std::string tempName = "Temperature"; - std::string fluxName = "Heat-Flux"; - - gsMatrix<> bbox(2,2); - bbox<<0.9,1.1,-0.1,1.1; - interface.setMeshAccessRegion(writeMeshName,bbox); - - std::vector writeIDs; - gsMatrix<> writePoints; - interface.getMeshVertexIDsAndCoordinates(writeMeshName,writeIDs,writePoints); - -/* - TO DO: - * Write the evaluation on the other coordinates on the other IDs -*/ -// ---------------------------------------------------------------------------------------------- - - gsBoundaryConditions<> bcInfo; - gsPreCICEFunction g_CD(&interface,readMeshName,(side==0 ? tempName : fluxName),patches,1); - gsPreCICEFunction g_CN(&interface,readMeshName,(side==0 ? tempName : fluxName),patches,1); - gsFunction<> * g_C = &u_ex; - - bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &u_ex, 0, false, 0); - bcInfo.addCondition(0, boundary::north, condition_type::dirichlet, &u_ex, 0, false, 0); - if (side==0) - { - bcInfo.addCondition(0, couplingSide,condition_type::dirichlet , g_C, 0, false, 0); - bcInfo.addCondition(0, couplingSide.opposite(),condition_type::dirichlet, &u_ex, 0, false, 0); - } - else - { - bcInfo.addCondition(0, couplingSide,condition_type::neumann , &g_CN); - bcInfo.addCondition(0, couplingSide.opposite(),condition_type::dirichlet, &u_ex, 0, false, 0); - } - - bcInfo.setGeoMap(patches); - // gsDebugVar(bcInfo); - -// ---------------------------------------------------------------------------------------------- - - // Generate system matrix and load vector - // gsInfo << "Assembling mass and stiffness...\n"; - - // Set the geometry map - geometryMap G = A.getMap(patches); - - // Set the discretization space - space u = A.getSpace(bases); - - // Set the source term - auto ff = A.getCoeff(f, G); - - // Set the solution - gsMatrix<> solVector, solVector_ini; - solution u_sol = A.getSolution(u, solVector); - solution u_ini = A.getSolution(u, solVector); - - u.setup(bcInfo, dirichlet::homogeneous, 0); - A.initSystem(); - gsDebugVar(A.numDofs()); - A.assemble( u * u.tr() * meas(G)); - gsSparseMatrix<> M = A.matrix(); - - - // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner - gsSparseSolver<>::CGDiagonal solver; - - real_t dt = 0.01; - - // Project u_wall as initial condition (violates Dirichlet side on precice interface) - auto uex = A.getCoeff(u_ex, G); - // RHS of the projection - u.setup(bcInfo, dirichlet::l2Projection, 0); - A.initSystem(); - A.assemble( u * u.tr() * meas(G), u * uex * meas(G) ); - solver.compute(A.matrix()); - solVector_ini = solVector = solver.solve(A.rhs()); - - gsMatrix<> result(1,uv.cols()), tmp2; - // Write initial data - if (interface.requiresWritingCheckpoint()) - { - for (index_t k=0; k!=uv.cols(); k++) - { - if (side==0) - { - gsWarn<<"Write the flux here!!!\n"; - tmp2 = ev.eval( - jac(u_sol) * nv(G).normalized(),writePoints.col(k)); - result(0,k) = tmp2.at(0); - } - else - { - tmp2 = ev.eval(u_sol,writePoints.col(k)); - result(0,k) = tmp2.at(0); - } - } - // gsDebugVar(result); - interface.writeData(writeMeshName,writeName=="Dirichlet" ? fluxName : tempName,writeIDs,result); - } - // interface.initialize_data(); - - // interface.readBlockScalarData(readMeshName,side==0 ? tempName : fluxName,xy,result); - // gsDebugVar(result); - - // Initialize the RHS for assembly - if (side==0) - g_C = &g_CD; - gsDebugVar(bcInfo); - A.initSystem(); - A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * uex * meas(G) ); - gsSparseMatrix<> K = A.matrix(); - - // Assemble the RHS - gsVector<> F = dt*A.rhs() + M*solVector; - gsVector<> F0 = F; - gsVector<> F_checkpoint = F; - gsMatrix<> solVector_checkpoint = solVector; - - gsParaviewCollection collection("solution"); - gsParaviewCollection exact_collection("exact_solution"); - - index_t timestep = 0; - index_t timestep_checkpoint = 0; - real_t t_checkpoint = 0; - if (plot) - { - std::string fileName = "solution_" + util::to_string(timestep); - ev.options().setSwitch("plot.elements", true); - ev.options().setInt("plot.npts", 1000); - ev.writeParaview( u_sol , G, fileName); - for (size_t p=0; p!=patches.nPatches(); p++) - { - fileName = "solution_" + util::to_string(timestep) + std::to_string(p); - collection.addTimestep(fileName,time,".vts"); - } - - // fileName = "exact_solution_" + util::to_string(timestep); - // ev.writeParaview( uex , G, fileName); - // for (size_t p=0; p!=patches.nPatches(); p++) - // { - // fileName = "exact_solution_" + util::to_string(timestep) + std::to_string(p); - // exact_collection.addTimestep(fileName,time,".vts"); - // } - - ev.writeParaview( u_sol , G, "initial_solution"); - - } - - // time += precice_dt; - while (interface.isCouplingOngoing()) - { - // read temperature from interface - // if (interface.isReadDataAvailable()) - // { - u_ex = gsFunctionExpr<>("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); - u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: - - A.initSystem(); - A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * ff * meas(G) ); - K = A.matrix(); - auto g_Neumann = A.getBdrFunction(G); - A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); - F = dt*A.rhs() + M*solVector; - - // gsMatrix<> result; - // interface.readBlockScalarData(readMeshName,side==0 ? tempName : fluxName,xy,result); - // gsDebugVar(result); - // } - - // save checkpoint - if (interface.requiresWritingCheckpoint()) - { - gsDebugVar("Write checkpoint"); - F_checkpoint = F0; - t_checkpoint = time; - timestep_checkpoint = timestep; - solVector_checkpoint = solVector; - } - - // potentially adjust non-matching timestep sizes - dt = std::min(dt,precice_dt); - - // solve gismo timestep - gsInfo << "Solving timestep " << timestep*dt << "..."; - solVector = solver.compute(M + dt*K).solve(F); - gsInfo<<"Finished\n"; - // write heat fluxes to interface - gsMatrix<> result(1,uv.cols()), tmp; - for (index_t k=0; k!=uv.cols(); k++) - { - // gsDebugVar(ev.eval(nv(G),uv.col(k))); - // tmp = ev.eval(k_temp * igrad(u_sol,G),uv.col(k)); - // Only exchange y component - // result(0,k) = -tmp.at(1); - // - // - if (side==0) - { - gsWarn<<"Write the flux here!!!\n"; - tmp = ev.eval(jac(u_sol) * nv(G).normalized(),uv.col(k)); - result(0,k) = tmp.at(0); - } - else - { - tmp = ev.eval(u_sol,uv.col(k)); - result(0,k) = tmp.at(0); - } - } - // TODO - interface.writeData(readMeshName,side==0 ? fluxName : tempName,xy,result); - - // do the coupling - precice_dt = interface.advance(dt); - - // advance variables - time += dt; - timestep += 1; - F0 = F; - - if (interface.requiresReadingCheckpoint()) - { - gsDebugVar("Read checkpoint"); - F0 = F_checkpoint; - time = t_checkpoint; - timestep = timestep_checkpoint; - solVector = solVector_checkpoint; - } - else - { - if (timestep % plotmod==0 && plot) - { - std::string fileName = "solution_" + util::to_string(timestep); - ev.options().setSwitch("plot.elements", true); - ev.options().setInt("plot.npts", 1000); - ev.writeParaview( u_sol , G, fileName); - for (size_t p=0; p!=patches.nPatches(); p++) - { - fileName = "solution_" + util::to_string(timestep) + std::to_string(p); - collection.addTimestep(fileName,time,".vts"); - } - - // fileName = "exact_solution_" + util::to_string(timestep); - // ev.writeParaview( uex , G, fileName); - // for (size_t p=0; p!=patches.nPatches(); p++) - // { - // fileName = "exact_solution_" + util::to_string(timestep) + std::to_string(p); - // exact_collection.addTimestep(fileName,time,".vts"); - // } - } - } - } - - if (plot) - { - collection.save(); - exact_collection.save(); - } - - - return EXIT_SUCCESS; -} diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/.DS_Store b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/.DS_Store deleted file mode 100644 index fcf242967c60dc3550512971a6603daeac0efcef..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8196 zcmeHMO^?z*7=C9#w`8*slQmgAHSxMeST?d76VpY!day5cO{uoa_^Rd_f+v;juqBF_NGwr<3d_0{G-hqf%V;oe8_KCOi0-Udn-@j_qkvJsC}0#Y3T%P`II~5t&N=s0sY#6jMuGoQ0kJ>W zhyoprbA@v2Kq6fLU=gonKpuX8C4p{#oAhpM{-CtGxG+$0IL{J3~nb zgAPnyJ`dx7Pa1q22Z`jms=z8)rB?aobXu)cDhGD$@OXA$Pfre!KdRlE%}Umd+js9j zX}yX1F@LW{#`Nw$yGtYekUz{IX-7V~EtolqMBS>}=KVpcksL;!H=9~RtFU+L^Mlr> z)BJY&Jh%XB(EHWiy>*9i#N#e{Bm^we;@jgMr$g@$+cd_ljz_+@-?{7=aa23Tgs9iW zgs4{QCl{R&xdTnliCp6)6?Z9yZey$<54jjhE>yj2{r=R7${1o=xeT*4D=EV|y^3a! zzkSiGD3kGH1^paj&Faahf9ad2HUyJ33T(In1=YWO;r~B6|Nei&M`5~W6fg?>jsmRM zbej#-aQ5qYxiaBv+la3bMG$s#g_43qy5m6Ujsq8e7$R@ORB#-PbA@<=@YgQ_PSI$_ PiEhq+`B_Zzk1OyUr^#ER diff --git a/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/.DS_Store b/examples/partitioned-heat-conduction/partitioned-heat-conduction-IGA-IGA/dirichlet-gismo/.DS_Store deleted file mode 100644 index 5008ddfcf53c02e82d7eee2e57c38e5672ef89f6..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0HcWPxsBv7*Y&)5JY4KX5Z}Y z%xv~Sc6R{Cs;}>X0>F~0C~C}z4$n?qxbuQ2IY)&jG}z-6cQc9p;*fs(4Qf2&NS%KF zf&CL+(AB4^*>zpREq>?vwmG&<*Bw4)NkceC5yLMe^#*pRzloG>kSF2nNm> zxU_De_x}Tbnb{)$JS2L-Krryn7?8z!vtDsiez$(vp5C>I`bbq#zd;QO?XyP!J9>`X e=1H4RYSXWD91Jy!j^}VH1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0 Date: Mon, 9 Dec 2024 20:57:52 +0100 Subject: [PATCH 13/14] Revise CI --- .github/workflows/ci.yml | 58 +++++++++++++++------------------------- 1 file changed, 22 insertions(+), 36 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index c8ad725..38daf06 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -9,11 +9,12 @@ on: - main concurrency: - group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }} - cancel-in-progress: true + group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.sha }} + cancel-in-progress: true env: BUILD_TYPE: RelWithDebInfo + PRECICE_VERSION: v3.1.2 jobs: build: @@ -23,54 +24,39 @@ jobs: matrix: os: [macos-latest, ubuntu-latest] steps: - - uses: actions/checkout@master + - uses: actions/checkout@v3 with: repository: gismo/gismo ref: stable path: ./gismo - - uses: actions/checkout@master + - uses: actions/checkout@v3 with: path: ./gismo/optional/${{ github.event.repository.name }} - #token: ${{ secrets.GH_PAT }} - # Step to install preCICE - name: Install preCICE - if: runner.os == 'Linux' run: | - sudo apt update - sudo apt install -y build-essential cmake libeigen3-dev libxml2-dev libboost-all-dev petsc-dev python3-dev python3-numpy - wget https://github.com/precice/precice/archive/v3.1.2.tar.gz - tar -xzvf v3.1.2.tar.gz - cd precice-3.1.2 - cmake --preset=production # Configure using the production preset - cd build + if [[ "$RUNNER_OS" == "Linux" ]]; then + sudo apt update + sudo apt install -y build-essential cmake libeigen3-dev libxml2-dev libboost-all-dev petsc-dev python3-dev python3-numpy + elif [[ "$RUNNER_OS" == "macOS" ]]; then + brew install cmake eigen libxml2 boost petsc openmpi python3 numpy + fi + wget https://github.com/precice/precice/archive/$PRECICE_VERSION.tar.gz + tar -xzvf $PRECICE_VERSION.tar.gz + cd precice-${PRECICE_VERSION/v/} + cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=/usr/local .. make -j 10 sudo make install - precice_DIR=$(pwd) + export precice_DIR=$(pwd)/.. - # Step to install preCICE - - name: Install preCICE for mac - if: runner.os == 'macOS' - run: | - brew install cmake eigen libxml2 boost petsc openmpi python3 numpy - wget https://github.com/precice/precice/archive/v3.1.2.tar.gz - tar -xzvf v3.1.2.tar.gz - cd precice-3.1.2 - cmake --preset=production # Configure using the production preset - cd build - make -j 10 - sudo make install - precice_DIR=$(pwd) - - - name: "Run for ${{ matrix.os }}" + - name: Run Tests shell: bash - working-directory: ${{ runner.workspace }} + working-directory: ${{ github.workspace }} run: | - # Use the precice_DIR set from the environment variable - ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake \ - -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" \ + ctest -S ${{ github.repository }}/gismo/cmake/ctest_script.cmake \ + -D CTEST_BUILD_NAME="${{ github.repository }}_actions_$GITHUB_RUN_NUMBER" \ -D CTEST_SITE="${{ matrix.os }}_[actions]" \ - -D LABELS_FOR_SUBPROJECTS="gsPreCICE-examples" \ + -D LABELS_FOR_SUBPROJECTS="gsPreCICE-examples" \ -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON;-Dprecice_DIR=$precice_DIR" \ - -D GISMO_OPTIONAL="gsElasticity\\;${{ github.event.repository.name }}" -Q + -D GISMO_OPTIONAL="gsElasticity\\;${{ github.repository }}" From 4308f161d15f2cd158599616917fad4e5162afa6 Mon Sep 17 00:00:00 2001 From: Jingya Li <70752306+Crazy-Rich-Meghan@users.noreply.github.com> Date: Mon, 9 Dec 2024 21:01:38 +0100 Subject: [PATCH 14/14] Export preCICE_DIR for CI --- .github/workflows/ci.yml | 60 +++++++++++++++++++++++++--------------- 1 file changed, 38 insertions(+), 22 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 38daf06..cb13697 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -9,12 +9,11 @@ on: - main concurrency: - group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.sha }} - cancel-in-progress: true + group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }} + cancel-in-progress: true env: BUILD_TYPE: RelWithDebInfo - PRECICE_VERSION: v3.1.2 jobs: build: @@ -24,39 +23,56 @@ jobs: matrix: os: [macos-latest, ubuntu-latest] steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@master with: repository: gismo/gismo ref: stable path: ./gismo - - uses: actions/checkout@v3 + - uses: actions/checkout@master with: path: ./gismo/optional/${{ github.event.repository.name }} + #token: ${{ secrets.GH_PAT }} + # Step to install preCICE - name: Install preCICE + if: runner.os == 'Linux' run: | - if [[ "$RUNNER_OS" == "Linux" ]]; then - sudo apt update - sudo apt install -y build-essential cmake libeigen3-dev libxml2-dev libboost-all-dev petsc-dev python3-dev python3-numpy - elif [[ "$RUNNER_OS" == "macOS" ]]; then - brew install cmake eigen libxml2 boost petsc openmpi python3 numpy - fi - wget https://github.com/precice/precice/archive/$PRECICE_VERSION.tar.gz - tar -xzvf $PRECICE_VERSION.tar.gz - cd precice-${PRECICE_VERSION/v/} - cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=/usr/local .. + sudo apt update + sudo apt install -y build-essential cmake libeigen3-dev libxml2-dev libboost-all-dev petsc-dev python3-dev python3-numpy + wget https://github.com/precice/precice/archive/v3.1.2.tar.gz + tar -xzvf v3.1.2.tar.gz + cd precice-3.1.2 + cmake --preset=production # Configure using the production preset + cd build make -j 10 sudo make install - export precice_DIR=$(pwd)/.. + precice_DIR=$(pwd) + export precice_DIR=$(pwd) - - name: Run Tests + # Step to install preCICE + - name: Install preCICE for mac + if: runner.os == 'macOS' + run: | + brew install cmake eigen libxml2 boost petsc openmpi python3 numpy + wget https://github.com/precice/precice/archive/v3.1.2.tar.gz + tar -xzvf v3.1.2.tar.gz + cd precice-3.1.2 + cmake --preset=production # Configure using the production preset + cd build + make -j 10 + sudo make install + precice_DIR=$(pwd) + export precice_DIR=$(pwd) + + - name: "Run for ${{ matrix.os }}" shell: bash - working-directory: ${{ github.workspace }} + working-directory: ${{ runner.workspace }} run: | - ctest -S ${{ github.repository }}/gismo/cmake/ctest_script.cmake \ - -D CTEST_BUILD_NAME="${{ github.repository }}_actions_$GITHUB_RUN_NUMBER" \ + # Use the precice_DIR set from the environment variable + ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake \ + -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" \ -D CTEST_SITE="${{ matrix.os }}_[actions]" \ - -D LABELS_FOR_SUBPROJECTS="gsPreCICE-examples" \ + -D LABELS_FOR_SUBPROJECTS="gsPreCICE-examples" \ -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON;-Dprecice_DIR=$precice_DIR" \ - -D GISMO_OPTIONAL="gsElasticity\\;${{ github.repository }}" + -D GISMO_OPTIONAL="gsElasticity\\;${{ github.event.repository.name }}" -Q