1 /// This is the structural solver benchmark CSM3 from this paper:
2 /// "Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow"
3 /// Stefan Turek and Jaroslav Hron, <Fluid-Structure Interaction>, 2006
4 ///
5 /// Author: A.Shamanskiy (2016 - ...., TU Kaiserslautern)
6 #include <gismo.h>
7 #include <gsElasticity/gsElasticityAssembler.h>
8 #include <gsElasticity/gsMassAssembler.h>
9 #include <gsElasticity/gsElTimeIntegrator.h>
10 #include <gsElasticity/gsWriteParaviewMultiPhysics.h>
11 
12 using namespace gismo;
13 
writeLog(std::ofstream & file,const gsMultiPatch<> & displacement,real_t simTime,real_t compTime,real_t numIters)14 void writeLog(std::ofstream & file, const gsMultiPatch<> & displacement,
15               real_t simTime, real_t compTime, real_t numIters)
16 {
17     // evaluating displacement at the point A
18     gsMatrix<> point(2,1);
19     point << 1., 0.5;
20     gsMatrix<> dispA = displacement.patch(0).eval(point);
21 
22     // print: simTime dispAx dispAy compTime numIters
23     file << simTime << " " << dispA.at(0) << " " << dispA.at(1) << " "
24          << compTime << " " << numIters << std::endl;
25 }
26 
main(int argc,char * argv[])27 int main(int argc, char* argv[]){
28     gsInfo << "Benchmark CSM3: dynamic deflection of an elastic beam.\n";
29 
30     //=====================================//
31                 // Input //
32     //=====================================//
33 
34     std::string filename = ELAST_DATA_DIR"/flappingBeam_beam.xml";
35     real_t poissonsRatio = 0.4;
36     real_t youngsModulus = 1.4e6;
37     real_t density = 1.0e3;
38     real_t loading = 2.;
39     // space discretization
40     index_t numUniRef = 3;
41     index_t numDegElev = 0;
42     // time integration
43     real_t timeSpan = 2;
44     real_t timeStep = 0.01;
45     // output
46     index_t numPlotPoints = 1000;
47 
48     // minimalistic user interface for terminal
49     gsCmdLine cmd("Benchmark CSM3: dynamic deflection of an elastic beam.");
50     cmd.addReal("l","load","Gravitational loading acting on the beam",loading);
51     cmd.addInt("r","refine","Number of uniform refinement application",numUniRef);
52     cmd.addInt("d","degelev","Number of degree elevation application",numDegElev);
53     cmd.addReal("t","time","Time span, sec",timeSpan);
54     cmd.addReal("s","step","Time step, sec",timeStep);
55     cmd.addInt("p","points","Number of sampling points to plot to Paraview",numPlotPoints);
56     try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
57 
58     //=============================================//
59         // Scanning geometry and creating bases //
60     //=============================================//
61 
62     // scanning geometry
63     gsMultiPatch<> geometry;
64     gsReadFile<>(filename, geometry);
65 
66     // creating bases
67     gsMultiBasis<> basisDisplacement(geometry);
68     for (index_t i = 0; i < numDegElev; ++i)
69         basisDisplacement.degreeElevate();
70     for (index_t i = 0; i < numUniRef; ++i)
71         basisDisplacement.uniformRefine();
72 
73     //=============================================//
74         // Setting loads and boundary conditions //
75     //=============================================//
76 
77     // boundary conditions
78     gsBoundaryConditions<> bcInfo; // numbers are: patch, function pointer (nullptr) for displacement, displacement component
79     bcInfo.addCondition(0,boundary::west,condition_type::dirichlet,0,0);
80     bcInfo.addCondition(0,boundary::west,condition_type::dirichlet,0,1);
81 
82     // gravity, rhs
83     gsConstantFunction<> gravity(0.,loading*density,2);
84 
85     //=============================================//
86           // Setting assemblers and solvers //
87     //=============================================//
88 
89     // creating stiffness assembler
90     gsElasticityAssembler<real_t> assembler(geometry,basisDisplacement,bcInfo,gravity);
91     assembler.options().setReal("YoungsModulus",youngsModulus);
92     assembler.options().setReal("PoissonsRatio",poissonsRatio);
93     assembler.options().setInt("MaterialLaw",material_law::saint_venant_kirchhoff);
94     gsInfo << "Initialized system with " << assembler.numDofs() << " dofs.\n";
95 
96     // creating mass assembler
97     gsMassAssembler<real_t> massAssembler(geometry,basisDisplacement,bcInfo,gravity);
98     massAssembler.options().setReal("Density",density);
99 
100     // creating time integrator
101     gsElTimeIntegrator<real_t> timeSolver(assembler,massAssembler);
102     timeSolver.options().setInt("Scheme",time_integration::implicit_nonlinear);
103     timeSolver.options().setInt("Verbosity",solver_verbosity::none);
104 
105     //=============================================//
106             // Setting output & auxilary//
107     //=============================================//
108 
109     // displacement field
110     gsMultiPatch<> displacement;
111     // stress field
112     gsPiecewiseFunction<> stresses;
113     // constructing an IGA field (geometry + solution)
114     gsField<> dispField(geometry,displacement);
115     gsField<> stressField(assembler.patches(),stresses,true);
116     std::map<std::string,const gsField<> *> fields;
117     fields["Displacement"] = &dispField;
118     fields["von Mises"] = &stressField;
119 
120     std::ofstream logFile;
121     logFile.open("flappingBeam_CSM3.txt");
122     logFile << "# simTime dispAx dispAy compTime numIters\n";
123 
124     gsProgressBar bar;
125     gsStopwatch iterClock, totalClock;
126 
127     //=============================================//
128                    // Initial conditions //
129     //=============================================//
130 
131     // set initial conditions
132     timeSolver.setDisplacementVector(gsMatrix<>::Zero(assembler.numDofs(),1));
133     timeSolver.setVelocityVector(gsMatrix<>::Zero(assembler.numDofs(),1));
134 
135     timeSolver.constructSolution(displacement);
136     assembler.constructCauchyStresses(displacement,stresses,stress_components::von_mises);
137     writeLog(logFile,displacement,0.,0.,0);
138     // plotting initial displacement
139     gsParaviewCollection collection("flappingBeam_CSM3");
140     if (numPlotPoints > 0)
141         gsWriteParaviewMultiPhysicsTimeStep(fields,"flappingBeam_CSM3",collection,0,numPlotPoints);
142 
143     //=============================================//
144                   // Solving //
145     //=============================================//
146 
147     real_t simTime = 0.;
148     real_t numTimeStep = 0;
149     real_t compTime = 0.;
150 
151     gsInfo << "Running the simulation...\n";
152     totalClock.restart();
153     while (simTime < timeSpan)
154     {
155         bar.display(simTime/timeSpan);
156         iterClock.restart();
157 
158         timeSolver.makeTimeStep(timeStep);
159         timeSolver.constructSolution(displacement);
160 
161         //assembler.constructCauchyStresses(displacement,stresses,stress_components::von_mises);
162 
163         compTime += iterClock.stop();
164         simTime += timeStep;
165         numTimeStep++;
166 
167         if (numPlotPoints > 0)
168             gsWriteParaviewMultiPhysicsTimeStep(fields,"flappingBeam_CSM3",collection,numTimeStep,numPlotPoints);
169         writeLog(logFile,displacement,simTime,compTime,timeSolver.numberIterations());
170     }
171 
172     //=============================================//
173                 // Final touches //
174     //=============================================//
175 
176     gsInfo << "Simulation time: " + secToHMS(compTime) << " (total time: " + secToHMS(totalClock.stop()) + ")\n";
177 
178     if (numPlotPoints > 0)
179     {
180         collection.save();
181         gsInfo << "Open \"flappingBeam_CSM3.pvd\" in Paraview for visualization.\n";
182     }
183 
184     logFile.close();
185     gsInfo << "Log file created in \"flappingBeam_CSM3.txt\".\n";
186 
187     return 0;
188 }
189