| Line |
Branch |
Exec |
Source |
| 1 |
|
|
// ______ ______ _ _ _____ ______ |
| 2 |
|
|
// | ____| ____| | (_)/ ____| | ____| |
| 3 |
|
|
// | |__ | |__ | | _| (___ ___| |__ |
| 4 |
|
|
// | __| | __| | | | |\___ \ / __| __| |
| 5 |
|
|
// | | | |____| |____| |____) | (__| |____ |
| 6 |
|
|
// |_| |______|______|_|_____/ \___|______| |
| 7 |
|
|
// Finite Elements for Life Sciences and Engineering |
| 8 |
|
|
// |
| 9 |
|
|
// License: LGL2.1 License |
| 10 |
|
|
// FELiScE default license: LICENSE in root folder |
| 11 |
|
|
// |
| 12 |
|
|
// Main authors: M. A. Fernandez & F.M. Gerosa |
| 13 |
|
|
// |
| 14 |
|
|
|
| 15 |
|
|
/*! |
| 16 |
|
|
\file linearProblemNitscheXFEM.cpp |
| 17 |
|
|
\authors M. Fernandez & F.M. Gerosa |
| 18 |
|
|
\date 04/2018 |
| 19 |
|
|
\brief Model class for Navier-Stokes with Nitsche XFEM formulation |
| 20 |
|
|
*/ |
| 21 |
|
|
|
| 22 |
|
|
// System includes |
| 23 |
|
|
|
| 24 |
|
|
// External includes |
| 25 |
|
|
|
| 26 |
|
|
// Project includes |
| 27 |
|
|
#include "FiniteElement/elementVector.hpp" |
| 28 |
|
|
#include "Model/NitscheXFEMModel.hpp" |
| 29 |
|
|
|
| 30 |
|
|
namespace felisce |
| 31 |
|
|
{ |
| 32 |
|
|
|
| 33 |
|
✗ |
NitscheXFEMModel::NitscheXFEMModel():Model() |
| 34 |
|
|
{ |
| 35 |
|
✗ |
m_name = "Nitsche XFEM Navier-Stokes Model"; |
| 36 |
|
✗ |
m_isParFluidVelAllocated = false; |
| 37 |
|
✗ |
m_pLinFluid = nullptr; |
| 38 |
|
✗ |
m_feStruc = nullptr; |
| 39 |
|
✗ |
m_indexTime = 0; |
| 40 |
|
✗ |
m_countToInitOldObjects = 0; |
| 41 |
|
✗ |
m_reinitializeFluidBdf = false; |
| 42 |
|
✗ |
m_intfDisp = nullptr; |
| 43 |
|
|
} |
| 44 |
|
|
|
| 45 |
|
|
/***********************************************************************************/ |
| 46 |
|
|
/***********************************************************************************/ |
| 47 |
|
|
|
| 48 |
|
✗ |
NitscheXFEMModel::~NitscheXFEMModel() |
| 49 |
|
|
{ |
| 50 |
|
|
|
| 51 |
|
✗ |
if( m_isParFluidVelAllocated ) { |
| 52 |
|
✗ |
for(std::size_t i = 0; i < m_parFluidVel.size(); ++i) |
| 53 |
|
✗ |
m_parFluidVel[i].destroy(); |
| 54 |
|
✗ |
m_parFluidVel.clear(); |
| 55 |
|
|
} |
| 56 |
|
|
|
| 57 |
|
✗ |
if( m_feStruc != NULL ) |
| 58 |
|
✗ |
delete m_feStruc; |
| 59 |
|
|
} |
| 60 |
|
|
|
| 61 |
|
|
/***********************************************************************************/ |
| 62 |
|
|
/***********************************************************************************/ |
| 63 |
|
|
|
| 64 |
|
✗ |
void NitscheXFEMModel::initializeDerivedLinearProblem() |
| 65 |
|
|
{ |
| 66 |
|
|
// type conversion of the linear problems |
| 67 |
|
✗ |
m_pLinFluid = dynamic_cast<LinearProblemNitscheXFEM*>(m_linearProblem[0]); |
| 68 |
|
|
|
| 69 |
|
✗ |
std::cout << " Reading interface mesh... " ; |
| 70 |
|
✗ |
std::string inputDirectory1 = "./"; |
| 71 |
|
✗ |
std::string inputFile1 = ""; |
| 72 |
|
✗ |
std::string inputMesh1 = FelisceParam::instance().meshFileImmersedStruct; |
| 73 |
|
✗ |
std::string outputMesh1 = "toto.geo"; |
| 74 |
|
✗ |
std::string meshDir1 = FelisceParam::instance().meshDirImmersedStruct; |
| 75 |
|
✗ |
std::string resultDir1 = ""; |
| 76 |
|
✗ |
std::string prefixName1 = ""; |
| 77 |
|
|
|
| 78 |
|
✗ |
IO io(inputDirectory1, inputFile1, inputMesh1, outputMesh1, meshDir1, resultDir1, prefixName1); |
| 79 |
|
|
|
| 80 |
|
✗ |
io.readMesh(m_interfMesh, 1.0); |
| 81 |
|
✗ |
std::cout << " done. " << std::endl; // TODO D.C. the model can already have more than one mesh... |
| 82 |
|
|
|
| 83 |
|
✗ |
m_pLinFluid->setInterfaceMesh(&m_interfMesh); |
| 84 |
|
|
|
| 85 |
|
✗ |
if( !FelisceParam::instance().duplicateSupportDof ) |
| 86 |
|
✗ |
FEL_ERROR("NXFEM model requires duplicateSupportDof = true"); |
| 87 |
|
|
|
| 88 |
|
|
// build the edges or faces of the fluid mesh |
| 89 |
|
✗ |
if( m_mesh[0]->domainDim() == 2 ) |
| 90 |
|
✗ |
m_mesh[0]->buildEdges(); |
| 91 |
|
|
else |
| 92 |
|
✗ |
m_mesh[0]->buildFaces(); |
| 93 |
|
|
|
| 94 |
|
|
// We intersect the two meshes here (because this function is called before computeDof) |
| 95 |
|
|
// The meshes can be found in the linear problem |
| 96 |
|
✗ |
m_intersectMeshes.setMeshesName("fluidMesh", "strucMesh"); |
| 97 |
|
✗ |
m_intersectMeshes.setVerbose(FelisceParam::verbose()); |
| 98 |
|
✗ |
m_intersectMeshes.initAndIntersectMeshes(*m_mesh[0], m_interfMesh); |
| 99 |
|
|
|
| 100 |
|
|
// give access to the intersection to each linear problems |
| 101 |
|
✗ |
m_pLinFluid->setDuplicateSupportObject(&m_intersectMeshes); |
| 102 |
|
|
|
| 103 |
|
|
// compute normals |
| 104 |
|
✗ |
m_interfMesh.computeElementNormal(m_interfMesh.listElementNormals()); |
| 105 |
|
|
} |
| 106 |
|
|
|
| 107 |
|
|
/***********************************************************************************/ |
| 108 |
|
|
/***********************************************************************************/ |
| 109 |
|
|
|
| 110 |
|
✗ |
void NitscheXFEMModel::initializeDerivedModel() |
| 111 |
|
|
{ |
| 112 |
|
|
|
| 113 |
|
✗ |
if (FelisceParam::instance().orderBdfNS > 2) |
| 114 |
|
✗ |
FEL_ERROR("BDF not yet implemented for order greater than 2 with Nitsche-XFEM."); |
| 115 |
|
|
|
| 116 |
|
|
// Definition and initialization of bdf schemes |
| 117 |
|
✗ |
m_bdfFluid.defineOrder(FelisceParam::instance().orderBdfNS); |
| 118 |
|
✗ |
m_pLinFluid->initializeTimeScheme(&m_bdfFluid); |
| 119 |
|
|
|
| 120 |
|
|
// Finite element for the linear problems |
| 121 |
|
|
// We assume that there is only one type of volume element for the structure |
| 122 |
|
✗ |
const std::vector<GeometricMeshRegion::ElementType>& bagElementTypeDomainStruc = m_interfMesh.bagElementTypeDomain(); |
| 123 |
|
✗ |
FEL_ASSERT(bagElementTypeDomainStruc.size() == 1); |
| 124 |
|
|
|
| 125 |
|
✗ |
const int typeOfFiniteElement = 0; |
| 126 |
|
✗ |
const GeoElement *geoEleStruc = GeometricMeshRegion::eltEnumToFelNameGeoEle[bagElementTypeDomainStruc[0]].second; |
| 127 |
|
✗ |
const RefElement *refEleStruc = geoEleStruc->defineFiniteEle(bagElementTypeDomainStruc[0], typeOfFiniteElement, m_interfMesh); |
| 128 |
|
✗ |
m_feStruc = new CurvilinearFiniteElement(*refEleStruc, *geoEleStruc, DegreeOfExactness_2); // WARNING: Provide the value with the fluid or by data file |
| 129 |
|
|
|
| 130 |
|
✗ |
FEL_ASSERT(m_feStruc); |
| 131 |
|
✗ |
m_pLinFluid->setStrucFiniteElement(m_feStruc); // TODO D.C. remove this part once merged multimesh in master |
| 132 |
|
|
|
| 133 |
|
|
// We assume that there is only one type of volume element for the fluid |
| 134 |
|
✗ |
FEL_ASSERT(m_pLinFluid->meshLocal()->bagElementTypeDomain().size() == 1); |
| 135 |
|
|
} |
| 136 |
|
|
|
| 137 |
|
|
/***********************************************************************************/ |
| 138 |
|
|
/***********************************************************************************/ |
| 139 |
|
|
|
| 140 |
|
✗ |
void NitscheXFEMModel::preAssemblingMatrixRHS(std::size_t /*iProblem*/) |
| 141 |
|
|
{ |
| 142 |
|
|
|
| 143 |
|
|
// Allocate the extrapolation of the velocity if necessary |
| 144 |
|
✗ |
if( !m_isParFluidVelAllocated ) { |
| 145 |
|
|
// There are orderBdfNS term in the velocity extrapolation |
| 146 |
|
✗ |
m_parFluidVel.resize(FelisceParam::instance().orderBdfNS); |
| 147 |
|
|
|
| 148 |
|
|
// If we use a projection of the old solution, all the terms have the same size |
| 149 |
|
✗ |
if( FelisceParam::instance().useProjectionForNitscheXFEM || m_fstransient->iteration == 0 ) { |
| 150 |
|
✗ |
for(std::size_t i = 0; i < m_parFluidVel.size(); ++i) |
| 151 |
|
✗ |
m_parFluidVel[i].duplicateFrom(m_pLinFluid->vector()); |
| 152 |
|
|
} else { |
| 153 |
|
✗ |
for(std::size_t i = 0; i < m_parFluidVel.size(); ++i) |
| 154 |
|
✗ |
m_parFluidVel[i].duplicateFrom(m_pLinFluid->solutionOld(i)); |
| 155 |
|
|
} |
| 156 |
|
|
|
| 157 |
|
|
// set the extrapolation term to zero |
| 158 |
|
✗ |
for(std::size_t i = 0; i < m_parFluidVel.size(); ++i) |
| 159 |
|
✗ |
m_parFluidVel[i].set( 0.); |
| 160 |
|
|
|
| 161 |
|
✗ |
m_isParFluidVelAllocated = true; |
| 162 |
|
|
} |
| 163 |
|
|
} |
| 164 |
|
|
|
| 165 |
|
|
/***********************************************************************************/ |
| 166 |
|
|
/***********************************************************************************/ |
| 167 |
|
|
|
| 168 |
|
✗ |
void NitscheXFEMModel::forward() |
| 169 |
|
|
{ |
| 170 |
|
✗ |
prepareForward(); |
| 171 |
|
✗ |
solve(); |
| 172 |
|
|
} |
| 173 |
|
|
|
| 174 |
|
|
/***********************************************************************************/ |
| 175 |
|
|
/***********************************************************************************/ |
| 176 |
|
|
|
| 177 |
|
✗ |
void NitscheXFEMModel::prepareForward() |
| 178 |
|
|
{ |
| 179 |
|
|
|
| 180 |
|
|
// Write solution for postprocessing (if required) |
| 181 |
|
✗ |
if ( m_fstransient->iteration > 0 ){ |
| 182 |
|
✗ |
writeSolutionAndMeshes(); |
| 183 |
|
✗ |
if( FelisceParam::verbose() > 4 ) |
| 184 |
|
✗ |
m_pLinFluid->printStructStress(m_fstransient->iteration, m_fstransient->timeStep); |
| 185 |
|
|
} |
| 186 |
|
|
|
| 187 |
|
|
// update the position of the interface and recompute the intersection |
| 188 |
|
✗ |
if ( ( FelisceParam::instance().useDynamicStructure ) and (m_fstransient->iteration > 0) ) { |
| 189 |
|
✗ |
FEL_CHECK(m_intfDisp != 0, "Error: the pointer of the interface force has not been initialized yet."); |
| 190 |
|
|
|
| 191 |
|
✗ |
m_intersectMeshes.updateInterfacePosition( *m_intfDisp ); |
| 192 |
|
|
|
| 193 |
|
✗ |
m_interfMesh.moveMesh( *m_intfDisp, 1.0); |
| 194 |
|
✗ |
m_interfMesh.computeElementNormal(m_interfMesh.listElementNormals()); |
| 195 |
|
|
|
| 196 |
|
✗ |
m_intersectMeshes.setIteration(m_fstransient->iteration); |
| 197 |
|
✗ |
m_intersectMeshes.intersectMeshes(*m_mesh[0], m_interfMesh); |
| 198 |
|
|
} |
| 199 |
|
|
|
| 200 |
|
|
// save the solution with the original support Dof and delete the matrices, vectors and mapping of the fluid |
| 201 |
|
✗ |
if ( m_fstransient->iteration > 0 ) { |
| 202 |
|
|
|
| 203 |
|
✗ |
if( m_intersectMeshes.hasIntersectionChange() ) { |
| 204 |
|
|
|
| 205 |
|
|
// delete dynamic parts |
| 206 |
|
✗ |
m_pLinFluid->updateOldSolution(1); |
| 207 |
|
✗ |
m_pLinFluid->deleteDynamicData(); |
| 208 |
|
|
|
| 209 |
|
✗ |
if(m_isParFluidVelAllocated) { |
| 210 |
|
✗ |
for(std::size_t i=0; i<m_parFluidVel.size(); ++i) |
| 211 |
|
✗ |
m_parFluidVel[i].destroy(); |
| 212 |
|
|
|
| 213 |
|
✗ |
m_isParFluidVelAllocated = false; |
| 214 |
|
|
} |
| 215 |
|
|
|
| 216 |
|
|
// set to true only to remember that the intersection has changed after this time step |
| 217 |
|
|
// If the intersection is not changing at the next time step, we need to reinitialize the old objects |
| 218 |
|
✗ |
m_countToInitOldObjects = FelisceParam::instance().orderBdfNS; |
| 219 |
|
|
} else { |
| 220 |
|
|
|
| 221 |
|
✗ |
m_pLinFluid->updateOldSolution(m_countToInitOldObjects); |
| 222 |
|
|
|
| 223 |
|
✗ |
if(m_countToInitOldObjects > 0) { |
| 224 |
|
|
|
| 225 |
|
✗ |
if(m_isParFluidVelAllocated) { |
| 226 |
|
✗ |
for(std::size_t i=0; i<m_parFluidVel.size(); ++i) |
| 227 |
|
✗ |
m_parFluidVel[i].destroy(); |
| 228 |
|
|
|
| 229 |
|
✗ |
m_isParFluidVelAllocated = false; |
| 230 |
|
|
} |
| 231 |
|
|
|
| 232 |
|
|
// pre assembling to reallocate the extrapolations |
| 233 |
|
✗ |
for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) { |
| 234 |
|
✗ |
preAssemblingMatrixRHS(ipb); |
| 235 |
|
|
} |
| 236 |
|
|
|
| 237 |
|
✗ |
m_countToInitOldObjects -= 1; |
| 238 |
|
✗ |
m_reinitializeFluidBdf = true; |
| 239 |
|
|
} |
| 240 |
|
|
} |
| 241 |
|
|
} |
| 242 |
|
|
|
| 243 |
|
|
// Note that hasIntersectionChange() is true only if the duplicated elements or the signs of the original fluid vertices (for instance, when the structure intersects a new edge or face) are different |
| 244 |
|
✗ |
if( m_intersectMeshes.hasIntersectionChange() ) { |
| 245 |
|
|
// reallocate everything for the fluid |
| 246 |
|
|
// supportDofMeshes, duplication, and dofs |
| 247 |
|
✗ |
m_pLinFluid->computeDof(MpiInfo::numProc(), MpiInfo::rankProc()); |
| 248 |
|
|
|
| 249 |
|
|
// new pattern |
| 250 |
|
✗ |
m_pLinFluid->userChangePattern(MpiInfo::numProc(), MpiInfo::rankProc()); |
| 251 |
|
|
|
| 252 |
|
|
// cut the mesh |
| 253 |
|
✗ |
m_pLinFluid->cutMesh(); |
| 254 |
|
|
|
| 255 |
|
|
// allocate matrices and vectors |
| 256 |
|
✗ |
m_pLinFluid->allocateMatrix(); |
| 257 |
|
|
|
| 258 |
|
|
// determine dof for the boundary conditions |
| 259 |
|
✗ |
m_pLinFluid->determineDofAssociateToLabel(); |
| 260 |
|
|
|
| 261 |
|
✗ |
m_pLinFluid->finalizeEssBCConstantInT(); |
| 262 |
|
|
|
| 263 |
|
|
// pre assembling to reallocate the extrapolations |
| 264 |
|
✗ |
for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) |
| 265 |
|
✗ |
preAssemblingMatrixRHS(ipb); |
| 266 |
|
|
|
| 267 |
|
|
// re build the solver (not possible to change the matrix size in the ksp without destroying it ?) |
| 268 |
|
✗ |
m_pLinFluid->buildSolver(); |
| 269 |
|
|
} |
| 270 |
|
|
|
| 271 |
|
✗ |
if( m_fstransient->iteration == 0 ) { |
| 272 |
|
|
// initialize the old objects |
| 273 |
|
✗ |
m_pLinFluid->initOldObjects(); |
| 274 |
|
|
|
| 275 |
|
✗ |
writeSolutionAndMeshes(); |
| 276 |
|
✗ |
if( FelisceParam::verbose() > 4) |
| 277 |
|
✗ |
m_pLinFluid->printStructStress(m_fstransient->iteration, m_fstransient->timeStep); |
| 278 |
|
|
} |
| 279 |
|
|
|
| 280 |
|
|
// if ( MpiInfo::rankProc() == 0 ){ |
| 281 |
|
|
// // m_ios[ipb]->writeMesh(m_mesh[0]); |
| 282 |
|
|
// m_ios[0]->ensight().writerGeoGold(m_mesh[0], FelisceParam::instance().resultDir+"fluid.geo"); |
| 283 |
|
|
// m_pLinFluid->printMeshPartition(m_fstransient->iteration, m_fstransient->timeStep); |
| 284 |
|
|
// } |
| 285 |
|
|
|
| 286 |
|
|
// Advance time step. |
| 287 |
|
✗ |
updateTime(); |
| 288 |
|
|
|
| 289 |
|
|
// Print time information |
| 290 |
|
✗ |
printNewTimeIterationBanner(); |
| 291 |
|
|
|
| 292 |
|
|
/////////////////////// |
| 293 |
|
|
// Update Bdf scheme // |
| 294 |
|
|
/////////////////////// |
| 295 |
|
✗ |
if( m_intersectMeshes.hasIntersectionChange() || m_reinitializeFluidBdf || m_fstransient->iteration == 1 ) { |
| 296 |
|
|
// initialize the bdf for the fluid |
| 297 |
|
✗ |
PetscVector zero = PetscVector::null(); |
| 298 |
|
✗ |
felInt order = FelisceParam::instance().orderBdfNS; |
| 299 |
|
✗ |
if( order == 1 ) |
| 300 |
|
✗ |
m_bdfFluid.reinitialize(FelisceParam::instance().orderBdfNS, m_pLinFluid->solutionOld(0), zero, zero); |
| 301 |
|
✗ |
else if( order == 2 ) |
| 302 |
|
✗ |
m_bdfFluid.reinitialize(FelisceParam::instance().orderBdfNS, m_pLinFluid->solutionOld(1), m_pLinFluid->solutionOld(0), zero); |
| 303 |
|
|
else |
| 304 |
|
✗ |
FEL_ERROR("This bdf order is not implemented for this problem"); |
| 305 |
|
|
|
| 306 |
|
|
// Compute the extrapolation of the velocity and initialize the one in the solver |
| 307 |
|
✗ |
m_bdfFluid.extrapolateByPart(m_parFluidVel); |
| 308 |
|
✗ |
m_pLinFluid->gatherSeqVelExtrapol(m_parFluidVel); |
| 309 |
|
|
|
| 310 |
|
|
// Compute the rhs coming from the time scheme and initialize the one in the solver |
| 311 |
|
✗ |
m_bdfFluid.computeRHSTimeByPart(m_fstransient->timeStep, m_parFluidVel); |
| 312 |
|
✗ |
m_pLinFluid->gatherSeqBdfRHS(m_parFluidVel); |
| 313 |
|
|
|
| 314 |
|
|
// set flag to false |
| 315 |
|
✗ |
m_reinitializeFluidBdf = false; |
| 316 |
|
✗ |
} else { |
| 317 |
|
✗ |
m_bdfFluid.update(m_pLinFluid->solution()); |
| 318 |
|
|
|
| 319 |
|
✗ |
m_bdfFluid.extrapolateByPart(m_parFluidVel); |
| 320 |
|
✗ |
m_pLinFluid->gatherSeqVelExtrapol(m_parFluidVel); |
| 321 |
|
|
|
| 322 |
|
✗ |
m_bdfFluid.computeRHSTimeByPart(m_fstransient->timeStep, m_parFluidVel); |
| 323 |
|
✗ |
m_pLinFluid->gatherSeqBdfRHS(m_parFluidVel); |
| 324 |
|
|
} |
| 325 |
|
|
|
| 326 |
|
|
// clear the matrices at each sub iteration and computation of the matrix and rhs |
| 327 |
|
✗ |
m_pLinFluid->clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs); |
| 328 |
|
✗ |
m_pLinFluid->vector(1).zeroEntries(); |
| 329 |
|
|
|
| 330 |
|
✗ |
if( FelisceParam::instance().confIntersection ){ |
| 331 |
|
|
|
| 332 |
|
✗ |
m_pLinFluid->assembleMatrixFrontPoints(); |
| 333 |
|
|
} else { |
| 334 |
|
✗ |
m_pLinFluid->assembleMatrixTip(); |
| 335 |
|
✗ |
m_pLinFluid->assembleMatrixTipFaces(); |
| 336 |
|
|
} |
| 337 |
|
|
|
| 338 |
|
✗ |
if(FelisceParam::instance().useGhostPenalty) |
| 339 |
|
✗ |
m_pLinFluid->assembleMatrixGhostPenalty(); |
| 340 |
|
|
|
| 341 |
|
✗ |
if(FelisceParam::instance().NSStabType == 1) |
| 342 |
|
✗ |
m_pLinFluid->assembleMatrixFaceOriented(); |
| 343 |
|
|
|
| 344 |
|
✗ |
if( FelisceParam::instance().contactDarcy > 0 ) |
| 345 |
|
✗ |
m_pLinFluid->assembleMatrixBCDarcy(); |
| 346 |
|
|
|
| 347 |
|
|
// Main assembly loop and assembly of the matrix |
| 348 |
|
|
// The pattern is reduced to where values have been set till here, that is why this loop is called in last |
| 349 |
|
✗ |
m_pLinFluid->assembleMatrixRHS(MpiInfo::rankProc()); |
| 350 |
|
|
|
| 351 |
|
✗ |
m_pLinFluid->addMatrixRHS(); |
| 352 |
|
|
} |
| 353 |
|
|
|
| 354 |
|
|
/***********************************************************************************/ |
| 355 |
|
|
/***********************************************************************************/ |
| 356 |
|
|
|
| 357 |
|
✗ |
void NitscheXFEMModel::solve() |
| 358 |
|
|
{ |
| 359 |
|
|
|
| 360 |
|
|
// Apply boundary conditions. |
| 361 |
|
✗ |
m_pLinFluid->finalizeEssBCTransient(); |
| 362 |
|
✗ |
m_pLinFluid->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc()); |
| 363 |
|
|
|
| 364 |
|
|
// Solve the linear system. |
| 365 |
|
✗ |
PetscPrintf(PETSC_COMM_WORLD, "*** Fluid solver ***\n"); |
| 366 |
|
✗ |
m_pLinFluid->solve(MpiInfo::rankProc(), MpiInfo::numProc()); |
| 367 |
|
|
|
| 368 |
|
|
// Gather the solution to use it in the solid problem |
| 369 |
|
✗ |
m_pLinFluid->gatherSolution(); |
| 370 |
|
|
|
| 371 |
|
|
// evaluation stresses |
| 372 |
|
✗ |
if ( FelisceParam::verbose() ) |
| 373 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Evaluating stress \n"); |
| 374 |
|
|
|
| 375 |
|
✗ |
m_pLinFluid->computeInterfaceStress(); |
| 376 |
|
|
} |
| 377 |
|
|
|
| 378 |
|
|
/***********************************************************************************/ |
| 379 |
|
|
/***********************************************************************************/ |
| 380 |
|
|
|
| 381 |
|
✗ |
void NitscheXFEMModel::solveLinearized(bool linearFlag) |
| 382 |
|
|
{ |
| 383 |
|
|
|
| 384 |
|
✗ |
m_pLinFluid->printSkipVolume(false); |
| 385 |
|
|
|
| 386 |
|
✗ |
m_linearProblem[0]->linearizedFlag() = linearFlag; |
| 387 |
|
|
|
| 388 |
|
✗ |
m_linearProblem[0]->vector().zeroEntries(); |
| 389 |
|
✗ |
m_pLinFluid->vector(1).zeroEntries(); |
| 390 |
|
|
|
| 391 |
|
✗ |
m_linearProblem[0]->assembleMatrixRHS(MpiInfo::rankProc(),FlagMatrixRHS::only_rhs); |
| 392 |
|
|
|
| 393 |
|
✗ |
if ( !linearFlag ) |
| 394 |
|
✗ |
m_linearProblem[0]->addMatrixRHS(); |
| 395 |
|
|
|
| 396 |
|
✗ |
if ( linearFlag && FelisceParam::verbose() ) |
| 397 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"NS linearized solver \n"); |
| 398 |
|
✗ |
else if ( FelisceParam::verbose() ) |
| 399 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"NS solver \n"); |
| 400 |
|
|
|
| 401 |
|
✗ |
if ( FelisceParam::verbose() ) |
| 402 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Applying BC \n"); |
| 403 |
|
|
|
| 404 |
|
|
//Apply essential transient boundary conditions. |
| 405 |
|
✗ |
m_linearProblem[0]->finalizeEssBCTransient(); |
| 406 |
|
|
|
| 407 |
|
✗ |
if ( linearFlag ) // only Dirichlet BC is applied (RHS) |
| 408 |
|
✗ |
m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, |
| 409 |
|
✗ |
MpiInfo::rankProc(), FlagMatrixRHS::only_rhs, FlagMatrixRHS::only_rhs, |
| 410 |
|
|
0, true, ApplyNaturalBoundaryConditions::no); |
| 411 |
|
|
else // all BC applied (RHS): |
| 412 |
|
✗ |
m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, |
| 413 |
|
✗ |
MpiInfo::rankProc(), FlagMatrixRHS::only_rhs, FlagMatrixRHS::only_rhs, |
| 414 |
|
|
0, true, ApplyNaturalBoundaryConditions::yes); |
| 415 |
|
|
|
| 416 |
|
✗ |
if ( linearFlag && FelisceParam::verbose() ) |
| 417 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Solving linearized NS system (same matrix) \n"); |
| 418 |
|
✗ |
else if ( FelisceParam::verbose() ) |
| 419 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Solving NS system (same matrix) \n"); |
| 420 |
|
|
|
| 421 |
|
✗ |
m_linearProblem[0]->solveWithSameMatrix(); |
| 422 |
|
|
|
| 423 |
|
✗ |
m_linearProblem[0]->gatherSolution(); |
| 424 |
|
|
|
| 425 |
|
✗ |
if ( FelisceParam::verbose() ) |
| 426 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Evaluating stress \n"); |
| 427 |
|
|
|
| 428 |
|
✗ |
m_pLinFluid->computeInterfaceStress(); |
| 429 |
|
|
|
| 430 |
|
✗ |
m_pLinFluid->printSkipVolume(true); |
| 431 |
|
|
} |
| 432 |
|
|
|
| 433 |
|
|
/***********************************************************************************/ |
| 434 |
|
|
/***********************************************************************************/ |
| 435 |
|
|
|
| 436 |
|
✗ |
void NitscheXFEMModel::writeSolutionAndMeshes() |
| 437 |
|
|
{ |
| 438 |
|
|
|
| 439 |
|
✗ |
if( (m_fstransient->iteration % FelisceParam::instance().frequencyWriteSolution == 0) ) { |
| 440 |
|
|
|
| 441 |
|
|
// write the mesh |
| 442 |
|
✗ |
if(FelisceParam::verbose()>1) |
| 443 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Write meshes\n"); |
| 444 |
|
|
|
| 445 |
|
✗ |
std::map<felInt, std::vector<felInt> > refToListOfIds; |
| 446 |
|
✗ |
std::string fileName; |
| 447 |
|
✗ |
std::string indexFile; |
| 448 |
|
|
|
| 449 |
|
✗ |
std::size_t ipb=0; |
| 450 |
|
✗ |
for(std::size_t iUnknown=0; iUnknown<m_linearProblem[ipb]->listUnknown().size(); ++iUnknown) { |
| 451 |
|
✗ |
std::stringstream oss; |
| 452 |
|
✗ |
oss << m_indexTime; |
| 453 |
|
✗ |
if(m_indexTime < 10) |
| 454 |
|
✗ |
indexFile = "0000" + oss.str(); |
| 455 |
|
✗ |
else if(m_indexTime < 100) |
| 456 |
|
✗ |
indexFile = "000" + oss.str(); |
| 457 |
|
✗ |
else if(m_indexTime < 1000) |
| 458 |
|
✗ |
indexFile = "00" + oss.str(); |
| 459 |
|
✗ |
else if(m_indexTime < 10000) |
| 460 |
|
✗ |
indexFile = "0" + oss.str(); |
| 461 |
|
|
else |
| 462 |
|
✗ |
indexFile = oss.str(); |
| 463 |
|
|
|
| 464 |
|
✗ |
fileName = FelisceParam::instance().resultDir+FelisceParam::instance().outputMesh[0]+"."+indexFile+".geo"; |
| 465 |
|
|
|
| 466 |
|
|
// fluid, write mesh from the support dof to see the duplication |
| 467 |
|
✗ |
if(iUnknown == 0) { |
| 468 |
|
|
// write the mesh only once (with the velocity, the pressure mesh is the same) |
| 469 |
|
✗ |
if(FelisceParam::instance().outputFileFormat == 0){ |
| 470 |
|
✗ |
m_linearProblem[ipb]->writeGeoForUnknown(iUnknown, fileName); |
| 471 |
|
|
} |
| 472 |
|
✗ |
else if(FelisceParam::instance().outputFileFormat == 1) { |
| 473 |
|
✗ |
m_linearProblem[ipb]->writeGeoForUnknownEnsightGold(iUnknown, refToListOfIds, fileName); |
| 474 |
|
✗ |
m_ios[0]->ensight().setRefToListOfIds(refToListOfIds); |
| 475 |
|
|
} |
| 476 |
|
|
else |
| 477 |
|
✗ |
FEL_ERROR("Unknown output file format"); |
| 478 |
|
|
} |
| 479 |
|
|
} |
| 480 |
|
✗ |
m_meshIsWritten = true; |
| 481 |
|
✗ |
++m_indexTime; |
| 482 |
|
|
|
| 483 |
|
|
// write the solutions |
| 484 |
|
✗ |
if(FelisceParam::verbose() > 1) |
| 485 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Write solutions\n"); |
| 486 |
|
|
|
| 487 |
|
✗ |
m_linearProblem[ipb]->writeSolution(MpiInfo::rankProc(), m_ios, m_fstransient->time, m_fstransient->iteration); |
| 488 |
|
|
|
| 489 |
|
✗ |
if (MpiInfo::rankProc() == 0) |
| 490 |
|
✗ |
m_pLinFluid->writeDuplication(MpiInfo::rankProc(), *m_ios[0], m_fstransient->time, m_fstransient->iteration); |
| 491 |
|
|
|
| 492 |
|
✗ |
if (MpiInfo::rankProc() == 0) { |
| 493 |
|
✗ |
for(std::size_t iio=0; iio<m_ios.size(); ++iio) { |
| 494 |
|
✗ |
if(m_ios[iio]->typeOutput == 1 ) // 1 : ensight |
| 495 |
|
✗ |
m_ios[iio]->postProcess(m_fstransient->time, 2); |
| 496 |
|
|
} |
| 497 |
|
|
} |
| 498 |
|
|
} |
| 499 |
|
|
|
| 500 |
|
✗ |
if(FelisceParam::verbose()>1) |
| 501 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"End Write meshes\n"); |
| 502 |
|
|
} |
| 503 |
|
|
|
| 504 |
|
|
} |
| 505 |
|
|
|