// ************************************************************************ // --- 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_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin.cc // // ******************************************************************** #include "GeorgesVisco_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin.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_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin::GeorgesVisco_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin() :TFGenerique(), aChampU(0), aChampV(0), aChampSigmaPrecedent(0), aChampS(0), aChampIndicateurEtat(0), aVectInfoTransformation(0), aVectInfoInterpolationNU(0), aVectInfoInterpolationNV(0), aVectInfoChampU(0), aVectInfoChampUPrecedent(0), aVectInfoChampSigmaPrecedent(0), aVectInfoChampS(0), aVectInfoChampIndicateurEtat(0) { _GIREF_CONSTRUCTEUR_INVARIANTS; _GIREF_TEST_INVARIANTS; } // ******************************************************************** // // Description: Constructeur par copie // // Entrée: pTF le terme de formulation à copier // // Sortie: Aucune // // ******************************************************************** GeorgesVisco_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin::GeorgesVisco_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin(const GeorgesVisco_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin& 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_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin::~GeorgesVisco_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin() { _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_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin& GeorgesVisco_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin::operator= (const GeorgesVisco_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin& 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 // ChampTensorielO2Sym& pChampSigmaPrecedent: Le champ contenant les deplacements au pas précédent // ChampTensorielO2Sym& pChampSigmaPrecedent: Le champ contenant le tenseur des contraintes au pas précédent pour la composante élastique // ChampTensorielO4Sym& pChampS: Le champ fonction de la valeur courante du tenseur d'élasticité // ChampScalaire& pChampIndicateurEtat: Le champ indiquant l'état du solide // 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_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin::asgnChamps(ChampGeometrique& pChampGeo, ChampVectoriel3D& pChampU, ChampVectoriel3D& pChampV, ChampVectoriel3D& pChampUPrecedent, ChampTensorielO2Sym& pChampSigmaPrecedent, ChampTensorielO4Sym& pChampS, ChampScalaire& pChampIndicateurEtat, const SchemaIntg& pSchemaIntg) { _GIREF_TEST_INVARIANTS; ERMsg lMsg = ERMsg(ERMsg::OK); asgnChampGeometrique(pChampGeo); aChampU = &pChampU; aChampV = &pChampV; aChampUPrecedent = &pChampUPrecedent; aChampSigmaPrecedent = &pChampSigmaPrecedent; aChampS = &pChampS; aChampIndicateurEtat = &pChampIndicateurEtat; 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 (aVectInfoChampS, pChampS, aVectInfoTransformation, pSchemaIntg); ajouteVectInfoChamp (aVectInfoChampUPrecedent, pChampUPrecedent, aVectInfoTransformation, pSchemaIntg); ajouteVectInfoChamp (aVectInfoChampSigmaPrecedent, pChampSigmaPrecedent, aVectInfoTransformation, pSchemaIntg); ajouteVectInfoChamp (aVectInfoChampIndicateurEtat, pChampIndicateurEtat, 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_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin::ajouteEvaluationsAFaire(bool pAssembleMatrice,bool pAssembleResidu) { _GIREF_TEST_INVARIANTS; _GIREF_PRECONDITION(0 != aChampU, "Les champs ne sont pas assignés"); Entier lEvaluationsAFaireChampGeo = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireNV = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireU = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireUPrecedent = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireSigmaPrecedent = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireS = EvalsChampLocal::Aucune; Entier lEvaluationsAFaireIndicateurEtat = EvalsChampLocal::Aucune; // Dans tous les cas on a lEvaluationsAFaireChampGeo |= EvalsChampLocal::EvalueDetJxPoid; lEvaluationsAFaireNV |= EvalsChampLocal::EvalueDeriveesPremieresFctBase; lEvaluationsAFaireS |= EvalsChampLocal::EvalueChamp; // Si on assemble la matrice, les évaluations suivantes seront nécessaires. if (pAssembleMatrice) { lEvaluationsAFaireU |= EvalsChampLocal::EvalueDeriveesPremieresFctBase; } // Si on assemble le résidu, les évaluations suivantes seront nécessaires. if (pAssembleResidu) { lEvaluationsAFaireU |= EvalsChampLocal::EvalueDeriveesPremieresChamp; lEvaluationsAFaireUPrecedent |= EvalsChampLocal::EvalueDeriveesPremieresChamp; lEvaluationsAFaireSigmaPrecedent |= EvalsChampLocal::EvalueChamp; lEvaluationsAFaireIndicateurEtat |= EvalsChampLocal::EvalueChamp; } reqChampGeometrique().ajouteEvaluationsTransformationAFaire(lEvaluationsAFaireChampGeo, reqSchemaIntg()); aChampU -> ajouteEvaluationsAFaire(lEvaluationsAFaireU, reqSchemaIntg()); aChampV -> ajouteEvaluationsAFaire(lEvaluationsAFaireNV, reqSchemaIntg()); aChampUPrecedent -> ajouteEvaluationsAFaire(lEvaluationsAFaireUPrecedent, reqSchemaIntg()); aChampSigmaPrecedent -> ajouteEvaluationsAFaire(lEvaluationsAFaireSigmaPrecedent, reqSchemaIntg()); aChampS -> ajouteEvaluationsAFaire(lEvaluationsAFaireS, reqSchemaIntg()); aChampIndicateurEtat -> ajouteEvaluationsAFaire(lEvaluationsAFaireIndicateurEtat, 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_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin::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 InfoLocalesInterpolation& lLocalNU = (*aVectInfoInterpolationNU)[k]; const VectDyn& ldNUdX = lLocalNU.reqDeriveeXFctBase(); const VectDyn& ldNUdY = lLocalNU.reqDeriveeYFctBase(); const VectDyn& ldNUdZ = lLocalNU.reqDeriveeZFctBase(); const Entier lNbDDLsU = ldNUdX.reqLongueur(); bool lUtiliseHooke = ((*aVectInfoChampIndicateurEtat)[k].reqValeurChamp() <= 0.); 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); if(lUtiliseHooke){ lContraintesPrecedentes.asgnComposantes(lSigmaXX,lSigmaYY,lSigmaZZ,lSigmaXY,lSigmaYZ,lSigmaXZ); } else { const InfoLocalesChamp& lLocalUPrec = (*aVectInfoChampUPrecedent)[k]; const Vecteur3D& ldUPrecdX = lLocalUPrec.reqValeurDeriveeXChamp(); const Vecteur3D& ldUPrecdY = lLocalUPrec.reqValeurDeriveeYChamp(); const Vecteur3D& ldUPrecdZ = lLocalUPrec.reqValeurDeriveeZChamp(); lEpsilonXX = ldUPrecdX.reqX(); lEpsilonYY = ldUPrecdY.reqY(); lEpsilonZZ = ldUPrecdZ.reqZ(); lGammaXZ = ldUPrecdX.reqZ() + ldUPrecdZ.reqX(); lGammaXY = ldUPrecdX.reqY() + ldUPrecdY.reqX(); lGammaYZ = ldUPrecdY.reqZ() + ldUPrecdZ.reqY(); // // On soustrait S:epsilon(UPrecedent) // lSigmaXZ -= (lS2000*lEpsilonXX+lS2011*lEpsilonYY+lS2220*lEpsilonZZ+lGammaYZ*lS2120+lGammaXZ*lS2020+lGammaXY*lS2010); lSigmaXY -= (lS1110*lEpsilonYY+lS1000*lEpsilonXX+lS2210*lEpsilonZZ+lGammaYZ*lS2110+lGammaXZ*lS2010+lGammaXY*lS1010); lSigmaYZ -= (lS2100*lEpsilonXX+lS2221*lEpsilonZZ+lS2111*lEpsilonYY+lGammaYZ*lS2121+lGammaXZ*lS2120+lGammaXY*lS2110); lSigmaXX -= (lS2200*lEpsilonZZ+lS0000*lEpsilonXX+lS1100*lEpsilonYY+lGammaYZ*lS2100+lGammaXZ*lS2000+lGammaXY*lS1000); lSigmaYY -= (lS2211*lEpsilonZZ+lS1111*lEpsilonYY+lS1100*lEpsilonXX+lGammaYZ*lS2111+lGammaXZ*lS2011+lGammaXY*lS1110); lSigmaZZ -= (lS2211*lEpsilonYY+lS2222*lEpsilonZZ+lS2200*lEpsilonXX+lGammaYZ*lS2221+lGammaXZ*lS2220+lGammaXY*lS2210); lContraintesPrecedentes.asgnComposantes(lSigmaXX,lSigmaYY,lSigmaZZ,lSigmaXY,lSigmaYZ,lSigmaXZ); lContraintesPrecedentes += (*aVectInfoChampSigmaPrecedent)[k].reqValeurChamp(); } } Entier i2 = 0; for(Entier i=0; i