// ************************************************************************ // --- Logiciel MEF++ // --- // --- Copyright 1996-2006 GIREF, Université Laval // --- TOUS DROITS RÉSERVÉS // --- ALL RIGHTS RESERVED // --- // --- Ce logiciel est couvert par la Loi sur le droit d'auteur. // --- L'utilisation, ou la modification de ce logiciel sous toutes ses // --- formes, que ce soit code source ou code compilé, à des fins // --- personnelles et non commerciales sont autorisées sans frais pour // --- autant que la présente notice de copyright ainsi que cette // --- permission apparaissent dans toutes les copies ainsi que dans la // --- documentation. // --- Le GIREF et l'Université Laval ne prétendent en aucune façon que // --- ce code convient à un emploi quelconque. Celui-ci est distribué // --- sans aucune garantie implicite ou explicite. // ************************************************************************ // ******************************************************************** // Nom du fichier: GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin.cc // // ******************************************************************** #include "GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin.h" #include "BlocMElem.h" #include "BlocVElem.h" #include "ChampGeometrique.h" #include "ChampTensorielO2Sym.h" #include "ChampTensorielO4Sym.h" #include "ChampTensorielO4SymMin.h" #include "ChampVectoriel3D.h" #include "FonctionsInterpolation.h" #include "fonctionsTenseurs.h" #include "InfoLocalesChamp.h" #include "InfoLocalesInterpolation.h" #include "InfoLocalesTransformation.h" #include "Maillage.h" #include "SchemaIntg.h" #include "TenseurO2Sym.h" #include "TenseurO4Sym.h" #include // ******************************************************************** // Description: Constructeur // // Entrée: Aucune // // Sortie: Aucune // // ******************************************************************** GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin::GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin() :TFGenerique(), aChampU(0), aChampV(0), aChampPseudoSigmaO4Precedent(0), aChampS(0), aChampT(0), aChampI(0), aVectInfoTransformation(0), aVectInfoInterpolationNU(0), aVectInfoInterpolationNV(0), aVectInfoChampU(0), aVectInfoChampUPrecedent(0), aVectInfoChampPseudoSigmaPrecedent(0), aVectInfoChampS(0), aVectInfoChampT(0), aVectInfoChampI(0) { _GIREF_CONSTRUCTEUR_INVARIANTS; _GIREF_TEST_INVARIANTS; } // ******************************************************************** // // Description: Constructeur par copie // // Entrée: pTF le terme de formulation à copier // // Sortie: Aucune // // ******************************************************************** GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin::GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin(const GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin& pTF) : TFGenerique(pTF) { _GIREF_CONSTRUCTEUR_INVARIANTS; _GIREF_APPEL_TEST_INVARIANTS(pTF); _GIREF_ASSERTION(false, "devrait jamais passer ici"); _GIREF_TEST_INVARIANTS; } // ******************************************************************** // // Description: Destructeur // // Entrée: Aucune // // Sortie: Aucune // // ******************************************************************** GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin::~GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin() { _GIREF_TEST_INVARIANTS; _GIREF_DESTRUCTEUR_INVARIANTS; } // ******************************************************************** // // Description: Opérateur d'assignation // // Entrée: pTF le terme de formulation à copier // // Sortie: une référence sur l'objet courant // // ******************************************************************** GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin& GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin::operator= (const GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin& pTF) { _GIREF_TEST_INVARIANTS; _GIREF_APPEL_TEST_INVARIANTS(pTF); _GIREF_ASSERTION(false, "devrait jamais passer ici"); if(this != &pTF) { } _GIREF_TEST_INVARIANTS; return (*this); } // ******************************************************************** // // Description: Assigne les champs pour les calculs dans ce terme de formulation // // Entrée: ChampGeometrique& pChampGeo: Le champ géométrique pour le calcul // ChampVectoriel3D& pChampU: Le champ inconnu dont on veut calculer // ChampVectoriel3D& pChampV: Le champ de fonctions tests V utilisé dans le calcul // ChampTensorielO4SymMin& pChampPseudoSigmaO4Precedent: Le champ contenant le tenseur des pseudo contraintes O4 au pas precedent // ChampTensorielO4Sym& pChampS: Le champ contenant l'integrale du produit du tenseur d'élasticité // ChampTensorielO4Sym& pChampT: Le champ contenant // et du tenseur obtenu par l'exponentielle (termes à termes) du temps de // relaxation (rapport elasticite/viscosite) // ChampTensorielO4Sym& pChampI: Le champ contenant l'exponentiel de l'intégrale du rapport elasticite viscosite // SchemaIntg& pSchemaIntg: Le schéma d'intégration pour ce terme de formulation // // Sortie: Message d'erreur si l'asgnCouplage ne fonctionne pas // // ******************************************************************** ERMsg GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin::asgnChamps(ChampGeometrique& pChampGeo, ChampVectoriel3D& pChampU, ChampVectoriel3D& pChampV, ChampVectoriel3D& pChampUPrecedent, ChampTensorielO4SymMin& pChampPseudoSigmaO4Precedent, ChampTensorielO4Sym& pChampS, ChampTensorielO4Sym& pChampT, ChampTensorielO4Sym& pChampI, const SchemaIntg& pSchemaIntg) { _GIREF_TEST_INVARIANTS; ERMsg lMsg = ERMsg(ERMsg::OK); asgnChampGeometrique(pChampGeo); aChampU = &pChampU; aChampV = &pChampV; aChampUPrecedent = &pChampUPrecedent; aChampPseudoSigmaO4Precedent = &pChampPseudoSigmaO4Precedent; aChampI = &pChampI; aChampS = &pChampS; aChampT = &pChampT; asgnSchemaIntg(pSchemaIntg); ajouteVectInfoTransformation(aVectInfoTransformation, pChampGeo, pSchemaIntg); ajouteVectInfoInterpolation (aVectInfoInterpolationNU, pChampU.reqFonctionsInterpolation(), aVectInfoTransformation, pChampGeo, pSchemaIntg); ajouteVectInfoInterpolation (aVectInfoInterpolationNV, pChampV.reqFonctionsInterpolation(), aVectInfoTransformation, pChampGeo, pSchemaIntg); ajouteVectInfoChamp (aVectInfoChampU, pChampU, aVectInfoTransformation, pSchemaIntg); ajouteVectInfoChamp (aVectInfoChampUPrecedent, pChampUPrecedent, aVectInfoTransformation, pSchemaIntg); ajouteVectInfoChamp (aVectInfoChampPseudoSigmaPrecedent, pChampPseudoSigmaO4Precedent, aVectInfoTransformation, pSchemaIntg); ajouteVectInfoChamp (aVectInfoChampS, pChampS, aVectInfoTransformation, pSchemaIntg); ajouteVectInfoChamp (aVectInfoChampT, pChampT, aVectInfoTransformation, pSchemaIntg); ajouteVectInfoChamp (aVectInfoChampI, pChampI, aVectInfoTransformation, pSchemaIntg); lMsg = asgnCouplage(pChampV, pChampU); _GIREF_TEST_INVARIANTS; return lMsg; } // ******************************************************************** // // Description: Ajoute aux champs les évaluations demandées pour le calcul // // Entrée: bool pAssembleMatrice: Booléen qui nous indique si on assemble la matrice // bool pAssembleResidu : Booléen qui nous indique si on assemble le résidu // // Sortie: Aucune // // ******************************************************************** void GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin::ajouteEvaluationsAFaire(bool pAssembleMatrice,bool pAssembleResidu) { _GIREF_TEST_INVARIANTS; Entier lEvaluationsAFaireChampGeo = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireNV = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireU = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireUPrecedent = EvalsChampLocal::Aucune; Entier lEvaluationsAFairePseudoSigmaO4 = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireS = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireT = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireI = EvalsChampLocal::Aucune; // Dans tous les cas on a lEvaluationsAFaireChampGeo |= EvalsChampLocal::EvalueDetJxPoid; lEvaluationsAFaireNV |= EvalsChampLocal::EvalueDeriveesPremieresFctBase; // Si on assemble la matrice, les évaluations suivantes seront nécessaires. if (pAssembleMatrice) { lEvaluationsAFaireS |= EvalsChampLocal::EvalueChamp; //lEvaluationsAFaireT |= EvalsChampLocal::EvalueChamp; lEvaluationsAFaireU |= EvalsChampLocal::EvalueDeriveesPremieresFctBase; } // Si on assemble le résidu, les évaluations suivantes seront nécessaires. if (pAssembleResidu) { lEvaluationsAFaireU |= EvalsChampLocal::EvalueDeriveesPremieresChamp; lEvaluationsAFaireUPrecedent |= EvalsChampLocal::EvalueDeriveesPremieresChamp; lEvaluationsAFairePseudoSigmaO4 |= EvalsChampLocal::EvalueChamp; lEvaluationsAFaireT |= EvalsChampLocal::EvalueChamp; lEvaluationsAFaireI |= EvalsChampLocal::EvalueChamp; } reqChampGeometrique().ajouteEvaluationsTransformationAFaire(lEvaluationsAFaireChampGeo, reqSchemaIntg()); aChampU -> ajouteEvaluationsAFaire(lEvaluationsAFaireU, reqSchemaIntg()); aChampV -> ajouteEvaluationsAFaire(lEvaluationsAFaireNV, reqSchemaIntg()); aChampUPrecedent -> ajouteEvaluationsAFaire(lEvaluationsAFaireUPrecedent, reqSchemaIntg()); aChampPseudoSigmaO4Precedent -> ajouteEvaluationsAFaire(lEvaluationsAFairePseudoSigmaO4, reqSchemaIntg()); aChampS -> ajouteEvaluationsAFaire(lEvaluationsAFaireS, reqSchemaIntg()); aChampT -> ajouteEvaluationsAFaire(lEvaluationsAFaireT, reqSchemaIntg()); aChampI -> ajouteEvaluationsAFaire(lEvaluationsAFaireI, reqSchemaIntg()); _GIREF_TEST_INVARIANTS; } // ******************************************************************** // // Description: Calcul élémentaire et assemblage pour un élément quelconque. // // Entrée: BlocMElem& pBlocMElemRigidite: La matrice rigidité élémentaire // BlocVElem& pBlocVElemResidu: Le vecteur résidu élémentaire. // // Sortie: Aucun // // ******************************************************************** void GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin::evalueIntegrale(BlocMElem& pBlocMElem, BlocVElem& pBlocVElem) { _GIREF_TEST_INVARIANTS; const bool lAssembleMatrice = reqAssembleMatrice(); const bool lAssembleResidu = reqAssembleResidu(); //La matrice élémentaire est-elle symétrique? //Si oui, on ne calcule que la moitié inférieure que l'on recopie dans la moitié supérieure. const bool lSymetrie = aChampV == aChampU; //Boucle sur les points d'intégration const Entier lNbPoints = aVectInfoTransformation->reqLongueur(); for (Entier k=0; k& ldNVdX = lLocalNV.reqDeriveeXFctBase(); const VectDyn& ldNVdY = lLocalNV.reqDeriveeYFctBase(); const VectDyn& ldNVdZ = lLocalNV.reqDeriveeZFctBase(); const Entier lNbFctBaseV = ldNVdX.reqLongueur(); const DReel& lDetJxPoid = (*aVectInfoTransformation)[k].reqDetJxPoid(); _GIREF_ASSERTION(lDetJxPoid != 0.0, "Le Poids*Jacobien de l'element est 0.0"); const TenseurO4Sym& lTensS = (*aVectInfoChampS)[k].reqValeurChamp(); const DReel lS0000 = lDetJxPoid*lTensS.reqComposante0000(); const DReel lS1111 = lDetJxPoid*lTensS.reqComposante1111(); const DReel lS2222 = lDetJxPoid*lTensS.reqComposante2222(); const DReel lS1000 = lDetJxPoid*lTensS.reqComposante1000(); const DReel lS2000 = lDetJxPoid*lTensS.reqComposante2000(); const DReel lS1010 = lDetJxPoid*lTensS.reqComposante1010(); const DReel lS2010 = lDetJxPoid*lTensS.reqComposante2010(); const DReel lS2020 = lDetJxPoid*lTensS.reqComposante2020(); const DReel lS1100 = lDetJxPoid*lTensS.reqComposante1100(); const DReel lS1110 = lDetJxPoid*lTensS.reqComposante1110(); const DReel lS2011 = lDetJxPoid*lTensS.reqComposante2011(); const DReel lS2100 = lDetJxPoid*lTensS.reqComposante2100(); const DReel lS2110 = lDetJxPoid*lTensS.reqComposante2110(); const DReel lS2120 = lDetJxPoid*lTensS.reqComposante2120(); const DReel lS2200 = lDetJxPoid*lTensS.reqComposante2200(); const DReel lS2210 = lDetJxPoid*lTensS.reqComposante2210(); const DReel lS2220 = lDetJxPoid*lTensS.reqComposante2220(); const DReel lS2111 = lDetJxPoid*lTensS.reqComposante2111(); const DReel lS2121 = lDetJxPoid*lTensS.reqComposante2121(); const DReel lS2211 = lDetJxPoid*lTensS.reqComposante2211(); const DReel lS2221 = lDetJxPoid*lTensS.reqComposante2221(); const TenseurO4Sym& lTensT = (*aVectInfoChampT)[k].reqValeurChamp(); const DReel lT0000 = lDetJxPoid*lTensT.reqComposante0000(); const DReel lT1111 = lDetJxPoid*lTensT.reqComposante1111(); const DReel lT2222 = lDetJxPoid*lTensT.reqComposante2222(); const DReel lT1000 = lDetJxPoid*lTensT.reqComposante1000(); const DReel lT2000 = lDetJxPoid*lTensT.reqComposante2000(); const DReel lT1010 = lDetJxPoid*lTensT.reqComposante1010(); const DReel lT2010 = lDetJxPoid*lTensT.reqComposante2010(); const DReel lT2020 = lDetJxPoid*lTensT.reqComposante2020(); const DReel lT1100 = lDetJxPoid*lTensT.reqComposante1100(); const DReel lT1110 = lDetJxPoid*lTensT.reqComposante1110(); const DReel lT2011 = lDetJxPoid*lTensT.reqComposante2011(); const DReel lT2100 = lDetJxPoid*lTensT.reqComposante2100(); const DReel lT2110 = lDetJxPoid*lTensT.reqComposante2110(); const DReel lT2120 = lDetJxPoid*lTensT.reqComposante2120(); const DReel lT2200 = lDetJxPoid*lTensT.reqComposante2200(); const DReel lT2210 = lDetJxPoid*lTensT.reqComposante2210(); const DReel lT2220 = lDetJxPoid*lTensT.reqComposante2220(); const DReel lT2111 = lDetJxPoid*lTensT.reqComposante2111(); const DReel lT2121 = lDetJxPoid*lTensT.reqComposante2121(); const DReel lT2211 = lDetJxPoid*lTensT.reqComposante2211(); const DReel lT2221 = lDetJxPoid*lTensT.reqComposante2221(); TenseurO2Sym lContraintesPrecedentes(0.); if(lAssembleResidu) { const InfoLocalesChamp& lLocalU = (*aVectInfoChampU)[k]; const Vecteur3D& ldUdX = lLocalU.reqValeurDeriveeXChamp(); const Vecteur3D& ldUdY = lLocalU.reqValeurDeriveeYChamp(); const Vecteur3D& ldUdZ = lLocalU.reqValeurDeriveeZChamp(); DReel lEpsilonXX = ldUdX.reqX(); DReel lEpsilonYY = ldUdY.reqY(); DReel lEpsilonZZ = ldUdZ.reqZ(); DReel lGammaXZ = ldUdX.reqZ() + ldUdZ.reqX(); DReel lGammaXY = ldUdX.reqY() + ldUdY.reqX(); DReel lGammaYZ = ldUdY.reqZ() + ldUdZ.reqY(); // // On calcule S:epsilon(U) // DReel lSigmaXZ = (lS2000*lEpsilonXX+lS2011*lEpsilonYY+lS2220*lEpsilonZZ+lGammaYZ*lS2120+lGammaXZ*lS2020+lGammaXY*lS2010); DReel lSigmaXY = (lS1110*lEpsilonYY+lS1000*lEpsilonXX+lS2210*lEpsilonZZ+lGammaYZ*lS2110+lGammaXZ*lS2010+lGammaXY*lS1010); DReel lSigmaYZ = (lS2100*lEpsilonXX+lS2221*lEpsilonZZ+lS2111*lEpsilonYY+lGammaYZ*lS2121+lGammaXZ*lS2120+lGammaXY*lS2110); DReel lSigmaXX = (lS2200*lEpsilonZZ+lS0000*lEpsilonXX+lS1100*lEpsilonYY+lGammaYZ*lS2100+lGammaXZ*lS2000+lGammaXY*lS1000); DReel lSigmaYY = (lS2211*lEpsilonZZ+lS1111*lEpsilonYY+lS1100*lEpsilonXX+lGammaYZ*lS2111+lGammaXZ*lS2011+lGammaXY*lS1110); DReel lSigmaZZ = (lS2211*lEpsilonYY+lS2222*lEpsilonZZ+lS2200*lEpsilonXX+lGammaYZ*lS2221+lGammaXZ*lS2220+lGammaXY*lS2210); const InfoLocalesChamp& lLocalUPrec = (*aVectInfoChampUPrecedent)[k]; const Vecteur3D& ldUPrecdX = lLocalUPrec.reqValeurDeriveeXChamp(); const Vecteur3D& ldUPrecdY = lLocalUPrec.reqValeurDeriveeYChamp(); const Vecteur3D& ldUPrecdZ = lLocalUPrec.reqValeurDeriveeZChamp(); DReel lEpsilonPrecXX = ldUPrecdX.reqX(); DReel lEpsilonPrecYY = ldUPrecdY.reqY(); DReel lEpsilonPrecZZ = ldUPrecdZ.reqZ(); DReel lGammaPrecXZ = ldUPrecdX.reqZ() + ldUPrecdZ.reqX(); DReel lGammaPrecXY = ldUPrecdX.reqY() + ldUPrecdY.reqX(); DReel lGammaPrecYZ = ldUPrecdY.reqZ() + ldUPrecdZ.reqY(); // // On soustrait T:epsilon(UPrecedent) // lSigmaXZ -= (lT2000*lEpsilonPrecXX+lT2011*lEpsilonPrecYY+lT2220*lEpsilonPrecZZ+lGammaPrecYZ*lT2120+lGammaPrecXZ*lT2020+lGammaPrecXY*lT2010); lSigmaXY -= (lT1110*lEpsilonPrecYY+lT1000*lEpsilonPrecXX+lT2210*lEpsilonPrecZZ+lGammaPrecYZ*lT2110+lGammaPrecXZ*lT2010+lGammaPrecXY*lT1010); lSigmaYZ -= (lT2100*lEpsilonPrecXX+lT2221*lEpsilonPrecZZ+lT2111*lEpsilonPrecYY+lGammaPrecYZ*lT2121+lGammaPrecXZ*lT2120+lGammaPrecXY*lT2110); lSigmaXX -= (lT2200*lEpsilonPrecZZ+lT0000*lEpsilonPrecXX+lT1100*lEpsilonPrecYY+lGammaPrecYZ*lT2100+lGammaPrecXZ*lT2000+lGammaPrecXY*lT1000); lSigmaYY -= (lT2211*lEpsilonPrecZZ+lT1111*lEpsilonPrecYY+lT1100*lEpsilonPrecXX+lGammaPrecYZ*lT2111+lGammaPrecXZ*lT2011+lGammaPrecXY*lT1110); lSigmaZZ -= (lT2211*lEpsilonPrecYY+lT2222*lEpsilonPrecZZ+lT2200*lEpsilonPrecXX+lGammaPrecYZ*lT2221+lGammaPrecXZ*lT2220+lGammaPrecXY*lT2210); lContraintesPrecedentes.asgnComposantes(lSigmaXX,lSigmaYY,lSigmaZZ,lSigmaXY,lSigmaYZ,lSigmaXZ); const TenseurO4Sym& lTenseurI = (*aVectInfoChampI)[k].reqValeurChamp(); const TenseurO4SymMin& lPseudoSigma = (*aVectInfoChampPseudoSigmaPrecedent)[k].reqValeurChamp(); TenseurO2Sym lTenseurPoids(lDetJxPoid); TenseurO4SymMin lIFoisSigma(produitSansSommation(lTenseurI,lPseudoSigma)); lContraintesPrecedentes += lIFoisSigma.contractionDoubleIndicesDeGauche(lTenseurPoids); } const InfoLocalesInterpolation& lLocalNU = (*aVectInfoInterpolationNU)[k]; const VectDyn& ldNUdX = lLocalNU.reqDeriveeXFctBase(); const VectDyn& ldNUdY = lLocalNU.reqDeriveeYFctBase(); const VectDyn& ldNUdZ = lLocalNU.reqDeriveeZFctBase(); const Entier lNbDDLsU = ldNUdX.reqLongueur(); Entier i2 = 0; for(Entier i=0; i