ElmerFoamFSI  2.0
ElmerFoamFSI is fluid-solid interaction simulation application built up from OpenFOAM CFD and Elmer CSM coupled through the IMPACT multiphysics software integration infrastructure.
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Groups Pages
fsiFOAM.C
Go to the documentation of this file.
1 #include "fsiFOAM.H"
2 
3 using namespace COM;
4 
5 int fsifoam_module::Initialize(int argc,char *argv[])
6 {
7  sumLocalContErr = 0;
8  globalContErr = 0;
9  cumulativeContErr = 0;
10  pRefCell = 0;
11  pRefValue = 0.0;
12 
13  argsPtr = new Foam::argList(argc,argv);
14  if(!argsPtr){
15  Info << "FsiFoam:Initialize: Error: Failed to create args object." << endl;
16  return(1);
17  }
18  if (!argsPtr->checkRootCase()){
19  Info << "FsiFoam:Initialize: Error: Root case setup failed." << endl;
20  return(1);
21  }
22  Foam::argList &args(*argsPtr);
23  runTimePtr = new Foam::Time
24  (
25  Foam::Time::controlDictName,
26  args.rootPath(),
27  args.caseName(),
28  "system",
29  "constant",
30  !args.optionFound("noFunctionObjects")
31  );
32  if(!runTimePtr){
33  Info << "FsiFoam:Initialize: Error: Failed to setup time object." << endl;
34  return(1);
35  }
36 
37  InitFluidMesh();
38  InitTransportProperties();
39  CreateFluidFields();
40  InitStructuresMesh();
41  CreateStructuresFields();
42  ReadCouplingProperties();
43  CreateInterZoneInterpolators();
44  FindGlobalFaceZones();
45 
46  // functions specific to IMPACT FSI
47  CreateFSISurfaceMesh();
48 
49  return(0);
50 };
51 
52 int fsifoam_module::InitFluidMesh(){
53  // Foam::argList &args(*argsPtr);
54  Foam::Time &runTime(*runTimePtr);
55  Info << "FsiFoam:Initialize: Create dynamic mesh for time = "
56  << runTime.timeName() << endl;
57  // Info << OutStream.str() << endl;
58  // OutStream.clear();
59  // OutStream.str("");
60 
61  meshPtr = dynamicFvMesh::New(IOobject(dynamicFvMesh::defaultRegion,
62  runTime.timeName(),
63  runTime,
64  IOobject::MUST_READ));
65  return(0);
66 };
67 
68 int fsifoam_module::InitTransportProperties(){
69  Foam::Time &runTime(*runTimePtr);
70  dynamicFvMesh &mesh = meshPtr();
71  Info << "FsiFoam:Initialize: Reading transportProperties" << endl;
72  transportPropertiesPtr = new IOdictionary(IOobject("transportProperties",
73  runTime.constant(),
74  mesh,
75  IOobject::MUST_READ,
76  IOobject::NO_WRITE));
77  return(0);
78 };
79 
80 int fsifoam_module::CreateFluidFields(){
81  // Foam::argList &args(*argsPtr);
82  Foam::Time &runTime(*runTimePtr);
83  dynamicFvMesh &fluidsMesh = meshPtr();
84  IOdictionary &transportProperties(*transportPropertiesPtr);
85 
86 
87  nuPtr = new dimensionedScalar(transportProperties.lookup("nu"));
88  rhoFluidPtr = new dimensionedScalar(transportProperties.lookup("rho"));
89 
90  Info << "FsiFoam:CreateFluidFields: Reading field p" << endl;
91  // Info << OutStream.str() << endl;
92  // OutStream.clear();
93  // OutStream.str("");
94 
95  pPtr = new volScalarField(IOobject("p",
96  runTime.timeName(),
97  fluidsMesh,
98  IOobject::MUST_READ,
99  IOobject::AUTO_WRITE),
100  fluidsMesh);
101 
102  volScalarField &p(*pPtr);
103 
104  Info << "FsiFoam:CreateFluidFields: Reading field U" << endl;
105  // Info << OutStream.str() << endl;
106  // OutStream.clear();
107  // OutStream.str("");
108 
109  UPtr = new volVectorField(IOobject("U",runTime.timeName(),fluidsMesh,
110  IOobject::MUST_READ,IOobject::AUTO_WRITE),fluidsMesh);
111 
112  volVectorField &U = *UPtr;
113 
114  //# include "createPhi.H"
115  Info << "FsiFoam:CreateFluidFields: Reading/calculating face flux field phi" << endl;
116 
117  phiPtr = new surfaceScalarField(IOobject("phi",runTime.timeName(),fluidsMesh,
118  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
119  linearInterpolate(U) & fluidsMesh.Sf());
120 
121 
122 
123  setRefCell(p, fluidsMesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
124  return(0);
125 };
126 
127 int fsifoam_module::InitStructuresMesh(){
128  Foam::Time &runTime(*runTimePtr);
129  stressMeshPtr = new fvMesh(IOobject("solid",runTime.timeName(),runTime,IOobject::MUST_READ));
130  return(stressMeshPtr != NULL);
131 };
132 
133 int fsifoam_module::CreateStructuresFields(){
134  Foam::Time &runTime(*runTimePtr);
135  fvMesh &structMesh(*stressMeshPtr);
136 
137  Info << "FsiFoam:CreateStructuresFields: Reading incremental displacement field DU" << endl;
138  // Info << OutStream.str() << endl;
139  // OutStream.clear();
140  // OutStream.str("");
141 
142  DUPtr = new volVectorField(IOobject("DU",runTime.timeName(),structMesh,
143  IOobject::MUST_READ,IOobject::AUTO_WRITE),
144  structMesh);
145 
146  volVectorField &DURef(*DUPtr);
147 
148  // Info << "DU(0): " << DU << endl;
149 
150  gradDUPtr = new volTensorField(fvc::grad(DURef));
151  volTensorField &gradDURef(*gradDUPtr);
152 
153  UsolidPtr = new volVectorField(IOobject("Usolid",runTime.timeName(),structMesh,
154  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
155  DURef);
156 
157 
158  Info << "FsiFoam:CreateStructuresFields: Reading incremental displacement field DV" << endl;
159  DVPtr = new volVectorField(IOobject("DV",runTime.timeName(),structMesh,
160  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
161  fvc::ddt(DURef));
162 
163 
164  Info << "FsiFoam:CreateStructuresFields: Reading accumulated velocity field V" << endl;
165  VsPtr = new volVectorField(IOobject("Vs",runTime.timeName(),structMesh,
166  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
167  structMesh,
168  dimensionedVector("zero",dimVelocity,vector::zero));
169 
170 
171  Info << "FsiFoam:CreateStructuresFields: Reading accumulated stress field sigma" << endl;
172  sigmaPtr = new volSymmTensorField(IOobject("sigma",runTime.timeName(),structMesh,
173  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
174  structMesh,
175  dimensionedSymmTensor("zero",dimForce/dimArea,symmTensor::zero));
176  volSymmTensorField &sigmaref(*sigmaPtr);
177 
178 
179  Info << "FsiFoam:CreateStructuresFields: Reading incremental stress field DSigma" << endl;
180  DSigmaPtr = new volSymmTensorField(IOobject("DSigma",runTime.timeName(),structMesh,
181  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
182  structMesh,
183  dimensionedSymmTensor("zero",dimForce/dimArea,symmTensor::zero));
184 
185 
186  rheologyPtr = new constitutiveModel(sigmaref,DURef);
187 
188  constitutiveModel &rheology(*rheologyPtr);
189 
190  // volScalarField &rhoTest(rheology.rho());
191 
192  rhoPtr = new volScalarField(rheology.rho());
193  muPtr = new volScalarField(rheology.mu());
194  lambdaPtr = new volScalarField(rheology.lambda());
195 
196  FPtr = new volTensorField(I + gradDURef.T());
197  volTensorField &FRef(*FPtr);
198  DFPtr = new volTensorField(FRef - I);
199  JPtr = new volScalarField(det(FRef));
200 
201  solidDdtSchemePtr = new word(structMesh.schemesDict().ddtScheme("ddt("+DURef.name()+')'));
202  return(0);
203 };
204 int fsifoam_module::ReadCouplingProperties(){
205  Foam::argList &args(*argsPtr);
206  Foam::Time &runTime(*runTimePtr);
207  dynamicFvMesh &fluidMesh(meshPtr());
208  fvMesh &structMesh(*stressMeshPtr);
209 
210  Info << "FsiFoam:ReadCouplingProperties: Reading coupling properties" << endl;
211 
212  couplingPropertiesPtr = new IOdictionary(IOobject("couplingProperties",runTime.constant(),
213  fluidMesh,IOobject::MUST_READ,IOobject::NO_WRITE));
214  IOdictionary &couplingProperties(*couplingPropertiesPtr);
215 
216  solidPatchName = word(couplingProperties.lookup("solidPatch"));
217  solidPatchID = structMesh.boundaryMesh().findPatchID(solidPatchName);
218 
219  if (solidPatchID < 0)
220  {
221  FatalErrorIn(args.executable())
222  << "FsiFoam:ReadCouplingProperties: Problem with finding solid patch"
223  << abort(FatalError);
224  }
225 
226  solidZoneName = word(couplingProperties.lookup("solidZone"));
227 
228 
229  solidZoneID =
230  structMesh.faceZones().findZoneID(solidZoneName);
231 
232  if (solidZoneID < 0)
233  {
234  FatalErrorIn(args.executable())
235  << "FsiFoam:ReadCouplingProperties: Problem with finding solid zone"
236  << abort(FatalError);
237  }
238 
239  fluidPatchName = word(couplingProperties.lookup("fluidPatch"));
240  fluidPatchID = fluidMesh.boundaryMesh().findPatchID(fluidPatchName);
241 
242 
243  if (fluidPatchID < 0)
244  {
245  FatalErrorIn(args.executable())
246  << "FsiFoam:ReadCouplingProperties: Problem with finding fluid patch"
247  << abort(FatalError);
248  }
249 
250  fluidZoneName = word(couplingProperties.lookup("fluidZone"));
251  fluidZoneID = fluidMesh.faceZones().findZoneID(fluidZoneName);
252  if (fluidZoneID < 0)
253  {
254  FatalErrorIn(args.executable())
255  << "FsiFoam:ReadCouplingProperties: Problem with finding fluid zone"
256  << abort(FatalError);
257  }
258 
259  feMotionSolver = fluidMesh.objectRegistry::foundObject<tetPointVectorField>("motionU");
260  fvMotionSolver = fluidMesh.objectRegistry::foundObject<pointVectorField>("pointMotionU");
261 
262  // Info << OutStream.str() << endl;
263  // OutStream.clear();
264  // OutStream.str("");
265 
266 
267  // Grab solid patch field
268  Info << "FsiFoam:ReadCouplingProperties:Solid patch ID: " << solidPatchID << endl;
269  volVectorField &DURef(*DUPtr);
270  if
271  (
272  DURef.boundaryField()[solidPatchID].type()
273  != solidTractionFvPatchVectorField::typeName
275  )
276  {
277  FatalErrorIn(args.executable())
278  << "FsiFoam:ReadCouplingProperties: Boundary condition on " << DURef.name()
279  << " is "
280  << DURef.boundaryField()[solidPatchID].type()
281  << "for fluid -solid interface patch, instead "
282  << solidTractionFvPatchVectorField::typeName
283  //<< tractionDisplacementIncrementFvPatchVectorField::typeName
284  << abort(FatalError);
285  }
286 
287  //tractionDisplacementIncrementFvPatchVectorField& tForce =
288  // refCast<tractionDisplacementIncrementFvPatchVectorField>
289  // tForcePtr = new solidTractionFvPatchVectorField(refCast<solidTractionFvPatchVectorField>(DU.boundaryField()[solidPatchID]));
290  // solidTractionFvPatchVectorField& tForce = *tForcePtr;
291 
292 
293  accumulatedFluidInterfaceDisplacementHeaderPtr = new IOobject("accumulatedFluidInterfaceDisplacement",
294  runTime.timeName(),
295  fluidMesh,
296  IOobject::MUST_READ);
297  IOobject &accumulatedFluidInterfaceDisplacementHeader(*accumulatedFluidInterfaceDisplacementHeaderPtr);
298 
299 
300 
301  if(accumulatedFluidInterfaceDisplacementHeader.headerOk())
302  {
303  Pout << "FsiFoam:ReadCouplingProperties: Reading accumulated fluid interface displacement" << endl;
304 
305  accumulatedFluidInterfaceDisplacementPtr =
306  new vectorIOField
307  (
308  IOobject
309  (
310  "accumulatedFluidInterfaceDisplacement",
311  runTime.timeName(),
312  fluidMesh,
313  IOobject::MUST_READ,
314  IOobject::AUTO_WRITE
315  )
316  );
317  }
318  else
319  {
320  accumulatedFluidInterfaceDisplacementPtr =
321  new vectorIOField
322  (
323  IOobject
324  (
325  "accumulatedFluidInterfaceDisplacement",
326  runTime.timeName(),
327  fluidMesh,
328  IOobject::NO_READ,
329  IOobject::AUTO_WRITE
330  ),
331  vectorField
332 
333  (
334  fluidMesh.boundaryMesh()[fluidPatchID].nPoints(),
335  vector::zero
336  )
337  );
338  }
339 
340 
341  // sumLocalContErr = 0;
342  // globalContErr = 0;
343  // cumulativeContErr = 0;
344 
345 
346 
347  return(0);
348 };
349 
350 int fsifoam_module::CreateInterZoneInterpolators()
351 {
352 
353  dynamicFvMesh &fluidMesh = meshPtr();
354  fvMesh &structMesh(*stressMeshPtr);
355 
356  if(!interpolatorFluidSolidPtr || !interpolatorSolidFluidPtr)
357  {
358  deleteDemandDrivenData(interpolatorFluidSolidPtr);
359  deleteDemandDrivenData(interpolatorSolidFluidPtr);
360 
361  Info << "FsiFoam:CreateInterZoneInterpolators: Create fluid-to-solid and solid-to-fluid interpolators" << endl;
362 
363  interpolatorFluidSolidPtr = new zoneToZoneInterpolation
364  (
365  fluidMesh.faceZones()[fluidZoneID](),
366  structMesh.faceZones()[solidZoneID](),
367  intersection::VISIBLE
368  );
369 
370  interpolatorSolidFluidPtr = new zoneToZoneInterpolation
371  (
372  structMesh.faceZones()[solidZoneID](),
373  fluidMesh.faceZones()[fluidZoneID](),
374  intersection::VISIBLE
375  );
376 
377  Info << "FsiFoam:CreateInterZoneInterpolators: Check fluid-to-solid and solid-to-fluid interpolators" << endl;
378 
379  {
380  vectorField fluidPatchFaceCentres =
381  vectorField(fluidMesh.boundaryMesh()[fluidPatchID].faceCentres());
382 
383  vectorField fluidZoneFaceCentres
384  (
385  fluidMesh.faceZones()[fluidZoneID].size(),
386  vector::zero
387  );
388 
389 
390  const label fluidPatchStart =
391  fluidMesh.boundaryMesh()[fluidPatchID].start();
392 
393  forAll (fluidPatchFaceCentres, i)
394  {
395  fluidZoneFaceCentres
396  [
397  fluidMesh.faceZones()[fluidZoneID].whichFace(fluidPatchStart + i)
398  ] =
399  fluidPatchFaceCentres[i];
400  }
401 
402  // Parallel data exchange: collect faceCentres field on all processors
403  reduce(fluidZoneFaceCentres, sumOp<vectorField>());
404 
405  vectorField solidZoneFaceCentres =
406  interpolatorFluidSolidPtr->faceInterpolate
407  (
408  fluidZoneFaceCentres
409  );
410 
411  vectorField solidPatchFaceCentres
412  (
413  structMesh.boundaryMesh()[solidPatchID].size(),
414  vector::zero
415  );
416 
417  const label solidPatchStart =
418  structMesh.boundaryMesh()[solidPatchID].start();
419 
420  forAll(solidPatchFaceCentres, i)
421  {
422  solidPatchFaceCentres[i] =
423  solidZoneFaceCentres
424  [
425  structMesh.faceZones()[solidZoneID]
426  .whichFace(solidPatchStart + i)
427  ];
428  }
429 
430  scalar maxDist = gMax
431  (
432  mag
433  (
434  solidPatchFaceCentres
435  - structMesh.boundaryMesh()[solidPatchID].faceCentres()
436  )
437  );
438 
439  Info << "FsiFoam:CreateInterZoneInterpolators: Fluid-to-solid face interpolation error: " << maxDist
440  << endl;
441  }
442 
443 
444 
445  {
446  vectorField solidPatchFaceCentres =
447  vectorField(structMesh.boundaryMesh()[solidPatchID].faceCentres());
448 
449  vectorField solidZoneFaceCentres
450  (
451  structMesh.faceZones()[solidZoneID].size(),
452  vector::zero
453  );
454 
455  const label solidPatchStart =
456  structMesh.boundaryMesh()[solidPatchID].start();
457 
458  forAll (solidPatchFaceCentres, i)
459  {
460  solidZoneFaceCentres
461  [
462  structMesh.faceZones()[solidZoneID]
463  .whichFace(solidPatchStart + i)
464  ] =
465  solidPatchFaceCentres[i];
466  }
467 
468  // Parallel data exchange: collect faceCentres field on all processors
469  reduce(solidZoneFaceCentres, sumOp<vectorField>());
470 
471  vectorField fluidZoneFaceCentres =
472  interpolatorSolidFluidPtr->faceInterpolate
473  (
474  solidZoneFaceCentres
475  );
476 
477  vectorField fluidPatchFaceCentres
478  (
479  fluidMesh.boundaryMesh()[fluidPatchID].size(),
480  vector::zero
481  );
482 
483  const label fluidPatchStart =
484  fluidMesh.boundaryMesh()[fluidPatchID].start();
485 
486  forAll(fluidPatchFaceCentres, i)
487  {
488  fluidPatchFaceCentres[i] =
489  fluidZoneFaceCentres
490  [
491  fluidMesh.faceZones()[fluidZoneID]
492  .whichFace(fluidPatchStart + i)
493  ];
494  }
495 
496  scalar maxDist = gMax
497  (
498  mag
499  (
500  fluidPatchFaceCentres
501  - fluidMesh.boundaryMesh()[fluidPatchID].faceCentres()
502  )
503  );
504 
505  Info << "FsiFoam:CreateInterZoneInterpolators: Solid-to-fluid face interpolation error: " << maxDist
506  << endl;
507  }
508 
509  }
510 
511  return(0);
512 };
513 
514 int fsifoam_module::FindGlobalFaceZones()
515 {
516 
517  fvMesh &structMesh(*stressMeshPtr);
518  Foam::argList &args = *argsPtr;
519 
520  {
521  SLList<label> globalFaceZonesSet;
522 
523  const faceZoneMesh& faceZones = structMesh.faceZones();
524 
525  forAll(faceZones, zoneI)
526  {
527  const faceZone& curFaceZone = faceZones[zoneI];
528 
529  forAll(curFaceZone, faceI)
530  {
531  // if unused face exist
532  if (curFaceZone[faceI] >= structMesh.nFaces())
533  {
534  globalFaceZonesSet.insert(zoneI);
535  break;
536  }
537  }
538  }
539 
540  globalFaceZones = labelList(globalFaceZonesSet);
541  }
542  globalToLocalFaceZonePointMap.resize(globalFaceZones.size());
543 
544  forAll(globalFaceZones, zoneI)
545  {
546  label curZoneID = globalFaceZones[zoneI];
547  labelList curMap(structMesh.faceZones()[curZoneID]().nPoints(), -1);
548  vectorField fzGlobalPoints =
549  structMesh.faceZones()[curZoneID]().localPoints();
550 
551  //- set all slave points to zero because only the master order is used
552  if(!Pstream::master())
553  {
554  fzGlobalPoints *= 0.0;
555  }
556 
557  //- pass points to all procs
558  reduce(fzGlobalPoints, sumOp<vectorField>());
559 
560  //- now every proc has the master's list of FZ points
561  //- every proc must now find the mapping from their local FZ points to
562  //- the global FZ points
563  const vectorField& fzLocalPoints =
564  structMesh.faceZones()[curZoneID]().localPoints();
565  const edgeList& fzLocalEdges =
566  structMesh.faceZones()[curZoneID]().edges();
567  const labelListList& fzPointEdges =
568  structMesh.faceZones()[curZoneID]().pointEdges();
569  scalarField minEdgeLength(fzLocalPoints.size(), GREAT);
570 
571  forAll(minEdgeLength, pI)
572  {
573  const labelList& curPointEdges = fzPointEdges[pI];
574  forAll(curPointEdges, eI)
575  {
576  scalar Le = fzLocalEdges[curPointEdges[eI]].mag(fzLocalPoints);
577  if (Le < minEdgeLength[pI])
578  {
579  minEdgeLength[pI] = Le;
580  }
581  }
582  }
583 
584  forAll(fzGlobalPoints, globalPointI)
585  {
586  // scalar minDist = GREAT;
587 
588  forAll(fzLocalPoints, procPointI)
589  {
590  scalar curDist =
591  mag
592  (
593  fzLocalPoints[procPointI]
594  - fzGlobalPoints[globalPointI]
595  );
596 
597  // if (curDist < minDist)
598  // {
599  // minDist = curDist;
600  // }
601 
602  if (curDist < 1e-4*minEdgeLength[procPointI])
603  {
604  curMap[globalPointI] = procPointI;
605  break;
606  }
607  }
608 
609  // if (curMap[globalPointI] == -1)
610  // {
611  // Pout << "minDist: " << minDist << endl;
612  // }
613  }
614 
615  forAll(curMap, globalPointI)
616  {
617  if (curMap[globalPointI] == -1)
618  {
619  FatalErrorIn(args.executable())
620  << "FsiFoam:FindGlobalFaceZones: local to global face zone point map is not correct"
621  << abort(FatalError);
622  }
623  }
624 
625  globalToLocalFaceZonePointMap[zoneI] = curMap;
626  }
627  return(0);
628 };
629 int fsifoam_module::ReadPISOControls()
630 {
631  dynamicFvMesh &fluidMesh(meshPtr());
632  dictionary &piso(fluidMesh.solutionDict().subDict("PISO"));
633 
634  nCorr = readInt(piso.lookup("nCorrectors"));
635  nNonOrthCorr = piso.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
636  momentumPredictor = piso.lookupOrDefault<Switch>("momentumPredictor", true);
637  transonic = piso.lookupOrDefault<Switch>("transonic", false);
638  nOuterCorr = piso.lookupOrDefault<int>("nOuterCorrectors", 1);
639 
640  return(0);
641 };
642 int fsifoam_module::ReadFSIControls()
643 {
644  IOdictionary &couplingProperties(*couplingPropertiesPtr);
645  // dynamicFvMesh &fluidMesh(meshPtr());
646  // Foam::Time &runTime(*runTimePtr);
647 
648  if (couplingProperties.found("couplingScheme"))
649  {
650  couplingScheme = word(couplingProperties.lookup("couplingScheme"));
651  if
652  (
653  (couplingScheme == "IQN-ILS")
654  || (couplingScheme == "Aitken")
655  || (couplingScheme == "FixedRelaxation")
656  || (couplingScheme == "NonIterative")
657  )
658  {
659  Info<< "FsiFoam:ReadFSIControl: Selecting coupling scheme " << couplingScheme << endl;
660  }
661  else
662  {
663  FatalErrorIn
664  (
665  "readFsiProperties"
666  ) << "couplingScheme: " << couplingScheme
667  << " is not a valid choice. "
668  << "Options are: IQN-ILS, Aitken, FixedRelaxation"
669  << abort(FatalError);
670  }
671  }
672 
673  interfaceDeformationLimit =
674  scalar(readScalar(couplingProperties.lookup("interfaceDeformationLimit")));
675 
676 
677  // if(dynamicMeshDictPtr)
678  // delete dynamicMeshDictPtr;
679  // dynamicMeshDictPtr = new IOdictionary(IOobject
680  // (
681  // "dynamicMeshDict",
682  // runTime.constant(),
683  // fluidMesh,
684  // IOobject::MUST_READ,
685  // IOobject::NO_WRITE,
686  // false
687  // ));
688 
689 
690 
691  fsiRelaxationFactorMin = scalar(readScalar(couplingProperties.lookup("fsiRelaxationFactor")));
692  fsiRelaxationFactor = fsiRelaxationFactorMin;
693  outerCorrTolerance = scalar(readScalar(couplingProperties.lookup("outerCorrTolerance")));
694  Info << "FsiFoam:ReadFSIControl: outerCorrTolerance = " << outerCorrTolerance << endl;
695  fsi = Switch(couplingProperties.lookup("fsi"));
696  if(fsi)
697  Info << "FsiFoam:ReadFSIControl: FSI IS ENABLED." << endl;
698  else
699  Info << "FsiFoam:ReadFSIControl: FSI IS DISABLED." << endl;
700  return(0);
701 };
702 
703 int fsifoam_module::Step(){
704  Foam::argList &args(this->ArgList());
705  Foam::Time &runTime(this->RunTime());
706  dynamicFvMesh &fluidsMesh(this->FluidMesh());
707  // IOdictionary &transportProperties(this->TransportProperties());
708  dimensionedScalar &nu(this->nu());
709  dimensionedScalar &rhoFluid(this->rhoFluid());
710  volScalarField &p(this->p());
711  volVectorField &U(this->U());
712  surfaceScalarField &phi(this->phi());
713  // label pRefCell = 0;
714  // scalar pRefValue = 0.0;
715  // setRefCell(p, fluidsMesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
716 
717  fvMesh &structuresMesh(this->StructuresMesh());
718 
719 
720  volVectorField &DU(this->DU());
721  volTensorField &gradDU(this->gradDU());
722  volVectorField &Usolid(this->Usolid());
723  volVectorField &DV(this->DV());
724  volVectorField &Vs(this->Vs());
725  volSymmTensorField &sigma(this->sigma());
726  volSymmTensorField &DSigma(this->DSigma());
727  // constitutiveModel &Rheology(){return(*rheologyPtr);};
728  volScalarField &rho(this->rhoSolid());
729  volScalarField &mu(this->mu());
730  volScalarField &lambda(this->lambda());
731  volTensorField &F(this->F());
732  volTensorField &DF(this->DF());
733  volScalarField &J(this->J());
734  word &solidDdtScheme(this->SolidDdtScheme());
735  // pointMesh &pStressMesh(){return(*pStressMeshPtr);};
736  // pointVectorField &PointDU(){return(*pointDUPtr);};
737 
738 
739  // IOdictionary &couplingProperties(this->CouplingProperties());
740  label solidPatchID(this->SolidPatchID());
741  // word solidPatchName(this->SolidPatchName());
742  label fluidPatchID(this->FluidPatchID());
743  // word fluidPatchName(this->FluidPatchName());
744  label fluidZoneID(this->FluidZoneID());
745  label solidZoneID(this->SolidZoneID());
746 
747  bool feMotionSolver(this->FEMotion());
748  bool fvMotionSolver(this->FVMotion());
749  vectorField &accumulatedFluidInterfaceDisplacement(this->AccumulatedFluidInterfaceDisplacements());
750  // Masoud:
751  Info << "FsiFoam:Step: accumulatedFluidInterfaceDisplacement = " << endl;
752  Info << accumulatedFluidInterfaceDisplacement << endl;
753  // Masoud: End
754  solidTractionFvPatchVectorField &tForce(this->tForce());
755 
756  zoneToZoneInterpolation &interpolatorFluidSolid(this->interpFluidSolid());
757  zoneToZoneInterpolation &interpolatorSolidFluid(this->interpSolidFluid());
758 
759  scalar &sumLocalContErr(this->LocalContErr());
760  scalar &globalContErr(this->GlobalContErr());
761  scalar &cumulativeContErr(this->CumulativeContErr());
762 
763  labelList &globalFaceZones(this->GlobalFaceZones());
764  labelListList &globalToLocalFaceZonePointMap(this->GlobalToLocalFaceZonePointMap());
765  runTime++;
766  if(!runTime.end()){
767  // for (runTime++; !runTime.end(); runTime++)
768  // {
769  // Info << "Time = " << runTime.timeName() << nl << endl;
770 
771  //# include "readPISOControls.H"
772  //# include "readFsiControls.H"
773  //# include "createStressPointMesh.H"
774 
775  this->ReadPISOControls();
776  this->ReadFSIControls();
777  word couplingScheme(this->CouplingScheme());
778  // label &outerCorr(this->OuterCorr());
779  scalar &fsiRelaxationFactor(this->FSIRelaxationFactor());
780  // scalar &fsiRelaxationFactorMin(this->FSIRelaxationFactorMin());
781  Switch fsi(this->FSIEnabled());
782  int &nCorrPISO(this->NCorrPISO());
783  int &nNonOrthCorr(this->NNonOrthCorr());
784  scalar &interfaceDeformationLimit(this->InterfaceDeformationLimit());
785  scalar &outerCorrTolerance(this->OuterCorrTolerance());
786  int &nOuterCorr(this->NOuterCorr());
787 
788  // this->CreateStructuresPointMesh();
789 
790 
791  // Create point stress mesh for interpolation
792  pointMesh StructuresPointMesh(structuresMesh);
793 
794  pointPatchInterpolation patchPointInterpolator(structuresMesh);
795 
796 
797  // Create point displacement field (check fixedValue patches)
798 
799  wordList types
800  (
801  StructuresPointMesh.boundary().size(),
802  calculatedFvPatchVectorField::typeName
803  );
804 
805  // wordList types = DU.boundaryField().types();
806 
807  forAll(DU.boundaryField().types(), patchI)
808  {
809  if
810  (
811  DU.boundaryField().types()[patchI]
812  == fixedValueFvPatchVectorField::typeName
813  )
814  {
815  types[patchI] = fixedValueFvPatchVectorField::typeName;
816  }
817  }
818 
819  pointVectorField pointDU
820  (
821  IOobject
822  (
823  "pointDU",
824  runTime.timeName(),
825  structuresMesh
826  ),
827  StructuresPointMesh,
828  dimensionedVector("zero", dimLength, vector::zero),
829  types
830  );
831  //# include "createInterfaceFields.H"
832  vectorField fluidPatchPointsDispl
833  (
834  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
835  vector::zero
836  );
837 
838  vectorField fluidPatchPointsDisplOld
839  (
840  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
841  vector::zero
842  );
843 
844  vectorField solidPatchPointsDispl
845  (
846  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
847  vector::zero
848  );
849 
850  vectorField fsiResidual
851  (
852  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
853  vector::zero
854  );
855 
856  vectorField fsiResidualOld
857  (
858  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
859  vector::zero
860  );
861 
862  vectorField solidPatchTraction
863  (
864  structuresMesh.boundary()[solidPatchID].size(),
865  vector::zero
866  );
867 
868  scalarField solidPatchPressure
869  (
870  structuresMesh.boundary()[solidPatchID].size(),
871  0.0
872  );
873 
874  scalar initialFsiResidualNorm = 0;
875  scalar fsiResidualNorm = 0;
876  label outerCorr=0;
877 
878  do
879  {
880  outerCorr++;
881 
882  //# include "setInterfaceDisplacement.H"
883  Info << "\nTime = " << runTime.timeName()
884  << ", iteration: " << outerCorr << endl;
885 
886  if (outerCorr < 3 || couplingScheme == "FixedRelaxation")
887  {
888  Info << "FsiFoam:Step: Current fsi under-relaxation factor: "
889  << fsiRelaxationFactor << endl;
890 
891  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
892 
893  fluidPatchPointsDispl += fsiRelaxationFactor*fsiResidual;
894  }
895  else
896  {
897  if (couplingScheme == "Aitken")
898  {
899  fsiRelaxationFactor =
900  -fsiRelaxationFactor
901  *(
902  gSum
903  (
904  fsiResidualOld
905  &(fsiResidual - fsiResidualOld)
906  )
907  /(
908  gSum
909  (
910  (fsiResidual - fsiResidualOld)
911  &(fsiResidual - fsiResidualOld)
912  )
913  )
914  );
915 
916  fsiRelaxationFactor = mag(fsiRelaxationFactor);
917 
918  Info << "FsiFoam:Step: Current fsi under-relaxation factor (Aitken): "
919  << fsiRelaxationFactor << endl;
920 
921  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
922 
923  fluidPatchPointsDispl +=
924  fsiRelaxationFactor*fsiResidual;
925  }
926  // else if (couplingScheme == "IQN-ILS")
927  // {
928 
929  // }
930  }
931 
932  //# include "moveFluidMesh.H"
933  {
934  // Move fluid mesh
935  const vectorField& n =
936  fluidsMesh.boundaryMesh()[fluidPatchID].pointNormals();
937 
938  primitivePatchInterpolation patchInterpolator
939  (
940  fluidsMesh.boundaryMesh()[fluidPatchID]
941  );
942 
943  scalarField pointDeltaCoeffs =
944  patchInterpolator.faceToPointInterpolate
945  (
946  fluidsMesh.boundary()[fluidPatchID].deltaCoeffs()
947  );
948 
949  scalar delta =
950  gMax
951  (
952  mag
953  (
954  n
955  & (
956  accumulatedFluidInterfaceDisplacement
957  + fluidPatchPointsDispl
958  - fluidPatchPointsDisplOld
959  )
960  )
961  *pointDeltaCoeffs
962  );
963 
964  Info << "FsiFoam:Step: Maximal accumulated displacement of interface points: "
965  << delta << endl;
966 
967  if(delta < interfaceDeformationLimit)
968  {
969  // Move only interface points
970  // Masoud
971  Info << "FsiFoam:Step: Moving only interface points..." << endl;
972  // Masoud: End
973  pointField newPoints = fluidsMesh.allPoints();
974 
975  const labelList& meshPoints =
976  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
977 
978  forAll (fluidPatchPointsDispl, pointI)
979  {
980  //Info << "1. meshPoints[" << pointI << "] = " << meshPoints[pointI]
981  //<< endl;
982  newPoints[meshPoints[pointI]] +=
983  fluidPatchPointsDispl[pointI]
984  - fluidPatchPointsDisplOld[pointI];
985  }
986 
987  twoDPointCorrector twoDCorrector(fluidsMesh);
988 
989  twoDCorrector.correctPoints(newPoints);
990 
991  fluidsMesh.movePoints(newPoints);
992 
993  // Accumulate interface points displacement
994  accumulatedFluidInterfaceDisplacement +=
995  fluidPatchPointsDispl
996  - fluidPatchPointsDisplOld;
997  }
998  else
999  {
1000  // Move whole fluid mesh
1001  // Masoud
1002  Info << "FsiFoam:Step: Moving the whole mesh .... " << endl;
1003  // Masoud: End
1004  pointField newPoints = fluidsMesh.allPoints();
1005 
1006  const labelList& meshPoints =
1007  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
1008 
1009  forAll (accumulatedFluidInterfaceDisplacement, pointI)
1010  {
1011  //Info << "2. meshPoints[" << pointI << "] = " << meshPoints[pointI]
1012  //<< endl;
1013  newPoints[meshPoints[pointI]] -=
1014  accumulatedFluidInterfaceDisplacement[pointI];
1015  }
1016 
1017  twoDPointCorrector twoDCorrector(fluidsMesh);
1018 
1019  twoDCorrector.correctPoints(newPoints);
1020 
1021  fluidsMesh.movePoints(newPoints);
1022 
1023  accumulatedFluidInterfaceDisplacement +=
1024  fluidPatchPointsDispl
1025  - fluidPatchPointsDisplOld;
1026 
1027 
1028  if (feMotionSolver)
1029  {
1030  tetPointVectorField& motionU =
1031  const_cast<tetPointVectorField&>
1032  (
1033  fluidsMesh.objectRegistry:: lookupObject<tetPointVectorField>
1034  (
1035  "motionU"
1036  )
1037  );
1038 
1039  fixedValueTetPolyPatchVectorField& motionUFluidPatch =
1040  refCast<fixedValueTetPolyPatchVectorField>
1041  (
1042  motionU.boundaryField()[fluidPatchID]
1043  );
1044 
1045  tetPolyPatchInterpolation tppi
1046  (
1047  refCast<const faceTetPolyPatch>(motionUFluidPatch.patch())
1048  );
1049 
1050  motionUFluidPatch ==
1051  tppi.pointToPointInterpolate
1052  (
1053  accumulatedFluidInterfaceDisplacement
1054  /runTime.deltaT().value()
1055  );
1056  }
1057  else if (fvMotionSolver)
1058  {
1059  pointVectorField& motionU =
1060  const_cast<pointVectorField&>
1061  (
1062  fluidsMesh.objectRegistry:: lookupObject<pointVectorField>
1063  (
1064  "pointMotionU"
1065  )
1066  );
1067 
1068  fixedValuePointPatchVectorField& motionUFluidPatch =
1069  refCast<fixedValuePointPatchVectorField>
1070  (
1071  motionU.boundaryField()[fluidPatchID]
1072  );
1073 
1074  motionUFluidPatch ==
1075  accumulatedFluidInterfaceDisplacement
1076  /runTime.deltaT().value();
1077  }
1078  else
1079  {
1080  FatalErrorIn(args.executable())
1081  << "FsiFoam:Step: Problem with mesh motion solver selection"
1082  << abort(FatalError);
1083  }
1084 
1085  fluidsMesh.update();
1086 
1087  accumulatedFluidInterfaceDisplacement = vector::zero;
1088  }
1089  }
1090 
1091  // Info << "DU after move fluid: " << DU << endl;
1092 
1093  //# include "solveFluid.H"
1094  if(fluidsMesh.moving())
1095  {
1096  // Make the fluxes relative
1097  phi -= fvc::meshPhi(U);
1098  }
1099 
1100 
1101  //# include "CourantNo.H"
1102  scalar CoNum = 0.0;
1103  scalar meanCoNum = 0.0;
1104  scalar velMag = 0.0;
1105 
1106  if (fluidsMesh.nInternalFaces())
1107  {
1108  surfaceScalarField magPhi = mag(phi);
1109 
1110  surfaceScalarField SfUfbyDelta =
1111  fluidsMesh.surfaceInterpolation::deltaCoeffs()*magPhi;
1112 
1113  const scalar deltaT = runTime.deltaT().value();
1114 
1115  CoNum = max(SfUfbyDelta/fluidsMesh.magSf()).value()*deltaT;
1116 
1117  meanCoNum = (sum(SfUfbyDelta)/sum(fluidsMesh.magSf())).value()*deltaT;
1118 
1119  velMag = max(magPhi/fluidsMesh.magSf()).value();
1120  }
1121 
1122  Info<< "FsiFoam:Step: Courant Number mean: " << meanCoNum
1123  << " max: " << CoNum
1124  << " velocity magnitude: " << velMag
1125  << endl;
1126 
1127  fvVectorMatrix UEqn
1128  (
1129  fvm::ddt(U)
1130  + fvm::div(phi, U)
1131  - fvm::laplacian(nu, U)
1132  );
1133 
1134  solve(UEqn == -fvc::grad(p));
1135 
1136  // --- PISO loop
1137  volScalarField rUA = 1.0/UEqn.A();
1138 
1139  for (int corr=0; corr<nCorrPISO; corr++)
1140  {
1141  U = rUA*UEqn.H();
1142  phi = (fvc::interpolate(U) & fluidsMesh.Sf());
1143 
1144  adjustPhi(phi, U, p);
1145 
1146  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
1147  {
1148  fvScalarMatrix pEqn
1149  (
1150  fvm::laplacian(rUA, p)
1151  == fvc::div(phi)
1152  );
1153 
1154  pEqn.setReference(pRefCell, pRefValue);
1155  pEqn.solve();
1156 
1157  if (nonOrth == nNonOrthCorr)
1158  {
1159  phi -= pEqn.flux();
1160  }
1161  }
1162 
1163  //# include "continuityErrs.H"
1164  {
1165  volScalarField contErr = fvc::div(phi);
1166 
1167  sumLocalContErr = runTime.deltaT().value()*
1168  mag(contErr)().weightedAverage(fluidsMesh.V()).value();
1169 
1170  globalContErr = runTime.deltaT().value()*
1171  contErr.weightedAverage(fluidsMesh.V()).value();
1172 
1173  cumulativeContErr += globalContErr;
1174 
1175  Info<< "FsiFoam:Step: time step continuity errors : sum local = " << sumLocalContErr
1176  << ", global = " << globalContErr
1177  << ", cumulative = " << cumulativeContErr
1178  << endl;
1179  }
1180 
1181  U -= rUA*fvc::grad(p);
1182  U.correctBoundaryConditions();
1183  }
1184 
1185  // Info << "DU after solve fluid: " << DU << endl;
1186 
1187  //# include "setInterfaceForce.H"
1188  {
1189  Info << "FsiFoam:Step: Setting traction on solid patch" << endl;
1190 
1191  // vectorField fluidPatchTraction =
1192  // -rhoFluid.value()*nu.value()
1193  // *U.boundaryField()[fluidPatchID].snGrad()
1194  // + rhoFluid.value()*p.boundaryField()[fluidPatchID]
1195  // *fluidsMesh.boundary()[fluidPatchID].nf();
1196 
1197  vectorField fluidPatchTraction =
1198  -rhoFluid.value()*nu.value()
1199  *U.boundaryField()[fluidPatchID].snGrad();
1200 
1201  // Info << "Fluid patch traction: " << fluidPatchTraction << endl;
1202 
1203  scalarField fluidPatchPressure =
1204  rhoFluid.value()*p.boundaryField()[fluidPatchID];
1205 
1206  // Info << "Fluid patch pressure: " << fluidPatchPressure << endl;
1207 
1208  vectorField fluidZoneTraction
1209  (
1210  fluidsMesh.faceZones()[fluidZoneID].size(),
1211  vector::zero
1212  );
1213 
1214  // Info << "Fluid zone traction: " << fluidZoneTraction << endl;
1215 
1216  const label fluidPatchStart =
1217  fluidsMesh.boundaryMesh()[fluidPatchID].start();
1218 
1219  forAll(fluidPatchTraction, i)
1220  {
1221  fluidZoneTraction
1222  [
1223  fluidsMesh.faceZones()[fluidZoneID].whichFace(fluidPatchStart + i)
1224  ] =
1225  fluidPatchTraction[i];
1226  }
1227 
1228  // Parallel data exchange: collect pressure field on all processors
1229  reduce(fluidZoneTraction, sumOp<vectorField>());
1230 
1231 
1232  scalarField fluidZonePressure
1233  (
1234  fluidsMesh.faceZones()[fluidZoneID].size(),
1235  0.0
1236  );
1237 
1238  forAll(fluidPatchPressure, i)
1239  {
1240  fluidZonePressure
1241  [
1242  fluidsMesh.faceZones()[fluidZoneID].whichFace(fluidPatchStart + i)
1243  ] =
1244  fluidPatchPressure[i];
1245  }
1246 
1247  // Parallel data exchange: collect pressure field on all processors
1248  reduce(fluidZonePressure, sumOp<scalarField>());
1249 
1250  vectorField solidZoneTraction =
1251  interpolatorFluidSolid.faceInterpolate
1252  (
1253  fluidZoneTraction
1254  );
1255 
1256  // Info << "Solid Zone Traction: " << solidZoneTraction << endl;
1257 
1258  scalarField solidZonePressure =
1259  interpolatorFluidSolid.faceInterpolate
1260  (
1261  fluidZonePressure
1262  );
1263 
1264  // Info << "Solid Zone Pressure: " << solidZonePressure << endl;
1265 
1266  const label solidPatchStart =
1267  structuresMesh.boundaryMesh()[solidPatchID].start();
1268 
1269  forAll(solidPatchTraction, i)
1270  {
1271  solidPatchTraction[i] =
1272  solidZoneTraction
1273  [
1274  structuresMesh.faceZones()[solidZoneID]
1275  .whichFace(solidPatchStart + i)
1276  ];
1277  }
1278 
1279  forAll(solidPatchPressure, i)
1280  {
1281  solidPatchPressure[i] =
1282  solidZonePressure
1283  [
1284  structuresMesh.faceZones()[solidZoneID]
1285  .whichFace(solidPatchStart + i)
1286  ];
1287  }
1288 
1289  if (fsi)
1290  {
1291  tForce.traction() = solidPatchTraction;
1292  tForce.pressure() = solidPatchPressure;
1293  }
1294 
1295  vector totalTractionForce =
1296  gSum
1297  (
1298  solidPatchTraction
1299  *structuresMesh.magSf().boundaryField()[solidPatchID]
1300  );
1301 
1302  Info << "FsiFoam:Step: Total traction force = "
1303  << totalTractionForce << endl;
1304  }
1305 
1306  // Info << "DU after setting interface force: " << DU << endl;
1307 
1308  //# include "solveSolid.H"
1309  if (solidDdtScheme == fv::EulerDdtScheme<vector>::typeName)
1310  {
1311  Info << "FsiFoam:Step: Solve Solid: Euler" << endl;
1312  //# include "solveSolidEuler.H"
1313  {
1314  const dictionary& stressControl =
1315  structuresMesh.solutionDict().subDict("solidMechanics");
1316 
1317  int nCorrStruct(readInt(stressControl.lookup("nCorrectors")));
1318  scalar convergenceTolerance(readScalar(stressControl.lookup("DU")));
1319  //# include "readSolidMechanicsControls.H"
1320 
1321  int iCorr = 0;
1322  lduMatrix::solverPerformance solverPerf;
1323  scalar initialResidual = 0;
1324 
1325  lduMatrix::debug = 0;
1326 
1327  //# include "EulerCoeffs.H"
1328  dimensionedScalar Cn("Cn", dimless/dimTime, 1.0/runTime.deltaT().value());
1329  dimensionedScalar Co("Co", dimless/dimTime, 1.0/runTime.deltaT().value());
1330 
1331  do
1332  {
1333  DU.storePrevIter();
1334 
1335  fvVectorMatrix DUEqn
1336  (
1337  Cn*rho*fvm::ddt(DU)
1338  - Co*rho*DV.oldTime()
1339  ==
1340  fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
1341  - fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
1342  + fvc::div
1343  (
1344  mu*gradDU.T()
1345  + lambda*(I*tr(gradDU))
1346  + mu*(gradDU&gradDU.T())
1347  + 0.5*lambda*(I*tr(gradDU & gradDU.T()))
1348  + (sigma & DF.T())
1349  + (DSigma & DF.T()),
1350  "div(sigma)"
1351  )
1352  );
1353 
1354  solverPerf = DUEqn.solve();
1355 
1356  DU.relax();
1357 
1358  if(iCorr == 0)
1359  {
1360  initialResidual = solverPerf.initialResidual();
1361  }
1362 
1363  gradDU = fvc::grad(DU);
1364 
1365  DF = gradDU.T();
1366 
1367  //# include "calculateDSigma.H"
1368  {
1369  volSymmTensorField Depsilon = symm(gradDU)
1370  + 0.5*symm(gradDU & gradDU.T());
1371 
1372  DSigma = 2*mu*Depsilon + I*(lambda*tr(Depsilon));
1373  }
1374  }
1375  while
1376  (
1377  solverPerf.initialResidual() > convergenceTolerance
1378  && ++iCorr < nCorrStruct
1379  );
1380 
1381  Info << "FsiFoam:Step: Solving for " << DU.name()
1382  << ", Initial residual = " << initialResidual
1383  << ", Final residual = " << solverPerf.initialResidual()
1384  << ", No outer iterations " << iCorr << endl;
1385 
1386  DV = fvc::ddt(DU);
1387 
1388  lduMatrix::debug = 1;
1389  }
1390  }
1391  else if (solidDdtScheme == fv::backwardDdtScheme<vector>::typeName)
1392  {
1393  Info << "FsiFoam:Step: Solve Solid: Backward" << endl;
1394  //# include "solveSolidBackward.H"
1395 
1396  {
1397  //# include "readSolidMechanicsControls.H"
1398  const dictionary& stressControl =
1399  structuresMesh.solutionDict().subDict("solidMechanics");
1400 
1401  int nCorrStruct(readInt(stressControl.lookup("nCorrectors")));
1402  scalar convergenceTolerance(readScalar(stressControl.lookup("DU")));
1403 
1404  int iCorr = 0;
1405  lduMatrix::solverPerformance solverPerf;
1406  scalar initialResidual = 0;
1407 
1408  lduMatrix::debug = 0;
1409 
1410  //# include "backwardCoeffs.H"
1411  dimensionedScalar Cn("Cn", dimless/dimTime, 0.0);
1412  dimensionedScalar Co("Co", dimless/dimTime, 0.0);
1413  dimensionedScalar Coo("Coo", dimless/dimTime, 0.0);
1414 
1415  scalar deltaT = runTime.deltaT().value();
1416  scalar deltaT0 = runTime.deltaT0().value();
1417 
1418  Cn.value() = 1 + deltaT/(deltaT + deltaT0);
1419  Coo.value() = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
1420  Co.value() = Cn.value() + Coo.value();
1421 
1422  if(runTime.timeIndex() == 1)
1423  {
1424  Cn.value() = 1.0;
1425  Co.value() = 1.0;
1426  Coo.value() = 0.0;
1427  }
1428 
1429  Cn.value() /= deltaT;
1430  Co.value() /= deltaT;
1431  Coo.value() /= deltaT;
1432 
1433  // Info << "ncorr = " << nCorr << endl;
1434  // Info << "convergence Tol: " << convergenceTolerance << endl;
1435  // Info << "DU:" << DU << endl;
1436  // Info << Cn << "," << Co << "," << Coo << endl;
1437  // Info << "RHO: " << rho << endl;
1438  // Info << "MU: " << mu << endl;
1439  // Info << "LAMBDA: " << lambda << endl;
1440  // Info << "gradDU:" << gradDU << endl;
1441  // Info << "gradDU.T(): " << gradDU.T() << endl;
1442  // Info << "DF: " << DF << endl;
1443  // Info << "DF.T(): " << DF.T() << endl;
1444  // Info << "SIGMA: " << sigma << endl;
1445  // Info << "DSigma: " << DSigma << endl;
1446  // Info << "DV: " << DV << endl;
1447  // Info << "DV.oldTime(): " << DV.oldTime() << endl;
1448  // Info << "DV.oldTime().oldTime(): " << DV.oldTime().oldTime() << endl;
1449 
1450  do
1451  {
1452  DU.storePrevIter();
1453 
1454  fvVectorMatrix DUEqn
1455  (
1456  Cn*rho*fvm::ddt(DU)
1457  - Co*rho*DV.oldTime()
1458  + Coo*rho*DV.oldTime().oldTime()
1459  ==
1460  fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
1461  - fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
1462  + fvc::div
1463  (
1464  mu*gradDU.T()
1465  + lambda*(I*tr(gradDU))
1466  + mu*(gradDU&gradDU.T())
1467  + 0.5*lambda*(I*tr(gradDU & gradDU.T()))
1468  + (sigma & DF.T())
1469  + (DSigma & DF.T()),
1470  "div(sigma)"
1471  )
1472  );
1473 
1474  solverPerf = DUEqn.solve();
1475 
1476  DU.relax();
1477 
1478  if(iCorr == 0)
1479  {
1480  initialResidual = solverPerf.initialResidual();
1481  }
1482 
1483  gradDU = fvc::grad(DU);
1484 
1485  DF = gradDU.T();
1486 
1487  //# include "calculateDSigma.H"
1488  {
1489  volSymmTensorField Depsilon = symm(gradDU)
1490  + 0.5*symm(gradDU & gradDU.T());
1491 
1492  DSigma = 2*mu*Depsilon + I*(lambda*tr(Depsilon));
1493  }
1494  // Info << "Solving for " << DU.name()
1495  // << ", Initial residual = " << initialResidual
1496  // << ", Final residual = " << solverPerf.initialResidual()
1497  // << ", No outer iterations " << iCorr << endl;
1498 
1499  }
1500  while
1501  (
1502  solverPerf.initialResidual() > convergenceTolerance
1503  && ++iCorr < nCorrStruct
1504  );
1505 
1506  Info << "FsiFoam:Step: Solving for " << DU.name()
1507  << ", Initial residual = " << initialResidual
1508  << ", Final residual = " << solverPerf.initialResidual()
1509  << ", No outer iterations " << iCorr << endl;
1510 
1511  DV = fvc::ddt(DU);
1512 
1513  lduMatrix::debug = 1;
1514  }
1515  }
1516  else
1517  {
1518  FatalErrorIn(args.executable())
1519  << "FsiFoam:Step: Wrong temporal (ddt) scheme for solid solver. "
1520  << "Possible schemes are: "
1521  << fv::EulerDdtScheme<vector>::typeName << " and "
1522  << fv::backwardDdtScheme<vector>::typeName
1523  << abort(FatalError);
1524  }
1525 
1526  //# include "calcFsiResidual.H"
1527  {
1528  const vectorField& solidPatchDisplacement =
1529  DU.boundaryField()[solidPatchID];
1530 
1531  vectorField solidZoneDisplacement
1532  (
1533  structuresMesh.faceZones()[solidZoneID]().size(),
1534  vector::zero
1535  );
1536 
1537  const label solidPatchStart =
1538  structuresMesh.boundaryMesh()[solidPatchID].start();
1539 
1540  forAll(solidPatchDisplacement, i)
1541  {
1542  solidZoneDisplacement
1543  [
1544  structuresMesh.faceZones()[solidZoneID]
1545  .whichFace(solidPatchStart + i)
1546  ] =
1547  solidPatchDisplacement[i];
1548  }
1549 
1550  // Parallel data exchange: collect displacement field on all processors
1551  reduce(solidZoneDisplacement, sumOp<vectorField>());
1552 
1553  vectorField fluidZoneDisplacement =
1554  interpolatorSolidFluid.faceInterpolate
1555  (
1556  solidZoneDisplacement
1557  );
1558 
1559  vectorField fluidPatchDisplacement
1560  (
1561  fluidsMesh.boundary()[fluidPatchID].size(),
1562  vector::zero
1563  );
1564 
1565  const label fluidPatchStart =
1566  fluidsMesh.boundaryMesh()[fluidPatchID].start();
1567 
1568  forAll(fluidPatchDisplacement, i)
1569  {
1570  fluidPatchDisplacement[i] =
1571  fluidZoneDisplacement
1572  [
1573  fluidsMesh.faceZones()[fluidZoneID].whichFace(fluidPatchStart + i)
1574  ];
1575  }
1576 
1577  primitivePatchInterpolation fluidPatchInterpolator
1578  (
1579  fluidsMesh.boundaryMesh()[fluidPatchID]
1580  );
1581 
1582  solidPatchPointsDispl =
1583  fluidPatchInterpolator.faceToPointInterpolate
1584  (
1585  fluidPatchDisplacement
1586  );
1587 
1588  fsiResidualOld = fsiResidual;
1589 
1590  fsiResidual = solidPatchPointsDispl - fluidPatchPointsDispl;
1591 
1592  // Masoud
1593  Info << "FsiFoam:Step: solidPatchPointsDispl = " << endl;
1594  Info << solidPatchPointsDispl << endl;
1595  Info << "FsiFoam:Step: fluidPatchPointsDispl = " << endl;
1596  Info << fluidPatchPointsDispl << endl;
1597  // Masoud : End
1598 
1599  // maxFsiResidual =
1600  // gMax
1601  // (
1602  // mag(fsiResidual)
1603  // /(mag(solidPatchPointsDispl) + SMALL)
1604  // );
1605 
1606  // Info << "Maximal fsi residual: " << maxFsiResidual << endl;
1607 
1608  fsiResidualNorm = ::sqrt(gSum(magSqr(fsiResidual)));
1609 
1610  if (outerCorr == 1)
1611  {
1612  initialFsiResidualNorm = fsiResidualNorm;
1613  }
1614 
1615  fsiResidualNorm /= initialFsiResidualNorm + SMALL;
1616 
1617  Info << "FsiFoam:Step: Current fsi residual norm: " << fsiResidualNorm << endl;
1618  }
1619  }
1620  while
1621  (
1622  (fsiResidualNorm > outerCorrTolerance)
1623  && (outerCorr < nOuterCorr)
1624  );
1625 
1626  Vs += DV;
1627 
1628  //# include "rotateSolidFields.H"
1629  {
1630  Info << "FsiFoam:Step: Rotating fields" << endl;
1631 
1632  F = I + DF;
1633 
1634  //U += DU;
1635  Usolid += DU;
1636 
1637  //epsilon += DEpsilon;
1638 
1639  sigma += DSigma;
1640 
1641  //volTensorField Finv = inv(F);
1642 
1643  J = det(F);
1644 
1645  rho = rho/J;
1646 
1647  //epsilon = symm(Finv.T() & epsilon & Finv);
1648 
1649  sigma = 1/J * symm(F & sigma & F.T());
1650  }
1651  //# include "moveSolidMeshLeastSquares.H"
1652  //--------------------------------------------------//
1653  //- move mesh
1654  //--------------------------------------------------//
1655  if(min(J.internalField()) > 0)
1656  {
1657  Info << "FsiFoam:Step: Moving mesh using least squares interpolation" << endl;
1658 
1659  leastSquaresVolPointInterpolation pointInterpolation(structuresMesh);
1660 
1661  // Create point mesh
1662  pointMesh pMesh(structuresMesh);
1663 
1664  wordList types
1665  (
1666  pMesh.boundary().size(),
1667  calculatedFvPatchVectorField::typeName
1668  );
1669 
1670  pointVectorField pointDU
1671  (
1672  IOobject
1673  (
1674  "pointDU",
1675  runTime.timeName(),
1676  structuresMesh
1677  ),
1678  pMesh,
1679  dimensionedVector("zero", dimLength, vector::zero),
1680  types
1681  );
1682 
1683  pointInterpolation.interpolate(DU, pointDU);
1684 
1685  const vectorField& pointDUI =
1686  pointDU.internalField();
1687 
1688  //- Move mesh
1689  vectorField newPoints = structuresMesh.allPoints();
1690 
1691  forAll (pointDUI, pointI)
1692  {
1693  newPoints[pointI] += pointDUI[pointI];
1694  }
1695 
1696  // Correct symmetryPlane points
1697 
1698  forAll(structuresMesh.boundaryMesh(), patchI)
1699  {
1700  if (isA<symmetryPolyPatch>(structuresMesh.boundaryMesh()[patchI]))
1701  {
1702  const labelList& meshPoints =
1703  structuresMesh.boundaryMesh()[patchI].meshPoints();
1704 
1705  vector avgN =
1706  gAverage(structuresMesh.boundaryMesh()[patchI].pointNormals());
1707 
1708  vector i(1, 0, 0);
1709  vector j(0, 1, 0);
1710  vector k(0, 0, 1);
1711 
1712  if (mag(avgN&i) > 0.95)
1713  {
1714  forAll(meshPoints, pI)
1715  {
1716  newPoints[meshPoints[pI]].x() = 0;
1717  }
1718  }
1719  else if (mag(avgN&j) > 0.95)
1720  {
1721  forAll(meshPoints, pI)
1722  {
1723  newPoints[meshPoints[pI]].y() = 0;
1724  }
1725  }
1726  else if (mag(avgN&k) > 0.95)
1727  {
1728  forAll(meshPoints, pI)
1729  {
1730  newPoints[meshPoints[pI]].z() = 0;
1731  }
1732  }
1733  }
1734  }
1735 
1736  //# include "calcUnusedNewPoints.H"
1737  forAll(globalFaceZones, zoneI)
1738  {
1739  const label curZoneID = globalFaceZones[zoneI];
1740 
1741  const labelList& curMap =
1742  globalToLocalFaceZonePointMap[zoneI];
1743 
1744  const labelList& curZoneMeshPoints =
1745  structuresMesh.faceZones()[curZoneID]().meshPoints();
1746 
1747  vectorField curGlobalZonePointDispl
1748  (
1749  curZoneMeshPoints.size(),
1750  vector::zero
1751  );
1752 
1753  //-Inter-proc points are shared by multiple procs
1754  // pointNumProc is the number of procs which a point lies on
1755  scalarField pointNumProcs(curZoneMeshPoints.size(), 0);
1756 
1757  forAll(curGlobalZonePointDispl, globalPointI)
1758  {
1759  label localPoint = curMap[globalPointI];
1760 
1761  if(curZoneMeshPoints[localPoint] < structuresMesh.nPoints())
1762  {
1763  label procPoint = curZoneMeshPoints[localPoint];
1764 
1765  curGlobalZonePointDispl[globalPointI] = pointDUI[procPoint];
1766 
1767  pointNumProcs[globalPointI] = 1;
1768  }
1769  }
1770 
1771  if (Pstream::parRun())
1772  {
1773  reduce(curGlobalZonePointDispl, sumOp<vectorField>());
1774  reduce(pointNumProcs, sumOp<scalarField>());
1775 
1776  //- now average the displacement between all procs
1777  curGlobalZonePointDispl /= pointNumProcs;
1778  }
1779 
1780  //- The curZonePointsDisplGlobal now contains the correct face zone
1781  // displacement in a global master processor order, now convert them
1782  // back into the local proc order
1783 
1784  vectorField curZonePointDispl
1785  (
1786  curZoneMeshPoints.size(),
1787  vector::zero
1788  );
1789 
1790  forAll(curGlobalZonePointDispl, globalPointI)
1791  {
1792  label localPoint = curMap[globalPointI];
1793 
1794  curZonePointDispl[localPoint] =
1795  curGlobalZonePointDispl[globalPointI];
1796  }
1797 
1798  forAll(curZonePointDispl, pointI)
1799  {
1800  // unused points
1801  if (curZoneMeshPoints[pointI] >= structuresMesh.nPoints())
1802  {
1803  newPoints[curZoneMeshPoints[pointI]] +=
1804  curZonePointDispl[pointI];
1805  }
1806  }
1807  }
1808 
1809  twoDPointCorrector twoDCorrector(structuresMesh);
1810  twoDCorrector.correctPoints(newPoints);
1811  structuresMesh.movePoints(newPoints);
1812  structuresMesh.V00();
1813  structuresMesh.moving(false);
1814  }
1815  else
1816  {
1817  FatalErrorIn(args.executable())
1818  << "FsiFoam:Step: Negative Jacobian"
1819  << exit(FatalError);
1820  }
1821  //# include "calculateStress.H"
1822  if (runTime.outputTime())
1823  {
1824  volScalarField sigmaEq
1825  (
1826  IOobject
1827  (
1828  "sigmaEq",
1829  runTime.timeName(),
1830  structuresMesh,
1831  IOobject::NO_READ,
1832  IOobject::AUTO_WRITE
1833  ),
1834  sqrt((3.0/2.0)*magSqr(dev(sigma)))
1835  );
1836 
1837  Info<< "FsiFoam:Step: Max sigmaEq = " << max(sigmaEq).value()
1838  << endl;
1839 
1840  runTime.write();
1841  }
1842 
1843  //# include "calculateLiftAndDrag.H"
1844 
1845  // Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
1846  // << " ClockTime = " << runTime.elapsedClockTime() << " s"
1847  // << endl << endl;
1848  }
1849  return 0;
1850 };
1851 
1852 int fsifoam_module::StepFluidAlone(){
1853 
1854  Foam::argList &args(this->ArgList());
1855  Foam::Time &runTime(this->RunTime());
1856  dynamicFvMesh &fluidsMesh(this->FluidMesh());
1857  // IOdictionary &transportProperties(this->TransportProperties());
1858  dimensionedScalar &nu(this->nu());
1859  dimensionedScalar &rhoFluid(this->rhoFluid());
1860  volScalarField &p(this->p());
1861  volVectorField &U(this->U());
1862  surfaceScalarField &phi(this->phi());
1863 
1864 
1865  // // IOdictionary &couplingProperties(this->CouplingProperties());
1866 
1867  label fluidPatchID(this->FluidPatchID());
1868  label fluidZoneID(this->FluidZoneID());
1869  label solidZoneID(this->SolidZoneID());
1870 
1871  bool feMotionSolver(this->FEMotion());
1872  bool fvMotionSolver(this->FVMotion());
1873  vectorField &accumulatedFluidInterfaceDisplacement(this->AccumulatedFluidInterfaceDisplacements());
1874  // Masoud : Checking AccumulatedFluidInterfaceDisplacement
1875  Info << "FsiFoam:StepFluidAlone: Checking AccumulatedFluidInterfaceDisplacementi = " << endl;
1876  Info << accumulatedFluidInterfaceDisplacement;
1877  // Masoud: End
1878  solidTractionFvPatchVectorField &tForce(this->tForce());
1879 
1880  zoneToZoneInterpolation &interpolatorFluidSolid(this->interpFluidSolid());
1881  zoneToZoneInterpolation &interpolatorSolidFluid(this->interpSolidFluid());
1882 
1883  scalar &sumLocalContErr(this->LocalContErr());
1884  scalar &globalContErr(this->GlobalContErr());
1885  scalar &cumulativeContErr(this->CumulativeContErr());
1886 
1887  labelList &globalFaceZones(this->GlobalFaceZones());
1888  labelListList &globalToLocalFaceZonePointMap(this->GlobalToLocalFaceZonePointMap());
1889  runTime++;
1890  if(!runTime.end()){
1891 
1892  this->ReadPISOControls();
1893  this->ReadFSIControls();
1894  word couplingScheme(this->CouplingScheme());
1895  // label &outerCorr(this->OuterCorr());
1896  scalar &fsiRelaxationFactor(this->FSIRelaxationFactor());
1897  // scalar &fsiRelaxationFactorMin(this->FSIRelaxationFactorMin());
1898  Switch fsi(this->FSIEnabled());
1899  int &nCorrPISO(this->NCorrPISO());
1900  int &nNonOrthCorr(this->NNonOrthCorr());
1901  scalar &interfaceDeformationLimit(this->InterfaceDeformationLimit());
1902  scalar &outerCorrTolerance(this->OuterCorrTolerance());
1903  int &nOuterCorr(this->NOuterCorr());
1904 
1905  vectorField fluidPatchPointsDispl
1906  (
1907  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
1908  vector::zero
1909  );
1910 
1911  vectorField fluidPatchPointsDisplOld
1912  (
1913  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
1914  vector::zero
1915  );
1916 
1917  vectorField solidPatchPointsDispl
1918  (
1919  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
1920  vector::zero
1921  );
1922 
1923  vectorField fsiResidual
1924  (
1925  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
1926  vector::zero
1927  );
1928 
1929  vectorField fsiResidualOld
1930  (
1931  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
1932  vector::zero
1933  );
1934 
1935 
1936  scalar initialFsiResidualNorm = 0;
1937  scalar fsiResidualNorm = 0;
1938  label outerCorr=0;
1939 
1940  do
1941  {
1942  outerCorr++;
1943  Info << "FsiFoam:StepFluidAlone: outerCorr = " << outerCorr << endl;
1944 
1945  //# include "setInterfaceDisplacement.H"
1946  Info << "\nFsiFoam:StepFluidAlone: Time = " << runTime.timeName()
1947  << ", iteration: " << outerCorr << endl;
1948 
1949  if (outerCorr < 3 || couplingScheme == "FixedRelaxation")
1950  {
1951  Info << "FsiFoam:StepFluidAlone: Current fsi under-relaxation factor: "
1952  << fsiRelaxationFactor << endl;
1953 
1954  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
1955 
1956  fluidPatchPointsDispl += fsiRelaxationFactor*fsiResidual;
1957  }
1958  else
1959  {
1960  if (couplingScheme == "Aitken")
1961  {
1962  fsiRelaxationFactor =
1963  -fsiRelaxationFactor
1964  *(
1965  gSum
1966  (
1967  fsiResidualOld
1968  &(fsiResidual - fsiResidualOld)
1969  )
1970  /(
1971  gSum
1972  (
1973  (fsiResidual - fsiResidualOld)
1974  &(fsiResidual - fsiResidualOld)
1975  )
1976  )
1977  );
1978 
1979  fsiRelaxationFactor = mag(fsiRelaxationFactor);
1980 
1981  Info << "FsiFoam:StepFluidAlone: Current fsi under-relaxation factor (Aitken): "
1982  << fsiRelaxationFactor << endl;
1983 
1984  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
1985 
1986  fluidPatchPointsDispl +=
1987  fsiRelaxationFactor*fsiResidual;
1988  }
1989  }
1990 
1991  {
1992 
1993  const vectorField& n =
1994  fluidsMesh.boundaryMesh()[fluidPatchID].pointNormals();
1995 
1996  primitivePatchInterpolation patchInterpolator
1997  (
1998  fluidsMesh.boundaryMesh()[fluidPatchID]
1999  );
2000 
2001  scalarField pointDeltaCoeffs =
2002  patchInterpolator.faceToPointInterpolate
2003  (
2004  fluidsMesh.boundary()[fluidPatchID].deltaCoeffs()
2005  );
2006 
2007  scalar delta =
2008  gMax
2009  (
2010  mag
2011  (
2012  n
2013  & (
2014  accumulatedFluidInterfaceDisplacement
2015  + fluidPatchPointsDispl
2016  - fluidPatchPointsDisplOld
2017  )
2018  )
2019  *pointDeltaCoeffs
2020  );
2021 
2022  Info << "FsiFoam:StepFluidAlone: Maximal accumulated displacement of interface points: "
2023  << delta << endl;
2024 
2025  if(delta < interfaceDeformationLimit)
2026  {
2027  // Move only interface points
2028  // Masoud
2029  Info << "FsiFoam:StepFluidAlone: Moving only interface...";
2030  // Masoud : End
2031  pointField newPoints = fluidsMesh.allPoints();
2032 
2033  const labelList& meshPoints =
2034  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2035 
2036  forAll (fluidPatchPointsDispl, pointI)
2037  {
2038  //Info << "1. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2039  //<< endl;
2040 
2041  newPoints[meshPoints[pointI]] +=
2042  fluidPatchPointsDispl[pointI]
2043  - fluidPatchPointsDisplOld[pointI];
2044  }
2045 
2046  twoDPointCorrector twoDCorrector(fluidsMesh);
2047 
2048  twoDCorrector.correctPoints(newPoints);
2049 
2050  fluidsMesh.movePoints(newPoints);
2051 
2052  // Accumulate interface points displacement
2053  accumulatedFluidInterfaceDisplacement +=
2054  fluidPatchPointsDispl
2055  - fluidPatchPointsDisplOld;
2056  }
2057  else
2058  {
2059  // Move whole fluid mesh
2060  // Masoud
2061  Info << "FsiFoam:StepFluidAlone: Moving the whole mesh...";
2062  // Masoud: End
2063  pointField newPoints = fluidsMesh.allPoints();
2064 
2065  const labelList& meshPoints =
2066  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2067 
2068  forAll (accumulatedFluidInterfaceDisplacement, pointI)
2069  {
2070  //Info << "2. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2071  //<< endl;
2072  newPoints[meshPoints[pointI]] -=
2073  accumulatedFluidInterfaceDisplacement[pointI];
2074  }
2075 
2076  twoDPointCorrector twoDCorrector(fluidsMesh);
2077 
2078  twoDCorrector.correctPoints(newPoints);
2079 
2080  fluidsMesh.movePoints(newPoints);
2081 
2082  accumulatedFluidInterfaceDisplacement +=
2083  fluidPatchPointsDispl
2084  - fluidPatchPointsDisplOld;
2085 
2086 
2087  if (feMotionSolver)
2088  {
2089  tetPointVectorField& motionU =
2090  const_cast<tetPointVectorField&>
2091  (
2092  fluidsMesh.objectRegistry:: lookupObject<tetPointVectorField>
2093  (
2094  "motionU"
2095  )
2096  );
2097 
2098  fixedValueTetPolyPatchVectorField& motionUFluidPatch =
2099  refCast<fixedValueTetPolyPatchVectorField>
2100  (
2101  motionU.boundaryField()[fluidPatchID]
2102  );
2103 
2104  tetPolyPatchInterpolation tppi
2105  (
2106  refCast<const faceTetPolyPatch>(motionUFluidPatch.patch())
2107  );
2108 
2109  motionUFluidPatch ==
2110  tppi.pointToPointInterpolate
2111  (
2112  accumulatedFluidInterfaceDisplacement
2113  /runTime.deltaT().value()
2114  );
2115  }
2116  else if (fvMotionSolver)
2117  {
2118  pointVectorField& motionU =
2119  const_cast<pointVectorField&>
2120  (
2121  fluidsMesh.objectRegistry:: lookupObject<pointVectorField>
2122  (
2123  "pointMotionU"
2124  )
2125  );
2126 
2127  fixedValuePointPatchVectorField& motionUFluidPatch =
2128  refCast<fixedValuePointPatchVectorField>
2129  (
2130  motionU.boundaryField()[fluidPatchID]
2131  );
2132 
2133  motionUFluidPatch ==
2134  accumulatedFluidInterfaceDisplacement
2135  /runTime.deltaT().value();
2136  }
2137  else
2138  {
2139  FatalErrorIn(args.executable())
2140  << "FsiFoam:StepFluidAlone: Problem with mesh motion solver selection"
2141  << abort(FatalError);
2142  }
2143 
2144  fluidsMesh.update();
2145 
2146  accumulatedFluidInterfaceDisplacement = vector::zero;
2147  }
2148  }
2149 
2150  if(fluidsMesh.moving())
2151  {
2152  // Make the fluxes relative
2153  phi -= fvc::meshPhi(U);
2154  }
2155 
2156 
2157  scalar CoNum = 0.0;
2158  scalar meanCoNum = 0.0;
2159  scalar velMag = 0.0;
2160 
2161  if (fluidsMesh.nInternalFaces())
2162  {
2163  surfaceScalarField magPhi = mag(phi);
2164 
2165  surfaceScalarField SfUfbyDelta =
2166  fluidsMesh.surfaceInterpolation::deltaCoeffs()*magPhi;
2167 
2168  const scalar deltaT = runTime.deltaT().value();
2169 
2170  CoNum = max(SfUfbyDelta/fluidsMesh.magSf()).value()*deltaT;
2171 
2172  meanCoNum = (sum(SfUfbyDelta)/sum(fluidsMesh.magSf())).value()*deltaT;
2173 
2174  velMag = max(magPhi/fluidsMesh.magSf()).value();
2175  }
2176 
2177  Info<< "FsiFoam:StepFluidAlone: Courant Number mean: " << meanCoNum
2178  << " max: " << CoNum
2179  << " velocity magnitude: " << velMag
2180  << endl;
2181 
2182  fvVectorMatrix UEqn
2183  (
2184  fvm::ddt(U)
2185  + fvm::div(phi, U)
2186  - fvm::laplacian(nu, U)
2187  );
2188 
2189  solve(UEqn == -fvc::grad(p));
2190 
2191  // --- PISO loop
2192  volScalarField rUA = 1.0/UEqn.A();
2193  // Masoud
2194  Info << "FsiFoam:StepFluidAlone: Performing " << nCorrPISO << " Pressure corrections." << endl;
2195  // Masoud End
2196 
2197  for (int corr=0; corr<nCorrPISO; corr++)
2198  {
2199  U = rUA*UEqn.H();
2200  phi = (fvc::interpolate(U) & fluidsMesh.Sf());
2201 
2202  adjustPhi(phi, U, p);
2203 
2204  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
2205  {
2206  fvScalarMatrix pEqn
2207  (
2208  fvm::laplacian(rUA, p)
2209  == fvc::div(phi)
2210  );
2211 
2212  pEqn.setReference(pRefCell, pRefValue);
2213  pEqn.solve();
2214 
2215  // Added by Masoud
2216  //Info << "p = " << (this-> p()) << endl;
2217 
2218  if (nonOrth == nNonOrthCorr)
2219  {
2220  phi -= pEqn.flux();
2221  }
2222  }
2223 
2224  {
2225  volScalarField contErr = fvc::div(phi);
2226 
2227  sumLocalContErr = runTime.deltaT().value()*
2228  mag(contErr)().weightedAverage(fluidsMesh.V()).value();
2229 
2230  globalContErr = runTime.deltaT().value()*
2231  contErr.weightedAverage(fluidsMesh.V()).value();
2232 
2233  cumulativeContErr += globalContErr;
2234 
2235  Info<< "FsiFoam:StepFluidAlone: time step continuity errors : sum local = " << sumLocalContErr
2236  << ", global = " << globalContErr
2237  << ", cumulative = " << cumulativeContErr
2238  << endl;
2239  }
2240 
2241  U -= rUA*fvc::grad(p);
2242  U.correctBoundaryConditions();
2243  }
2244 
2245  {
2246  Info << "FsiFoam:StepFluidAlone: Not Setting traction on solid patch" << endl;
2247 
2248  {
2249 
2250  fsiResidualOld = fsiResidual;
2251  // Original: a fixed displacement delta will be used for the solid
2252  //this->UpdateFSISurface(solidPatchPointsDispl);
2253  // Oridinal : End
2254 
2255  // Masoud : in the initial step displacement delta will be added
2256  // and for the rest of the process, it will be consider-
2257  // as zeros vector.
2258  if (outerCorr == 1)
2259  {
2260  this->UpdateFSISurface(solidPatchPointsDispl);
2261  } else {
2262  solidPatchPointsDispl = vector::zero;
2263  //fluidPatchPointsDispl = vector::zero;
2264  }
2265  // Masoud : End
2266 
2267  // Masoud
2268  Info << "FsiFoam:StepFluidAlone: solidPatchPointsDispl = " << solidPatchPointsDispl << endl;
2269  //Masoud end
2270 
2271  // Original (solidPatchPointsDispl won't change by iteration
2272  // therefore huge pressure build-up happens)
2273  fsiResidual = (solidPatchPointsDispl - fluidPatchPointsDispl);
2274  // Original end
2275 
2276 
2277  fsiResidualNorm = ::sqrt(gSum(magSqr(fsiResidual)));
2278 
2279  if (outerCorr == 1)
2280  {
2281  initialFsiResidualNorm = fsiResidualNorm;
2282  }
2283 
2284  fsiResidualNorm /= initialFsiResidualNorm + SMALL;
2285 
2286  Info << "FsiFoam:StepFluidAlone: Current fsi residual norm: " << fsiResidualNorm << endl;
2287  }
2288  }
2289  }
2290  while
2291  (
2292  (fsiResidualNorm > outerCorrTolerance)
2293  && (outerCorr < nOuterCorr) //Masoud changed from outerCorr < nOuterCorr
2294  );
2295 
2296 
2297  if (runTime.outputTime())
2298  {
2299  runTime.write();
2300  }
2301  }
2302  return(0);
2303 }
2304 
2305 // Masoud: This is simple iteration-less scheme
2306 int fsifoam_module::StepFluidNonItr(){
2307 
2308  Foam::argList &args(this->ArgList());
2309  Foam::Time &runTime(this->RunTime());
2310  dynamicFvMesh &fluidsMesh(this->FluidMesh());
2311  // IOdictionary &transportProperties(this->TransportProperties());
2312  dimensionedScalar &nu(this->nu());
2313  dimensionedScalar &rhoFluid(this->rhoFluid());
2314  volScalarField &p(this->p());
2315  volVectorField &U(this->U());
2316  surfaceScalarField &phi(this->phi());
2317 
2318 
2319  // // IOdictionary &couplingProperties(this->CouplingProperties());
2320 
2321  label fluidPatchID(this->FluidPatchID());
2322  label fluidZoneID(this->FluidZoneID());
2323  label solidZoneID(this->SolidZoneID());
2324 
2325  bool feMotionSolver(this->FEMotion());
2326  bool fvMotionSolver(this->FVMotion());
2327  vectorField &accumulatedFluidInterfaceDisplacement(this->AccumulatedFluidInterfaceDisplacements());
2328  //Info << "Checking AccumulatedFluidInterfaceDisplacementi = " << endl;
2329  //Info << accumulatedFluidInterfaceDisplacement;
2330  solidTractionFvPatchVectorField &tForce(this->tForce());
2331 
2332  zoneToZoneInterpolation &interpolatorFluidSolid(this->interpFluidSolid());
2333  zoneToZoneInterpolation &interpolatorSolidFluid(this->interpSolidFluid());
2334 
2335  scalar &sumLocalContErr(this->LocalContErr());
2336  scalar &globalContErr(this->GlobalContErr());
2337  scalar &cumulativeContErr(this->CumulativeContErr());
2338 
2339  labelList &globalFaceZones(this->GlobalFaceZones());
2340  labelListList &globalToLocalFaceZonePointMap(this->GlobalToLocalFaceZonePointMap());
2341  runTime++;
2342  if(!runTime.end()){
2343 
2344  this->ReadPISOControls();
2345  this->ReadFSIControls();
2346  word couplingScheme(this->CouplingScheme());
2347  // label &outerCorr(this->OuterCorr());
2348  scalar &fsiRelaxationFactor(this->FSIRelaxationFactor());
2349  // scalar &fsiRelaxationFactorMin(this->FSIRelaxationFactorMin());
2350  Switch fsi(this->FSIEnabled());
2351  int &nCorrPISO(this->NCorrPISO());
2352  int &nNonOrthCorr(this->NNonOrthCorr());
2353  scalar &interfaceDeformationLimit(this->InterfaceDeformationLimit());
2354  scalar &outerCorrTolerance(this->OuterCorrTolerance());
2355  int &nOuterCorr(this->NOuterCorr());
2356 
2357  vectorField fluidPatchPointsDispl
2358  (
2359  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2360  vector::zero
2361  );
2362 
2363  vectorField fluidPatchPointsDisplOld
2364  (
2365  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2366  vector::zero
2367  );
2368 
2369  vectorField solidPatchPointsDispl
2370  (
2371  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2372  vector::zero
2373  );
2374 
2375  vectorField fsiResidual
2376  (
2377  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2378  vector::zero
2379  );
2380 
2381  vectorField fsiResidualOld
2382  (
2383  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2384  vector::zero
2385  );
2386 
2387 
2388  scalar initialFsiResidualNorm = 0;
2389  scalar fsiResidualNorm = 0;
2390  label outerCorr=0;
2391 
2392  do
2393  {
2394  outerCorr++;
2395  Info << "FsiFoam:StepFluidNonItr: outerCorr = " << outerCorr << endl;
2396 
2397  //# include "setInterfaceDisplacement.H"
2398  Info << "FsiFoam:StepFluidNonItr: Time = " << runTime.timeName()
2399  << ", iteration: " << outerCorr << endl;
2400  // Updating fluid patch displacement
2401  this->UpdateFSISurface(fluidPatchPointsDispl);
2402  //Info << " Updating fluid interface coordinates with : " << endl;
2403  //Info << fluidPatchPointsDispl << endl;
2404  {
2405 
2406  const vectorField& n =
2407  fluidsMesh.boundaryMesh()[fluidPatchID].pointNormals();
2408 
2409  primitivePatchInterpolation patchInterpolator
2410  (
2411  fluidsMesh.boundaryMesh()[fluidPatchID]
2412  );
2413 
2414  scalarField pointDeltaCoeffs =
2415  patchInterpolator.faceToPointInterpolate
2416  (
2417  fluidsMesh.boundary()[fluidPatchID].deltaCoeffs()
2418  );
2419 
2420  scalar delta =
2421  gMax
2422  (
2423  mag
2424  (
2425  n
2426  & (
2427  accumulatedFluidInterfaceDisplacement
2428  + fluidPatchPointsDispl
2429  )
2430  )
2431  *pointDeltaCoeffs
2432  );
2433 
2434  Info << "FsiFoam:StepFluidNonItr: Maximal accumulated displacement of interface points: "
2435  << delta << endl;
2436 
2437  if(delta < interfaceDeformationLimit)
2438  {
2439  // Move only interface points
2440  Info << "FsiFoam:StepFluidNonItr: Moving only interface...";
2441  pointField newPoints = fluidsMesh.allPoints();
2442 
2443  const labelList& meshPoints =
2444  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2445 
2446  forAll (fluidPatchPointsDispl, pointI)
2447  {
2448  //Info << "1. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2449  //<< endl;
2450 
2451  newPoints[meshPoints[pointI]] +=
2452  fluidPatchPointsDispl[pointI];
2453  }
2454 
2455  twoDPointCorrector twoDCorrector(fluidsMesh);
2456 
2457  twoDCorrector.correctPoints(newPoints);
2458 
2459  fluidsMesh.movePoints(newPoints);
2460 
2461  // Accumulate interface points displacement
2462  accumulatedFluidInterfaceDisplacement +=
2463  fluidPatchPointsDispl;
2464  }
2465  else
2466  {
2467  // Move whole fluid mesh
2468  Info << "FsiFoam:StepFluidNonItr: Moving the whole mesh...";
2469  pointField newPoints = fluidsMesh.allPoints();
2470 
2471  const labelList& meshPoints =
2472  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2473 
2474  forAll (accumulatedFluidInterfaceDisplacement, pointI)
2475  {
2476  //Info << "2. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2477  //<< endl;
2478  newPoints[meshPoints[pointI]] -=
2479  accumulatedFluidInterfaceDisplacement[pointI];
2480  }
2481 
2482  twoDPointCorrector twoDCorrector(fluidsMesh);
2483 
2484  twoDCorrector.correctPoints(newPoints);
2485 
2486  fluidsMesh.movePoints(newPoints);
2487 
2488  accumulatedFluidInterfaceDisplacement +=
2489  fluidPatchPointsDispl;
2490 
2491  if (feMotionSolver)
2492  {
2493  tetPointVectorField& motionU =
2494  const_cast<tetPointVectorField&>
2495  (
2496  fluidsMesh.objectRegistry:: lookupObject<tetPointVectorField>
2497  (
2498  "motionU"
2499  )
2500  );
2501 
2502  fixedValueTetPolyPatchVectorField& motionUFluidPatch =
2503  refCast<fixedValueTetPolyPatchVectorField>
2504  (
2505  motionU.boundaryField()[fluidPatchID]
2506  );
2507 
2508  tetPolyPatchInterpolation tppi
2509  (
2510  refCast<const faceTetPolyPatch>(motionUFluidPatch.patch())
2511  );
2512 
2513  motionUFluidPatch ==
2514  tppi.pointToPointInterpolate
2515  (
2516  accumulatedFluidInterfaceDisplacement
2517  /runTime.deltaT().value()
2518  );
2519  }
2520  else if (fvMotionSolver)
2521  {
2522  pointVectorField& motionU =
2523  const_cast<pointVectorField&>
2524  (
2525  fluidsMesh.objectRegistry:: lookupObject<pointVectorField>
2526  (
2527  "pointMotionU"
2528  )
2529  );
2530 
2531  fixedValuePointPatchVectorField& motionUFluidPatch =
2532  refCast<fixedValuePointPatchVectorField>
2533  (
2534  motionU.boundaryField()[fluidPatchID]
2535  );
2536 
2537  motionUFluidPatch ==
2538  accumulatedFluidInterfaceDisplacement
2539  /runTime.deltaT().value();
2540  }
2541  else
2542  {
2543  FatalErrorIn(args.executable())
2544  << "FsiFoam:StepFluidNonItr: Problem with mesh motion solver selection"
2545  << abort(FatalError);
2546  }
2547 
2548 
2549  Info << "FsiFoam:StepFluidNonItr: runTime.deltaT() = "
2550  << runTime.deltaT().value() << endl;
2551 
2552 
2553  fluidsMesh.update();
2554 
2555  accumulatedFluidInterfaceDisplacement = vector::zero;
2556  }
2557  }
2558 
2559  if(fluidsMesh.moving())
2560  {
2561  // Make the fluxes relative
2562  Info << "FsiFoam:StepFluidNonItr: Compensating for mesh movment !" << endl;
2563  phi -= fvc::meshPhi(U);
2564  }
2565 
2566 
2567  scalar CoNum = 0.0;
2568  scalar meanCoNum = 0.0;
2569  scalar velMag = 0.0;
2570 
2571  if (fluidsMesh.nInternalFaces())
2572  {
2573  surfaceScalarField magPhi = mag(phi);
2574 
2575  surfaceScalarField SfUfbyDelta =
2576  fluidsMesh.surfaceInterpolation::deltaCoeffs()*magPhi;
2577 
2578  const scalar deltaT = runTime.deltaT().value();
2579 
2580  CoNum = max(SfUfbyDelta/fluidsMesh.magSf()).value()*deltaT;
2581 
2582  meanCoNum = (sum(SfUfbyDelta)/sum(fluidsMesh.magSf())).value()*deltaT;
2583 
2584  velMag = max(magPhi/fluidsMesh.magSf()).value();
2585  }
2586 
2587  Info<< "FsiFoam:StepFluidNonItr: Courant Number mean: " << meanCoNum
2588  << " max: " << CoNum
2589  << " velocity magnitude: " << velMag
2590  << endl;
2591 
2592  fvVectorMatrix UEqn
2593  (
2594  fvm::ddt(U)
2595  + fvm::div(phi, U)
2596  - fvm::laplacian(nu, U)
2597  );
2598 
2599  solve(UEqn == -fvc::grad(p));
2600 
2601  // --- PISO loop
2602  volScalarField rUA = 1.0/UEqn.A();
2603  //Info << "nCorrPISO = " << nCorrPISO << endl;
2604  //Info << "nNonOrthCorr = " << nNonOrthCorr << endl;
2605  for (int corr=0; corr<nCorrPISO; corr++)
2606  {
2607  U = rUA*UEqn.H();
2608  phi = (fvc::interpolate(U) & fluidsMesh.Sf());
2609 
2610  adjustPhi(phi, U, p);
2611 
2612  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
2613  {
2614  fvScalarMatrix pEqn
2615  (
2616  fvm::laplacian(rUA, p)
2617  == fvc::div(phi)
2618  );
2619 
2620  pEqn.setReference(pRefCell, pRefValue);
2621  pEqn.solve();
2622 
2623  //Info << "p = " << p << endl;
2624 
2625  if (nonOrth == nNonOrthCorr)
2626  {
2627  phi -= pEqn.flux();
2628  }
2629  }
2630 
2631  {
2632  volScalarField contErr = fvc::div(phi);
2633 
2634  sumLocalContErr = runTime.deltaT().value()*
2635  mag(contErr)().weightedAverage(fluidsMesh.V()).value();
2636 
2637  globalContErr = runTime.deltaT().value()*
2638  contErr.weightedAverage(fluidsMesh.V()).value();
2639 
2640  cumulativeContErr += globalContErr;
2641 
2642  Info<< "FsiFoam:StepFluidNonItr: time step continuity errors : sum local = " << sumLocalContErr
2643  << ", global = " << globalContErr
2644  << ", cumulative = " << cumulativeContErr
2645  << endl;
2646  }
2647 
2648  U -= rUA*fvc::grad(p);
2649  U.correctBoundaryConditions();
2650  }
2651 
2652  }
2653  while
2654  ( outerCorr < 1 );
2655 
2656 
2657  if (runTime.outputTime())
2658  {
2659  runTime.write();
2660  }
2661  }
2662  return(0);
2663 }
2664 //Masoud End (StepFluidNonItr)
2665 
2666 
2667 //Masoud: StepFluidItr
2668 int fsifoam_module::StepFluidItr(){
2669 
2670  Foam::argList &args(this->ArgList());
2671  Foam::Time &runTime(this->RunTime());
2672  dynamicFvMesh &fluidsMesh(this->FluidMesh());
2673  // IOdictionary &transportProperties(this->TransportProperties());
2674  dimensionedScalar &nu(this->nu());
2675  dimensionedScalar &rhoFluid(this->rhoFluid());
2676  volScalarField &p(this->p());
2677  volVectorField &U(this->U());
2678  surfaceScalarField &phi(this->phi());
2679 
2680 
2681  // // IOdictionary &couplingProperties(this->CouplingProperties());
2682 
2683  label fluidPatchID(this->FluidPatchID());
2684  label fluidZoneID(this->FluidZoneID());
2685  label solidZoneID(this->SolidZoneID());
2686 
2687  bool feMotionSolver(this->FEMotion());
2688  bool fvMotionSolver(this->FVMotion());
2689  vectorField &accumulatedFluidInterfaceDisplacement(this->AccumulatedFluidInterfaceDisplacements());
2690  // Masoud : Checking AccumulatedFluidInterfaceDisplacement
2691  Info << "FsiFoam:StepFluidItr: Checking AccumulatedFluidInterfaceDisplacementi = " << endl;
2692  Info << accumulatedFluidInterfaceDisplacement;
2693  // Masoud: End
2694  solidTractionFvPatchVectorField &tForce(this->tForce());
2695 
2696  zoneToZoneInterpolation &interpolatorFluidSolid(this->interpFluidSolid());
2697  zoneToZoneInterpolation &interpolatorSolidFluid(this->interpSolidFluid());
2698 
2699  scalar &sumLocalContErr(this->LocalContErr());
2700  scalar &globalContErr(this->GlobalContErr());
2701  scalar &cumulativeContErr(this->CumulativeContErr());
2702 
2703  labelList &globalFaceZones(this->GlobalFaceZones());
2704  labelListList &globalToLocalFaceZonePointMap(this->GlobalToLocalFaceZonePointMap());
2705  runTime++;
2706  if(!runTime.end()){
2707 
2708  this->ReadPISOControls();
2709  this->ReadFSIControls();
2710  word couplingScheme(this->CouplingScheme());
2711  // label &outerCorr(this->OuterCorr());
2712  scalar &fsiRelaxationFactor(this->FSIRelaxationFactor());
2713  // scalar &fsiRelaxationFactorMin(this->FSIRelaxationFactorMin());
2714  Switch fsi(this->FSIEnabled());
2715  int &nCorrPISO(this->NCorrPISO());
2716  int &nNonOrthCorr(this->NNonOrthCorr());
2717  scalar &interfaceDeformationLimit(this->InterfaceDeformationLimit());
2718  scalar &outerCorrTolerance(this->OuterCorrTolerance());
2719  int &nOuterCorr(this->NOuterCorr());
2720 
2721  vectorField fluidPatchPointsDispl
2722  (
2723  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2724  vector::zero
2725  );
2726 
2727  vectorField fluidPatchPointsDisplOld
2728  (
2729  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2730  vector::zero
2731  );
2732 
2733  vectorField solidPatchPointsDispl
2734  (
2735  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2736  vector::zero
2737  );
2738 
2739  vectorField fsiResidual
2740  (
2741  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2742  vector::zero
2743  );
2744 
2745  vectorField fsiResidualOld
2746  (
2747  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2748  vector::zero
2749  );
2750 
2751 
2752  scalar initialFsiResidualNorm = 0;
2753  scalar fsiResidualNorm = 0;
2754  label outerCorr=0;
2755 
2756  do
2757  {
2758  outerCorr++;
2759  Info << "FsiFoam:StepFluidItr: outerCorr = " << outerCorr << endl;
2760 
2761  //# include "setInterfaceDisplacement.H"
2762  Info << "\nFsiFoam:StepFluidItr: Time = " << runTime.timeName()
2763  << ", iteration: " << outerCorr << endl;
2764 
2765  if (outerCorr < 3 || couplingScheme == "FixedRelaxation")
2766  {
2767  Info << "FsiFoam:StepFluidItr: Current fsi under-relaxation factor: "
2768  << fsiRelaxationFactor << endl;
2769 
2770  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
2771 
2772  fluidPatchPointsDispl += fsiRelaxationFactor*fsiResidual;
2773  }
2774  else
2775  {
2776  if (couplingScheme == "Aitken")
2777  {
2778  fsiRelaxationFactor =
2779  -fsiRelaxationFactor
2780  *(
2781  gSum
2782  (
2783  fsiResidualOld
2784  &(fsiResidual - fsiResidualOld)
2785  )
2786  /(
2787  gSum
2788  (
2789  (fsiResidual - fsiResidualOld)
2790  &(fsiResidual - fsiResidualOld)
2791  )
2792  )
2793  );
2794 
2795  fsiRelaxationFactor = mag(fsiRelaxationFactor);
2796 
2797  Info << "FsiFoam:StepFluidItr: Current fsi under-relaxation factor (Aitken): "
2798  << fsiRelaxationFactor << endl;
2799 
2800  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
2801 
2802  fluidPatchPointsDispl +=
2803  fsiRelaxationFactor*fsiResidual;
2804  }
2805  }
2806 
2807  {
2808 
2809  const vectorField& n =
2810  fluidsMesh.boundaryMesh()[fluidPatchID].pointNormals();
2811 
2812  primitivePatchInterpolation patchInterpolator
2813  (
2814  fluidsMesh.boundaryMesh()[fluidPatchID]
2815  );
2816 
2817  scalarField pointDeltaCoeffs =
2818  patchInterpolator.faceToPointInterpolate
2819  (
2820  fluidsMesh.boundary()[fluidPatchID].deltaCoeffs()
2821  );
2822 
2823  scalar delta =
2824  gMax
2825  (
2826  mag
2827  (
2828  n
2829  & (
2830  accumulatedFluidInterfaceDisplacement
2831  + fluidPatchPointsDispl
2832  - fluidPatchPointsDisplOld
2833  )
2834  )
2835  *pointDeltaCoeffs
2836  );
2837 
2838  Info << "FsiFoam:StepFluidItr: Maximal accumulated displacement of interface points: "
2839  << delta << endl;
2840 
2841  if(delta < interfaceDeformationLimit)
2842  {
2843  // Move only interface points
2844  // Masoud
2845  Info << "FsiFoam:StepFluidItr: Moving only interface...";
2846  // Masoud : End
2847  pointField newPoints = fluidsMesh.allPoints();
2848 
2849  const labelList& meshPoints =
2850  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2851 
2852  forAll (fluidPatchPointsDispl, pointI)
2853  {
2854  //Info << "1. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2855  //<< endl;
2856 
2857  newPoints[meshPoints[pointI]] +=
2858  fluidPatchPointsDispl[pointI]
2859  - fluidPatchPointsDisplOld[pointI];
2860  }
2861 
2862  twoDPointCorrector twoDCorrector(fluidsMesh);
2863 
2864  twoDCorrector.correctPoints(newPoints);
2865 
2866  fluidsMesh.movePoints(newPoints);
2867 
2868  // Accumulate interface points displacement
2869  accumulatedFluidInterfaceDisplacement +=
2870  fluidPatchPointsDispl
2871  - fluidPatchPointsDisplOld;
2872  }
2873  else
2874  {
2875  // Move whole fluid mesh
2876  // Masoud
2877  Info << "FsiFoam:StepFluidItr: Moving the whole mesh...";
2878  // Masoud: End
2879  pointField newPoints = fluidsMesh.allPoints();
2880 
2881  const labelList& meshPoints =
2882  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2883 
2884  forAll (accumulatedFluidInterfaceDisplacement, pointI)
2885  {
2886  //Info << "2. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2887  //<< endl;
2888  newPoints[meshPoints[pointI]] -=
2889  accumulatedFluidInterfaceDisplacement[pointI];
2890  }
2891 
2892  twoDPointCorrector twoDCorrector(fluidsMesh);
2893 
2894  twoDCorrector.correctPoints(newPoints);
2895 
2896  fluidsMesh.movePoints(newPoints);
2897 
2898  accumulatedFluidInterfaceDisplacement +=
2899  fluidPatchPointsDispl
2900  - fluidPatchPointsDisplOld;
2901 
2902 
2903  if (feMotionSolver)
2904  {
2905  tetPointVectorField& motionU =
2906  const_cast<tetPointVectorField&>
2907  (
2908  fluidsMesh.objectRegistry:: lookupObject<tetPointVectorField>
2909  (
2910  "motionU"
2911  )
2912  );
2913 
2914  fixedValueTetPolyPatchVectorField& motionUFluidPatch =
2915  refCast<fixedValueTetPolyPatchVectorField>
2916  (
2917  motionU.boundaryField()[fluidPatchID]
2918  );
2919 
2920  tetPolyPatchInterpolation tppi
2921  (
2922  refCast<const faceTetPolyPatch>(motionUFluidPatch.patch())
2923  );
2924 
2925  motionUFluidPatch ==
2926  tppi.pointToPointInterpolate
2927  (
2928  accumulatedFluidInterfaceDisplacement
2929  /runTime.deltaT().value()
2930  );
2931  }
2932  else if (fvMotionSolver)
2933  {
2934  pointVectorField& motionU =
2935  const_cast<pointVectorField&>
2936  (
2937  fluidsMesh.objectRegistry:: lookupObject<pointVectorField>
2938  (
2939  "pointMotionU"
2940  )
2941  );
2942 
2943  fixedValuePointPatchVectorField& motionUFluidPatch =
2944  refCast<fixedValuePointPatchVectorField>
2945  (
2946  motionU.boundaryField()[fluidPatchID]
2947  );
2948 
2949  motionUFluidPatch ==
2950  accumulatedFluidInterfaceDisplacement
2951  /runTime.deltaT().value();
2952  }
2953  else
2954  {
2955  FatalErrorIn(args.executable())
2956  << "FsiFoam:StepFluidItr: Problem with mesh motion solver selection"
2957  << abort(FatalError);
2958  }
2959 
2960  fluidsMesh.update();
2961 
2962  accumulatedFluidInterfaceDisplacement = vector::zero;
2963  }
2964  }
2965 
2966  if(fluidsMesh.moving())
2967  {
2968  // Make the fluxes relative
2969  phi -= fvc::meshPhi(U);
2970  }
2971 
2972 
2973  scalar CoNum = 0.0;
2974  scalar meanCoNum = 0.0;
2975  scalar velMag = 0.0;
2976 
2977  if (fluidsMesh.nInternalFaces())
2978  {
2979  surfaceScalarField magPhi = mag(phi);
2980 
2981  surfaceScalarField SfUfbyDelta =
2982  fluidsMesh.surfaceInterpolation::deltaCoeffs()*magPhi;
2983 
2984  const scalar deltaT = runTime.deltaT().value();
2985 
2986  CoNum = max(SfUfbyDelta/fluidsMesh.magSf()).value()*deltaT;
2987 
2988  meanCoNum = (sum(SfUfbyDelta)/sum(fluidsMesh.magSf())).value()*deltaT;
2989 
2990  velMag = max(magPhi/fluidsMesh.magSf()).value();
2991  }
2992 
2993  Info<< "FsiFoam:StepFluidItr: Courant Number mean: " << meanCoNum
2994  << " max: " << CoNum
2995  << " velocity magnitude: " << velMag
2996  << endl;
2997 
2998  fvVectorMatrix UEqn
2999  (
3000  fvm::ddt(U)
3001  + fvm::div(phi, U)
3002  - fvm::laplacian(nu, U)
3003  );
3004 
3005  solve(UEqn == -fvc::grad(p));
3006 
3007  // --- PISO loop
3008  volScalarField rUA = 1.0/UEqn.A();
3009  // Masoud
3010  Info << "FsiFoam:StepFluidItr: Performing " << nCorrPISO << " Pressure corrections." << endl;
3011  // Masoud End
3012 
3013  for (int corr=0; corr<nCorrPISO; corr++)
3014  {
3015  U = rUA*UEqn.H();
3016  phi = (fvc::interpolate(U) & fluidsMesh.Sf());
3017 
3018  adjustPhi(phi, U, p);
3019 
3020  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
3021  {
3022  fvScalarMatrix pEqn
3023  (
3024  fvm::laplacian(rUA, p)
3025  == fvc::div(phi)
3026  );
3027 
3028  pEqn.setReference(pRefCell, pRefValue);
3029  pEqn.solve();
3030 
3031  // Added by Masoud
3032  //Info << "p = " << (this-> p()) << endl;
3033 
3034  if (nonOrth == nNonOrthCorr)
3035  {
3036  phi -= pEqn.flux();
3037  }
3038  }
3039 
3040  {
3041  volScalarField contErr = fvc::div(phi);
3042 
3043  sumLocalContErr = runTime.deltaT().value()*
3044  mag(contErr)().weightedAverage(fluidsMesh.V()).value();
3045 
3046  globalContErr = runTime.deltaT().value()*
3047  contErr.weightedAverage(fluidsMesh.V()).value();
3048 
3049  cumulativeContErr += globalContErr;
3050 
3051  Info<< "FsiFoam:StepFluidItr: time step continuity errors : sum local = " << sumLocalContErr
3052  << ", global = " << globalContErr
3053  << ", cumulative = " << cumulativeContErr
3054  << endl;
3055  }
3056 
3057  U -= rUA*fvc::grad(p);
3058  U.correctBoundaryConditions();
3059  }
3060 
3061  {
3062  Info << "FsiFoam:StepFluidItr: Setting tractionis on solid patch" << endl;
3063 
3064  {
3065 
3066  fsiResidualOld = fsiResidual;
3067  this->UpdateFSISurface(solidPatchPointsDispl);
3068  Info << "FsiFoam:StepFluidItr: solidPatchPointsDispl = " << solidPatchPointsDispl << endl;
3069 
3070  fsiResidual = (solidPatchPointsDispl - fluidPatchPointsDispl);
3071  fsiResidualNorm = ::sqrt(gSum(magSqr(fsiResidual)));
3072 
3073  if (outerCorr == 1)
3074  {
3075  initialFsiResidualNorm = fsiResidualNorm;
3076  }
3077 
3078  fsiResidualNorm /= initialFsiResidualNorm + SMALL;
3079 
3080  Info << "FsiFoam:StepFluidItr: Current fsi residual norm: " << fsiResidualNorm << endl;
3081  }
3082  }
3083  }
3084  while
3085  (
3086  (fsiResidualNorm > outerCorrTolerance)
3087  && (outerCorr < nOuterCorr) //Masoud changed from outerCorr < nOuterCorr
3088  );
3089 
3090 
3091  if (runTime.outputTime())
3092  {
3093  runTime.write();
3094  }
3095  }
3096  return(0);
3097 }
3098 //Masoud End StepFluidItr
3099 
3100 int fsifoam_module::Dump(){
3101  Foam::Time &runTime(this->RunTime());
3102  fvMesh &structuresMesh(this->StructuresMesh());
3103  volSymmTensorField &sigma(this->sigma());
3104  volScalarField sigmaEq
3105  (
3106  IOobject
3107  (
3108  "sigmaEq",
3109  runTime.timeName(),
3110  structuresMesh,
3111  IOobject::NO_READ,
3112  IOobject::AUTO_WRITE
3113  ),
3114  sqrt((3.0/2.0)*magSqr(dev(sigma)))
3115  );
3116 
3117  Info<< "FsiFoam:Dump: Max sigmaEq = " << max(sigmaEq).value()
3118  << endl;
3119 
3120  runTime.write();
3121  return 0;
3122 };
3123 
3128 //static void fsifoam_module::Load(const std::string &name){
3129 void fsifoam_module::Load(const std::string &name){
3130  std::cout << "FsiFoam:Load: Loading FsiFoamModule with name " << name
3131  << "." << std::endl;
3132 
3134  fsifoam_module *module_pointer = new fsifoam_module();
3135  COM_new_window(name, MPI_COMM_NULL);
3136  module_pointer->my_window_name = name;
3137  std::string global_name(name+".global");
3138  COM_new_dataitem(global_name.c_str(),'w',COM_VOID,1,"");
3139  COM_set_object(global_name.c_str(),0,module_pointer);
3140 
3141 
3143  std::vector<COM_Type> types(13,COM_INT);
3144 
3145  types[0] = COM_RAWDATA;
3146  types[2] = COM_VOID;
3147  COM_set_member_function( (name+".InitFoam").c_str(),
3148  (Member_func_ptr)(&fsifoam_module::InitFoam),
3149  global_name.c_str(), "biii", &types[0]);
3150 
3151  COM_set_member_function( (name+".RunFoam").c_str(),
3152  (Member_func_ptr)(&fsifoam_module::RunFoam),
3153  global_name.c_str(), "b", &types[0]);
3154 
3155 
3156  COM_set_member_function( (name+".StepFoam").c_str(),
3157  (Member_func_ptr)(&fsifoam_module::StepFoam),
3158  global_name.c_str(), "b", &types[0]);
3159 
3160  COM_set_member_function( (name+".StepFluid").c_str(),
3161  (Member_func_ptr)(&fsifoam_module::StepFluid),
3162  global_name.c_str(), "b", &types[0]);
3163 
3164 /*
3165  types[1] = COM_DOUBLE;
3166  COM_set_member_function( (name+".ModifyEndTime").c_str(),
3167  (Member_func_ptr)(&fsifoam_module::ModifyEndTime),
3168  global_name.c_str(), "bi", &types[0]);
3169 */
3170 
3171 
3172  COM_window_init_done(name);
3173 }
3174 
3175 
3176 void fsifoam_module::InitFoam(int *pargc, void **pargv, int *verbIn)
3177 {
3178  int argc = *pargc;
3179  char** argv = (char**)(pargv);
3180  verbosity = *verbIn;
3181 
3182  std::cout << "OF Module verbosity = " << verbosity << std::endl;
3183 
3184 /*
3185  // register function status variables
3186  std::string initStatusName(my_window_name+".initStatus");
3187  COM_new_dataitem(initStatusName.c_str(),'w', COM_INT, 1, "");
3188  COM_set_array(initStatusName.c_str(), 11, &initStatus);
3189 */
3190 
3191  Initialize(argc, argv);
3192 
3194 
3195  // set up solution meta data
3196  Solution().Meta().AddField("time", 's', 1, 8, "s");
3197  Solution().Meta().AddField("endTime", 's', 1, 8, "s");
3198  Solution().Meta().AddField("initStatus", 's', 1, 4, "");
3199  Solution().Meta().AddField("runStatus", 's', 1, 4, "");
3200  Solution().Meta().AddField("dt", 's', 1, 8, "s");
3201  Solution().Meta().AddField("solidDisplacement", 'n', 3, 8, "m");
3202  Solution().Meta().AddField("pressure", 'c', 1, 8, "Pa");
3203  Solution().Meta().AddField("traction", 'c', 3, 8, "N");
3204 
3205  Mesh().nc.init(numPointsSurface, &surfaceCoordinates[0]);
3206  // this assumes that our elements are quads...we can generalize this
3207  Mesh().con.AddElements(numElementsSurface, 4, surfaceConnectivity);
3208 
3209  // create buffers for the actual data
3210  CreateSoln();
3211 
3212  //unsigned int nnodes = Mesh().nc.NNodes();
3213  //unsigned int nelem = Mesh().con.Nelem();
3214 
3215  // allocate our own local buffers for the data
3216 
3217  time.resize(1, -1.0);
3218  endTime.resize(1, -1.0);
3219  UpdateTime();
3220 
3221  initStatus.resize(1, -1000);
3222  runStatus.resize(1, -1000);
3223 
3224  surfacePressure.resize(numElementsSurface, -1);
3225  surfaceTraction.resize(3*numElementsSurface, -1);
3226  solidDisplacement.resize(3*numPointsSurface,0.);
3227 
3228  // reset the buffers to be our own local buffers
3229  Solution().SetFieldBuffer("initStatus", initStatus);
3230  Solution().SetFieldBuffer("runStatus", runStatus);
3231  Solution().SetFieldBuffer("time", time);
3232  Solution().SetFieldBuffer("endTime", endTime);
3233  Solution().SetFieldBuffer("pressure", surfacePressure);
3234  Solution().SetFieldBuffer("traction", surfaceTraction);
3235  Solution().SetFieldBuffer("solidDisplacement", solidDisplacement);
3236 
3237  SolverUtils::RegisterSolverInto(my_window_name.c_str(), *this);
3238 
3239  initStatus[0] = 0;
3240  return;
3241 }
3242 
3243 
3248 //static void fsifoam_module::Unload(const std::string &name){
3249 void fsifoam_module::Unload(const std::string &name){
3250  std::cout << "FsiFoam:Unload: Unloading FsiFoamModule with name " << name
3251  << "." << std::endl;
3252  fsifoam_module *module_pointer = NULL;
3253  std::string global_name(name+".global");
3254  COM_get_object(global_name.c_str(),0,&module_pointer);
3255  COM_assertion_msg( module_pointer->validate_object()==0, "Invalid object");
3256  delete module_pointer;
3257  COM_delete_window(std::string(name));
3258 }
3259 
3263 void fsifoam_module::RunFoam(){
3264 
3265 
3266  Foam::Time &runTime(RunTime());
3267 
3268  Info << "\nFsiFoam:RunFoam: Starting time loop\n" << endl;
3269 
3270  while(!runTime.end()){
3271  Info << "FsiFoam:RunFoam: Time = " << runTime.timeName() << nl << endl;
3272  Step();
3273  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
3274  << " ClockTime = " << runTime.elapsedClockTime() << " s"
3275  << endl << endl;
3276  if (runTime.outputTime())
3277  Dump();
3278  }
3279 
3280  Info<< "FsiFoam:RunFoam: End\n" << endl;
3281  runStatus[0] = 0;
3282 
3283  return;
3284 }
3285 
3289 void fsifoam_module::StepFoam(){
3290 
3291 
3292  Foam::Time &runTime(RunTime());
3293 
3294  Info << "\nFsiFoam:StepFoam: Stepping time loop\n" << endl;
3295 
3296  // while(!runTime.end()){
3297  if(!runTime.end()){
3298  Info << "FsiFoam:StepFoam: Time = " << runTime.timeName() << nl << endl;
3299  Step();
3300  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
3301  << " ClockTime = " << runTime.elapsedClockTime() << " s"
3302  << endl << endl;
3303  if (runTime.outputTime())
3304  Dump();
3305 
3306  // update IMPACT time
3307  UpdateTime();
3308  UpdateFSISurfaceData();
3309  //Masoud
3310  UpdateFSISurfaceMesh();
3311  //Masoud : End
3312 
3313  // Masoud : cheking impact data
3314  Info << "FsiFoam:StepFoam: numPointSurface = " << numPointsSurface << endl;
3315  Info << "FsiFoam:StepFoam: numElementSurface = " << numElementsSurface << endl;
3316  Info << "FsiFoam:StepFoam: SurfaceCoordinates = " << endl;
3317  for (int i = 0; i < numPointsSurface; i++)
3318  {
3319  Info << surfaceCoordinates[i*3] << " "
3320  << surfaceCoordinates[i*3+1] << " "
3321  << surfaceCoordinates[i*3+2] << endl;
3322  }
3323 
3324  // Masoud : end
3325  }
3326  if(runTime.end()){
3327  Info<< "End\n" << endl;
3328  }
3329  runStatus[0] = 0;
3330  return;
3331 }
3332 
3336 void fsifoam_module::StepFluid(){
3337 
3338  // Masoud: probe output file and preps
3339  ofstream prb2DmpFile;
3340  dynamicFvMesh &fluidsMesh(this->FluidMesh());
3341  pointField foamCoords = fluidsMesh.allPoints();
3342  // Masoud: end
3343 
3344  Foam::Time &runTime(RunTime());
3345 
3346  Info << "FsiFoam:StepFluid: Stepping fluid solver..." << endl;
3347 
3348  // while(!runTime.end()){
3349  if(!runTime.end()){
3350  Info << "FsiFoam:StepFluid: Time = " << runTime.timeName() << endl;
3351  // Original: This step routine does perform iterations properly
3352  // StepFluidAlone();
3353  // Original: End
3354 
3355  // Masoud: Non-iterative stepping
3356  StepFluidNonItr();
3357  // Masoud: end
3358 
3359  // Masoud: Iterative stepping
3360  //StepFluidItr();
3361  // Masoud: end
3362 
3363  Info<< "FsiFoam:StepFluid: ExecutionTime = " << runTime.elapsedCpuTime() << " s"
3364  << " ClockTime = " << runTime.elapsedClockTime() << " s"
3365  << endl;
3366  if (runTime.outputTime())
3367  Dump();
3368 
3369  // update IMPACT data
3370  UpdateTime();
3371  UpdateFSISurfaceData(); // Masoud: updates tractions and pressures
3372  UpdateFSISurfaceMesh(); // Masoud: updates surface mesh
3373 
3374  // Masoud : cheking impact data
3375  //Info << "numPointSurface = " << numPointsSurface << endl;
3376  //Info << "numElementSurface = " << numElementsSurface << endl;
3377  //Info << "SurfaceCoordinates = " << endl;
3378  //for (int i = 0; i < numPointsSurface; i++)
3379  //{
3380  // Info << surfaceCoordinates[i*3] << " "
3381  // << surfaceCoordinates[i*3+1] << " "
3382  // << surfaceCoordinates[i*3+2] << endl;
3383  //}
3384  // Masoud : end
3385 
3386  // Masoud: output for probe
3387  prb2DmpFile.open("prb2.dat", std::ios::app);
3388  prb2DmpFile << runTime.timeName() << " " << \
3389  foamCoords[64].x() << " " << \
3390  foamCoords[64].y() << " " << \
3391  foamCoords[64].z() << "\n";
3392  prb2DmpFile.close();
3393  // Masoud: end
3394 
3395  }
3396  if(runTime.end()){
3397  Info<< "FsiFoam:StepFluid: End\n" << endl;
3398  }
3399  runStatus[0] = 0;
3400  return;
3401 }
3402 
3403 /*
3404 void fsifoam_module::ModifyEndTime(const double &endTime){
3405  RunTime().setEndTime(endTime);
3406 }
3407 */
3408 
3409 void fsifoam_module::CreateFSISurfaceMesh(){
3410 
3411  dynamicFvMesh &fluidsMesh(this->FluidMesh());
3412  numPointsSurface = fluidsMesh.boundaryMesh()[fluidPatchID].nPoints();
3413  //numElementsSurface = fluidsMesh.boundaryMesh()[fluidPatchID].faceCentres().nPoints();
3414 
3415  // face storage
3416  const faceList& faces = fluidsMesh.faces();
3417  const labelList& owner = fluidsMesh.faceOwner();
3418  const polyBoundaryMesh& boundaryMesh = fluidsMesh.boundaryMesh();
3419  const label start = boundaryMesh[fluidPatchID].start();
3420  const label end = start + boundaryMesh[fluidPatchID].size();
3421  numElementsSurface = end - start;
3422 
3423  // point storage
3424  const pointField points = fluidsMesh.points();
3425 
3426  // node 1 and element 1 live at index zero
3427  int nodeNumber = 1;
3428  for (label faceI = start; faceI < end; ++faceI) {
3429  const labelList& vrtList = faces[faceI];
3430  forAll(vrtList, i) {
3431  const label& vertex = vrtList[i];
3432 
3433  // create forward and reverse mapping for this vertex if not already in map
3434  if (!surfaceNodeMap[vertex]) {
3435  surfaceNodeMap[vertex] = nodeNumber;
3436  interfaceToFoamNodeMap[nodeNumber] = vertex;
3437  nodeNumber++;
3438  // add to node list
3439  surfaceCoordinates.push_back(points[vertex].x());
3440  surfaceCoordinates.push_back(points[vertex].y());
3441  surfaceCoordinates.push_back(points[vertex].z());
3442  }
3443  // build local element connectivity
3444  surfaceConnectivity.push_back(vertex);
3445  }
3446  }
3447 
3448  // build map for OpenFoam global vertices to OpenFoam patch vertices
3449  const labelList& meshPoints =
3450  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
3451  for (int i=0; i<numPointsSurface; ++i) {
3452  foamGlobalToPatchNodeMap[meshPoints[i]] = i;
3453  }
3454 
3455  // translate OpenFOAM global node numbers to local surface node numbers
3456  // node 1 and element 1 live at index zero
3457  for (unsigned int i=0; i<surfaceConnectivity.size(); ++i) {
3458  surfaceConnectivity[i] = surfaceNodeMap[surfaceConnectivity[i]];
3459  }
3460 }
3461 
3462 void fsifoam_module::UpdateFSISurfaceMesh(){
3463 
3464  dynamicFvMesh &fluidsMesh(this->FluidMesh());
3465  pointField foamCoords = fluidsMesh.allPoints();
3466 
3467  for (int i=0; i<numPointsSurface; ++i) {
3468  int vertex = interfaceToFoamNodeMap[i+1];
3469  surfaceCoordinates[3*i] = foamCoords[vertex].x();
3470  surfaceCoordinates[3*i+1] = foamCoords[vertex].y();
3471  surfaceCoordinates[3*i+2] = foamCoords[vertex].z();
3472  }
3473 }
3474 
3475 void fsifoam_module::UpdateTime(){
3476  Foam::Time &runTime(this->RunTime());
3477  time[0] = runTime.value();
3478  endTime[0] = runTime.endTime().value();
3479 }
3480 
3481 int fsifoam_module::UpdateFSISurface(Foam::vectorField &solidDispl){
3482  //Masoud : Testing displacement
3483  std::cout <<
3484  "FsiFoam:UpdateFSISurface: Reporting solid displacements to the fluid solver" << std::endl;
3485  //std::cout<< "Displacements passed to me : " << std::endl;
3486  // Masoud : End
3488  // final displacement of the fluid surface
3490  for(int i=0; i<numPointsSurface; ++i) {
3491  solidDispl[foamGlobalToPatchNodeMap[interfaceToFoamNodeMap[i+1]]].x() =
3492  solidDisplacement[3*i];
3493  solidDispl[foamGlobalToPatchNodeMap[interfaceToFoamNodeMap[i+1]]].y() =
3494  solidDisplacement[3*i+1];
3495  solidDispl[foamGlobalToPatchNodeMap[interfaceToFoamNodeMap[i+1]]].z() =
3496  solidDisplacement[3*i+2];
3497  //Masoud : Testing displacement
3498  //std::cout<< solidDisplacement[3*i] << " "
3499  // << solidDisplacement[3*i+1] << " "
3500  // << solidDisplacement[3*i+2] << std::endl;
3501  //Masoud: End
3502  }
3503  return(0);
3504 }
3505 
3506 void fsifoam_module::UpdateFSISurfaceData(){
3507 
3508  //Masoud
3509  double currTime((this->RunTime()).value());
3510  Info << "FsiFoam:UpdateFSISurface: Reporting surface loads to the solid solver " << endl;
3511 
3513  // pressure on the surface
3515  dimensionedScalar &rhoFluid(this->rhoFluid());
3516  volScalarField &p(this->p());
3517 
3518  scalarField foamSurfacePressure = rhoFluid.value()*p.boundaryField()[fluidPatchID];
3519 
3520  for(int i=0; i<numElementsSurface; ++i) {
3521  surfacePressure[i] = foamSurfacePressure[i];
3522  }
3523 
3525  // traction on the surface
3527  dimensionedScalar &nu(this->nu());
3528  volVectorField &U(this->U());
3529  // Original
3530  //vectorField foamSurfaceTraction = -rhoFluid.value()*nu.value()*
3531  // U.boundaryField()[fluidPatchID].snGrad();
3532  // Original end
3533 
3534  // Masoud change: based on following taken from StepFoam()
3535  // vectorField fluidPatchTraction =
3536  // -rhoFluid.value()*nu.value()
3537  // *U.boundaryField()[fluidPatchID].snGrad()
3538  // + rhoFluid.value()*p.boundaryField()[fluidPatchID]
3539  // *fluidsMesh.boundary()[fluidPatchID].nf();
3540 
3541  dynamicFvMesh &fluidsMesh(this->FluidMesh());
3542  vectorField foamSurfaceTraction = -rhoFluid.value()*nu.value()
3543  *U.boundaryField()[fluidPatchID].snGrad()
3544  + rhoFluid.value()*p.boundaryField()[fluidPatchID]
3545  *fluidsMesh.boundary()[fluidPatchID].nf();
3546  //Masoud: end
3547 
3548  //Masoud: sending tractions after some initial time
3549  if (currTime > -0.1) {
3550  for(int i=0; i<numElementsSurface; ++i) {
3551  surfaceTraction[3*i] = foamSurfaceTraction[i].x();
3552  surfaceTraction[3*i+1] = foamSurfaceTraction[i].y();
3553  surfaceTraction[3*i+2] = foamSurfaceTraction[i].z();
3554  }
3555  } else {
3556  for(int i=0; i<numElementsSurface; ++i) {
3557  surfaceTraction[3*i] = 0.0;
3558  surfaceTraction[3*i+1] = 0.0;
3559  surfaceTraction[3*i+2] = 0.0;
3560  }
3561  }
3562  //Masoud: end
3563 }
3564 
3568 extern "C" void OpenFoamFSI_load_module( const char *name) {
3569  fsifoam_module::Load(name);
3570 }
3574 extern "C" void OpenFoamFSI_unload_module( const char *name){
3575  fsifoam_module::Unload(name);
3576 }
3577 
fsifoam_module::UpdateFSISurfaceData
void UpdateFSISurfaceData()
Update the data registered on the FSI surface.
Definition: fsiFOAM.C:3514
fsifoam_module::StepFoam
void StepFoam()
the OpenFOAM stepping
Definition: fsiFOAM.C:3297
fsifoam_module::StepFluidItr
int StepFluidItr()
Definition: fsiFOAM.C:2674
fsifoam_module::InitFoam
void InitFoam(int *pargc, void **pargv, int *verbIn)
function to register through IMPACT
Definition: fsiFOAM.C:3184
fsifoam_module
Definition: fsiFOAM.H:38
fsifoam_module::InitFluidMesh
int InitFluidMesh()
Definition: fsiFOAM.C:52
fsifoam_module::UpdateFSISurface
int UpdateFSISurface(Foam::vectorField &solidPatchPointsDispl)
Update the nodal coordinates of the IMAPCT and OpenFoam FSI surfaces from the IMPACT displacement dat...
Definition: fsiFOAM.C:3489
fsifoam_module::UpdateTime
void UpdateTime()
Update the time control data registered with OpenFOAM i.e.
Definition: fsiFOAM.C:3483
fsifoam_module::Load
static void Load(const std::string &name)
&quot;Loads&quot; IcoFoamModule
Definition: fsiFOAM.C:3137
fsifoam_module::Unload
static void Unload(const std::string &name)
Unloads the IcoFoamModule.
Definition: fsiFOAM.C:3257
fsifoam_module::ReadPISOControls
int ReadPISOControls()
Definition: fsiFOAM.C:629
fsifoam_module::Dump
int Dump()
Definition: fsiFOAM.C:3108
fsifoam_module::RunFoam
void RunFoam()
the OpenFOAM main
Definition: fsiFOAM.C:3271
OpenFoamFSI_load_module
void OpenFoamFSI_load_module(const char *name)
C/C++ bindings to load IcoFoamModule.
Definition: fsiFOAM.C:3576
fsifoam_module::Initialize
int Initialize(int argc, char *argv[])
Definition: fsiFOAM.C:5
fsifoam_module::my_window_name
std::string my_window_name
Definition: fsiFOAM.H:41
fsifoam_module::ReadFSIControls
int ReadFSIControls()
Definition: fsiFOAM.C:642
fsifoam_module::StepFluid
void StepFluid()
the OpenFOAM stepping fluid alone
Definition: fsiFOAM.C:3344
fsifoam_module::UpdateFSISurfaceMesh
void UpdateFSISurfaceMesh()
Update of the surface mesh data coordinates.
Definition: fsiFOAM.C:3470
fsifoam_module::CreateInterZoneInterpolators
int CreateInterZoneInterpolators()
Definition: fsiFOAM.C:350
OpenFoamFSI_unload_module
void OpenFoamFSI_unload_module(const char *name)
C/C++ bindings to unload IcoFoamModule.
Definition: fsiFOAM.C:3582
fsifoam_module::StepFluidNonItr
int StepFluidNonItr()
Definition: fsiFOAM.C:2310
fsifoam_module::CreateFSISurfaceMesh
void CreateFSISurfaceMesh()
Initialization of the surface mesh data structures from the OpenFOAM mesh data structures.
Definition: fsiFOAM.C:3417
fsifoam_module::InitTransportProperties
int InitTransportProperties()
Definition: fsiFOAM.C:68
fsifoam_module::ReadCouplingProperties
int ReadCouplingProperties()
Definition: fsiFOAM.C:204
fsifoam_module::Step
int Step()
Definition: fsiFOAM.C:703
fsifoam_module::InitStructuresMesh
int InitStructuresMesh()
Definition: fsiFOAM.C:127
fsifoam_module::CreateFluidFields
int CreateFluidFields()
Definition: fsiFOAM.C:80
fsifoam_module::CreateStructuresFields
int CreateStructuresFields()
Definition: fsiFOAM.C:133
fsifoam_module::FindGlobalFaceZones
int FindGlobalFaceZones()
Definition: fsiFOAM.C:514
fsifoam_module::StepFluidAlone
int StepFluidAlone()
Definition: fsiFOAM.C:1854