| Directory: | ./ |
|---|---|
| File: | Solver/linearProblemHeat.cpp |
| Date: | 2024-04-14 07:32:34 |
| Exec | Total | Coverage | |
|---|---|---|---|
| Lines: | 46 | 126 | 36.5% |
| Branches: | 35 | 207 | 16.9% |
| 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: J. Foulon & J-F. Gerbeau & V. Martin | ||
| 13 | // | ||
| 14 | |||
| 15 | /*! | ||
| 16 | \file linearProblemHeat.cpp | ||
| 17 | \authors J. Foulon & J-F. Gerbeau & V. Martin | ||
| 18 | \date 05/01/2011 | ||
| 19 | \brief File where is implemented laplacian class. | ||
| 20 | */ | ||
| 21 | |||
| 22 | // System includes | ||
| 23 | |||
| 24 | // External includes | ||
| 25 | |||
| 26 | // Project includes | ||
| 27 | #include "Solver/linearProblemHeat.hpp" | ||
| 28 | #include "Core/felisceTransient.hpp" | ||
| 29 | #include "FiniteElement/elementVector.hpp" | ||
| 30 | #include "FiniteElement/elementMatrix.hpp" | ||
| 31 | #include "InputOutput/io.hpp" | ||
| 32 | #include "DegreeOfFreedom/boundaryCondition.hpp" | ||
| 33 | |||
| 34 | |||
| 35 | namespace felisce { | ||
| 36 | 76 | LinearProblemHeat::LinearProblemHeat(): | |
| 37 | LinearProblem("Heat equation",/*numberOfMatrices*/2), | ||
| 38 |
5/10✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 76 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 76 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 76 times.
✗ Branch 14 not taken.
✓ Branch 17 taken 76 times.
✗ Branch 18 not taken.
|
76 | m_sourceTermDynamicValue(nullptr) |
| 39 | { | ||
| 40 |
2/4✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 76 times.
✗ Branch 6 not taken.
|
76 | m_seqVecs.Init("oldTemperature"); |
| 41 | // When using cvgraph, the following vectors are needed. | ||
| 42 | // If we use a D2N or a N2D some of them may be unused | ||
| 43 | #ifdef FELISCE_WITH_CVGRAPH | ||
| 44 |
2/4✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 76 times.
|
76 | if ( FelisceParam::instance().withCVG ) { |
| 45 | ✗ | m_cvgDirichletVariable = "cvgraphTEMPERATURE"; | |
| 46 | ✗ | m_seqVecs.Init("cvgraphFLOWOFHEAT"); | |
| 47 | ✗ | m_seqVecs.Init("cvgraphTEMPERATURE"); | |
| 48 | } | ||
| 49 | #endif | ||
| 50 | 76 | } | |
| 51 | 2896 | void LinearProblemHeat::updateOldQuantities() { | |
| 52 |
3/6✓ Branch 2 taken 2896 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2896 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 2896 times.
✗ Branch 10 not taken.
|
2896 | m_seqVecs.Get("oldTemperature").copyFrom(sequentialSolution()); |
| 53 | 2896 | } | |
| 54 | ✗ | void LinearProblemHeat::readDataDisplacement(std::vector<IO::Pointer>& io, double iteration) { | |
| 55 | // Read displacement file (*.00000.vct) | ||
| 56 | ✗ | const int iMesh = m_listVariable[0].idMesh(); | |
| 57 | ✗ | m_vectorDisp.resize(m_mesh[iMesh]->numPoints()*3,0.); | |
| 58 | ✗ | io[iMesh]->readVariable(0, iteration,&m_vectorDisp[0], m_vectorDisp.size()); | |
| 59 | } | ||
| 60 | |||
| 61 | 76 | void LinearProblemHeat::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) { | |
| 62 |
1/2✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
|
76 | LinearProblem::initialize(mesh,comm, doUseSNES); |
| 63 | 76 | m_fstransient = fstransient; | |
| 64 |
1/2✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
|
76 | std::vector<PhysicalVariable> listVariable(1); |
| 65 |
1/2✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
|
76 | std::vector<std::size_t> listNumComp(1); |
| 66 | 76 | listVariable[0] = temperature; | |
| 67 | 76 | listNumComp[0] = 1; | |
| 68 | //define unknown of the linear system. | ||
| 69 |
1/2✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
|
76 | m_listUnknown.push_back(temperature); |
| 70 |
1/2✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
|
76 | definePhysicalVariable(listVariable,listNumComp); |
| 71 | |||
| 72 |
1/2✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
|
76 | m_iTemperature = m_listVariable.getVariableIdList(temperature); |
| 73 |
2/4✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 76 times.
✗ Branch 5 not taken.
|
76 | m_sourceTermDynamicValue = FelisceParam::instance().elementFieldDynamicValue("heatSourceTerm"); |
| 74 | 76 | } | |
| 75 | |||
| 76 | 2246 | void LinearProblemHeat::initPerElementType(ElementType /*eltType*/, FlagMatrixRHS /*flagMatrixRHS*/) { | |
| 77 | 2246 | m_feTemp = m_listCurrentFiniteElement[m_iTemperature]; | |
| 78 | |||
| 79 | // define elemField to put term u_n_1.phi_j in RHS. | ||
| 80 | // m_elemField contains u_n_1 values. | ||
| 81 | 2246 | m_elemField.initialize(DOF_FIELD,*m_feTemp); | |
| 82 | 2246 | } | |
| 83 | |||
| 84 | 2168 | void LinearProblemHeat::initPerElementTypeBoundaryCondition(ElementType& /*eltType*/, FlagMatrixRHS /*flagMatrixRHS*/) { | |
| 85 | #ifdef FELISCE_WITH_CVGRAPH | ||
| 86 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 2168 times.
|
2168 | if (FelisceParam::instance().withCVG ) { |
| 87 | ✗ | m_curvFeTemp = m_listCurvilinearFiniteElement[m_iTemperature]; | |
| 88 | ✗ | if ( slave()->thereIsAtLeastOneRobinCondition() ){ | |
| 89 | ✗ | m_robinAux.initialize( DOF_FIELD, *m_curvFeTemp, this->dimension() ); | |
| 90 | } | ||
| 91 | } | ||
| 92 | #endif | ||
| 93 | 2168 | } | |
| 94 | 645100 | void LinearProblemHeat::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/, felInt& iel, FlagMatrixRHS flagMatrixRHS) { | |
| 95 | 645100 | const double coef = 1./m_fstransient->timeStep; | |
| 96 | |||
| 97 | 645100 | m_feTemp->updateFirstDeriv(0, elemPoint); | |
| 98 | |||
| 99 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 645100 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
645100 | if (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_matrix) { |
| 100 |
2/2✓ Branch 1 taken 75212 times.
✓ Branch 2 taken 569888 times.
|
645100 | if ( m_fstransient->iteration == 0 ) { |
| 101 | 75212 | m_elementMat[1]->grad_phi_i_grad_phi_j(1.,*m_feTemp,0,0,1); | |
| 102 | 75212 | m_elementMat[1]->phi_i_phi_j(coef,*m_feTemp,0,0,1); | |
| 103 | } | ||
| 104 | } | ||
| 105 | |||
| 106 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 645100 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
645100 | if (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_rhs) { |
| 107 |
2/2✓ Branch 1 taken 569888 times.
✓ Branch 2 taken 75212 times.
|
645100 | if ( m_fstransient->iteration > 0 ) { |
| 108 | // vector operators. | ||
| 109 |
3/6✓ Branch 3 taken 569888 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 569888 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 569888 times.
✗ Branch 10 not taken.
|
569888 | m_elemField.setValue(m_seqVecs.Get("oldTemperature"), *m_feTemp, iel, m_iTemperature, m_ao, dof()); |
| 110 | 569888 | m_elementVector[0]->source(coef,*m_feTemp,m_elemField,0,1); | |
| 111 | |||
| 112 |
2/2✓ Branch 0 taken 8192 times.
✓ Branch 1 taken 561696 times.
|
569888 | if (m_sourceTermDynamicValue != nullptr) { |
| 113 |
1/2✓ Branch 3 taken 8192 times.
✗ Branch 4 not taken.
|
8192 | m_sourceTerm.setValue(*m_feTemp, *m_sourceTermDynamicValue, m_fstransient->time); |
| 114 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 8192 times.
|
8192 | assert(!m_elementVector.empty()); |
| 115 | 8192 | m_elementVector[0]->source(1.,*m_feTemp,m_sourceTerm,0,1); | |
| 116 | } | ||
| 117 | } | ||
| 118 | } | ||
| 119 | 645100 | } | |
| 120 | |||
| 121 | 70592 | void LinearProblemHeat::computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/, felInt& iel, FlagMatrixRHS /*flagMatrixRHS*/) { | |
| 122 | (void) iel; | ||
| 123 | (void) elemPoint; | ||
| 124 | #ifdef FELISCE_WITH_CVGRAPH | ||
| 125 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 70592 times.
|
70592 | if (FelisceParam::instance().withCVG ) { |
| 126 | ✗ | m_curvFeTemp->updateMeasNormal(0, elemPoint); | |
| 127 | ✗ | cvgraphNaturalBC(iel); | |
| 128 | } | ||
| 129 | #endif | ||
| 130 | 70592 | } | |
| 131 | #ifdef FELISCE_WITH_CVGRAPH | ||
| 132 | ✗ | void LinearProblemHeat::sendData() { | |
| 133 | ✗ | gatherSolution(); | |
| 134 | ✗ | for ( std::size_t iConn(0); iConn<slave()->numConnections(); ++iConn ) { | |
| 135 | ✗ | std::vector<PetscVector> vecs; | |
| 136 | ✗ | for (std::size_t cVar(0); cVar<slave()->numVarToSend(iConn); ++cVar) { | |
| 137 | ✗ | std::string varToSend=slave()->sendVariable(iConn,cVar); | |
| 138 | ✗ | if (varToSend=="FLOWOFHEAT") { | |
| 139 | ✗ | prepareResidual(iConn); | |
| 140 | ✗ | vecs.push_back(m_seqStressOut); | |
| 141 | ✗ | } else if(varToSend =="TEMPERATURE") { | |
| 142 | ✗ | vecs.push_back(this->sequentialSolution()); | |
| 143 | } else { | ||
| 144 | ✗ | std::cout<<"variable to send: "<<varToSend<<std::endl; | |
| 145 | ✗ | FEL_ERROR("Error in variable to read"); | |
| 146 | } | ||
| 147 | } | ||
| 148 | ✗ | slave()->sendData(vecs,iConn); | |
| 149 | } | ||
| 150 | } | ||
| 151 | ✗ | void LinearProblemHeat::readData() { | |
| 152 | ✗ | for ( std::size_t iConn(0); iConn<slave()->numConnections(); ++iConn ) { | |
| 153 | ✗ | std::vector<PetscVector> vecs; | |
| 154 | ✗ | for (std::size_t cVar(0); cVar<slave()->numVarToRead(iConn); ++cVar) { | |
| 155 | ✗ | std::string varToRead=slave()->readVariable(iConn,cVar); | |
| 156 | ✗ | if ( varToRead == "FLOWOFHEAT" ) { | |
| 157 | ✗ | vecs.push_back(m_seqVecs.Get("cvgraphFLOWOFHEAT")); | |
| 158 | ✗ | } else if ( varToRead == "TEMPERATURE" ) { | |
| 159 | ✗ | vecs.push_back(m_seqVecs.Get("cvgraphTEMPERATURE")); | |
| 160 | } else { | ||
| 161 | ✗ | std::cout<<"Connection: "<<iConn<<"variable to read: "<<varToRead<<std::endl; | |
| 162 | ✗ | FEL_ERROR("Error in variable to read"); | |
| 163 | } | ||
| 164 | } | ||
| 165 | ✗ | slave()->receiveData(vecs,iConn); | |
| 166 | } | ||
| 167 | } | ||
| 168 | /*! \brief Function to assemble the mass matrix | ||
| 169 | * Function to be called in the assemblyLoopBoundaryGeneral | ||
| 170 | */ | ||
| 171 | void | ||
| 172 | ✗ | LinearProblemHeat::initPerETMass() { | |
| 173 | ✗ | felInt numDofTotal = 0; | |
| 174 | //pay attention this numDofTotal is bigger than expected since it counts for all the VARIABLES | ||
| 175 | ✗ | for (std::size_t iFe = 0; iFe < this->m_listCurvilinearFiniteElement.size(); iFe++) {//this loop is on variables while it should be on unknown | |
| 176 | ✗ | numDofTotal += this->m_listCurvilinearFiniteElement[iFe]->numDof()*this->m_listVariable[iFe].numComponent(); | |
| 177 | } | ||
| 178 | ✗ | m_globPosColumn.resize(numDofTotal); m_globPosRow.resize(numDofTotal); m_matrixValues.resize(numDofTotal*numDofTotal); | |
| 179 | ✗ | m_curvFeTemp = this->m_listCurvilinearFiniteElement[ m_iTemperature ]; | |
| 180 | } | ||
| 181 | |||
| 182 | ✗ | void LinearProblemHeat::massMatrixComputer(felInt ielSupportDof) { | |
| 183 | ✗ | m_elementMatBD[0]->zero(); | |
| 184 | ✗ | m_elementMatBD[0]->phi_i_phi_j(/*coef*/1.,*m_curvFeTemp,/*iblock*/0,/*iblock*/0,1/*dim*/); | |
| 185 | ✗ | setValueMatrixBD(ielSupportDof); | |
| 186 | } | ||
| 187 | |||
| 188 | ✗ | void LinearProblemHeat::updateFE(const std::vector<Point*>& elemPoint, const std::vector<int>&) { | |
| 189 | ✗ | m_curvFeTemp->updateMeasNormal(0, elemPoint); | |
| 190 | } | ||
| 191 | |||
| 192 | ✗ | void LinearProblemHeat::cvgraphNaturalBC(felInt iel) { | |
| 193 | BoundaryCondition* BC; | ||
| 194 | // Loop over the cvgraph connections | ||
| 195 | ✗ | for ( std::size_t iConn(0); iConn<slave()->numConnections(); ++iConn ) { | |
| 196 | // Extract the labels of the current connection | ||
| 197 | ✗ | std::vector<int> labels = slave()->interfaceLabels(iConn); | |
| 198 | // Verify if the current label is part of the interface of this connection | ||
| 199 | ✗ | if ( std::find( labels.begin(), labels.end(), m_currentLabel ) != labels.end() ) { | |
| 200 | ✗ | switch ( slave()->numVarToRead(iConn) ) { | |
| 201 | ✗ | case 1: // only one var to read: neumann interface BC | |
| 202 | // Verify that we have to apply the stress in its strong form. | ||
| 203 | ✗ | if ( slave()->readVariable(iConn,/*cVar*/0) == "FLOWOFHEAT" && slave()->residualInterpolationType(iConn) == 0 ) { | |
| 204 | // Look through the different neumann conditions to find the correct one. | ||
| 205 | ✗ | for (std::size_t iNeumann=0; iNeumann<m_boundaryConditionList.numNeumannBoundaryCondition(); iNeumann++ ) { | |
| 206 | ✗ | BC = m_boundaryConditionList.Neumann(iNeumann);//get the the BC | |
| 207 | ✗ | const auto& list_bc = BC->listLabel(); | |
| 208 | ✗ | if ( list_bc.find(m_currentLabel) != list_bc.end()) { | |
| 209 | //compute the force term. | ||
| 210 | ✗ | m_elemFieldNeumann[iNeumann].setValue( m_seqVecs.Get("cvgraphFLOWOFHEAT"), *m_curvFeTemp , iel, m_iTemperature, m_ao, dof()); | |
| 211 | } | ||
| 212 | } | ||
| 213 | } | ||
| 214 | ✗ | break; | |
| 215 | ✗ | case 2: //two variables to read: robin interface BC | |
| 216 | // Look through the different robin conditions to find the correct one. | ||
| 217 | ✗ | for (std::size_t iRobin=0; iRobin<m_boundaryConditionList.numRobinBoundaryCondition(); iRobin++ ) { | |
| 218 | ✗ | BC = m_boundaryConditionList.Robin(iRobin); | |
| 219 | ✗ | const auto& list_bc = BC->listLabel(); | |
| 220 | ✗ | if ( list_bc.find(m_currentLabel) != list_bc.end()) { | |
| 221 | ✗ | if ( slave()->residualInterpolationType(iConn) == 0 ) { | |
| 222 | ✗ | m_elemFieldRobin[iRobin].setValue( m_seqVecs.Get("cvgraphFLOWOFHEAT"), *m_curvFeTemp , iel, m_iTemperature, m_ao, dof()); | |
| 223 | } else { | ||
| 224 | ✗ | m_elemFieldRobin[iRobin].val *=0 ; | |
| 225 | } | ||
| 226 | ✗ | m_robinAux.setValue(m_seqVecs.Get("cvgraphTEMPERATURE"),*m_curvFeTemp , iel, m_iTemperature, m_ao, dof()); | |
| 227 | ✗ | m_elemFieldRobin[iRobin].val += FelisceParam::instance().alphaRobin[iRobin]*m_robinAux.val; | |
| 228 | } | ||
| 229 | } | ||
| 230 | ✗ | break; | |
| 231 | ✗ | default: | |
| 232 | ✗ | FEL_ERROR("Three variables to read not yet implemented for this linear problem."); | |
| 233 | ✗ | break; | |
| 234 | } | ||
| 235 | } | ||
| 236 | } | ||
| 237 | } | ||
| 238 | ✗ | void LinearProblemHeat::assembleMassBoundaryAndInitKSP( std::size_t iConn ) { | |
| 239 | ✗ | LinearProblem::assembleMassBoundaryAndInitKSP(&LinearProblemHeat::massMatrixComputer, | |
| 240 | ✗ | this->slave()->interfaceLabels(iConn), | |
| 241 | &LinearProblemHeat::initPerETMass, | ||
| 242 | &LinearProblemHeat::updateFE, | ||
| 243 | iConn); | ||
| 244 | } | ||
| 245 | #endif | ||
| 246 | } | ||
| 247 |