geometrie/src/mg_volume.cpp

//------------------------------------------------------------
//------------------------------------------------------------
//  MAGiC
//  Jean Christophe Cuilli?re et Vincent FRANCOIS
//  D?partement de G?nie M?canique - UQTR
//------------------------------------------------------------
//  Le projet MAGIC est un projet de recherche du d?partement
//  de g?nie m?canique de l'Universit? du Qu?bec ?
//  Trois Rivi?res
//  Les librairies ne peuvent ?tre utilis?es sans l'accord
//  des auteurs (contact : francois@uqtr.ca)
//------------------------------------------------------------
//------------------------------------------------------------
//
//       mg_volume.cpp
//
//------------------------------------------------------------
//------------------------------------------------------------
//  COPYRIGHT 2000
//  Version du 02/03/2006 ? 11H22
//------------------------------------------------------------
//------------------------------------------------------------


#include "gestionversion.h"
#include "mg_volume.h"
#include "ot_mathematique.h"
#include "mg_maillage.h"
#include "vct_volume.h"
//#include "message.h"

MG_VOLUME::MG_VOLUME(std::string idori,unsigned long num):MG_ELEMENT_TOPOLOGIQUE(num,idori),num_materiau(-1),vect(NULL)
{
}

MG_VOLUME::MG_VOLUME(std::string idori):MG_ELEMENT_TOPOLOGIQUE(idori),num_materiau(-1),vect(NULL)
{
}

MG_VOLUME::MG_VOLUME(MG_VOLUME& mdd):MG_ELEMENT_TOPOLOGIQUE(mdd),lst_coquille(mdd.lst_coquille),vect(mdd.vect)
{
}

MG_VOLUME::~MG_VOLUME()
{
}

void MG_VOLUME::ajouter_mg_coquille(class MG_COQUILLE* mgcoq)
{
    lst_coquille.insert(lst_coquille.end(),mgcoq);
}

void MG_VOLUME::supprimer_mg_coquille(class MG_COQUILLE* mgcoq)
{
    std::vector<MG_COQUILLE*>::iterator i;
    for (i=lst_coquille.begin();i!=lst_coquille.end();i++)
    {
        if ((*i)==mgcoq)
        {
            lst_coquille.erase(i);
            return;
        }
    }
}


int MG_VOLUME::get_nb_mg_coquille(void)
{
    return lst_coquille.size();
}

MG_COQUILLE* MG_VOLUME::get_mg_coquille(int num)
{
    return lst_coquille[num];
}

int MG_VOLUME::get_dimension(void)
{
    return 3;
}
VCT& MG_VOLUME::get_vectorisation(void)
{
//VCT* p=NULL;
//return *p;
    if (vect==NULL) vect=new VCT_VOLUME(this);
    return *vect;
}
void MG_VOLUME::get_topologie_sousjacente(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lst)
{
    int nbcoq=lst_coquille.size();
    for (int i=0;i<nbcoq;i++)
    {
        MG_COQUILLE* coq=lst_coquille[i];
        int nbface=coq->get_nb_mg_coface();
        for (int j=0;j<nbface;j++)
        {
            MG_FACE* face=coq->get_mg_coface(j)->get_face();
            lst->ajouter(face);
            face->get_topologie_sousjacente(lst);
        }
    }
}


BOITE_3D MG_VOLUME::get_boite_3D(int pas_echantillon)
{
BOITE_3D boitevol;
int nb_coquille=get_nb_mg_coquille();
for (int i=0;i<nb_coquille;i++)
    {
    MG_COQUILLE* coq=get_mg_coquille(i);
    int nb_face=coq->get_nb_mg_coface();
    for (int j=0;j<nb_face;j++)
          {
          MG_FACE* face=coq->get_mg_coface(j)->get_face();
          BOITE_3D tmp=face->get_boite_3D(pas_echantillon);
          if ((i==0) && (j==0))  boitevol=tmp;
          else boitevol=tmp+boitevol;
          }
    }
return boitevol;
}

void MG_VOLUME::enregistrer(std::ostream& o,double version)
{
    o << "%" << get_id() << "=VOLUME(" << get_idoriginal() << ",(";
    for (unsigned int i=0;i<lst_coquille.size();i++)
    {
        o << "$" << lst_coquille[i]->get_id();
        if (i!=lst_coquille.size()-1) o << ",";
        else o << ")";
    }
    int nb=get_nb_ccf();
    if (version<2) 
    {
    o << "," << num_materiau << "," << nb;
    if (nb!=0)
    {
        o << ",(";
        for (int i=0;i<nb;i++)
        {
            char nom[3];
            get_type_ccf(i,nom);
            o << "(" <<  nom << "," << get_valeur_ccf(i) << ")";
            if (i!=nb-1) o << "," ;
        }
        o << ")";
    }
    o << ");" << std::endl;
    }
    else
    {
    o << "," ;
    enregistrer_ccf(o,version);  
    o << ");" << std::endl;
    }
}

void MG_VOLUME::change_num_materiau(int num)
{
    num_materiau=num;
}

int MG_VOLUME::get_num_materiau(void)
{
    return num_materiau;
}


bool MG_VOLUME::est_un_volume_element(void)
{
return false;
}
void MG_VOLUME::get_propriete_massique(class MG_MAILLAGE* mai,double& volume,class OT_VECTEUR_3D& cdm,class OT_MATRICE_3D& inertieglobale,class OT_MATRICE_3D& inertiecdm,double dens)
{
    double vol=0.;
    double xg=0.;
    double yg=0.;
    double zg=0.;
    double a=0.;
    double b=0.;
    double c=0.;
    double d=0.;
    double e=0.;
    double f=0.;
    int nbcoquille=get_nb_mg_coquille();
    for (int i=0;i<nbcoquille;i++)
    {
        MG_COQUILLE* coq=get_mg_coquille(i);
        int nbface=coq->get_nb_mg_coface();
        for (int j=0;j<nbface;j++)
        {
            MG_COFACE* coface=coq->get_mg_coface(j);
            MG_FACE* face=coface->get_face();
            TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
            for (MG_TRIANGLE* tri=(MG_TRIANGLE*)face->get_lien_maillage()->get_premier(it);tri!=NULL;tri=(MG_TRIANGLE*)face->get_lien_maillage()->get_suivant(it))
            {
                if (mai->get_mg_triangleid(tri->get_id())==NULL) continue;
                MG_NOEUD* noeud1=tri->get_noeud1();
                MG_NOEUD* noeud2=tri->get_noeud2();
                MG_NOEUD* noeud3=tri->get_noeud3();
                double *xyz1=noeud1->get_coord();
                double *xyz2=noeud2->get_coord();
                double *xyz3=noeud3->get_coord();
                OT_VECTEUR_3D vec1(xyz1,xyz2);
                OT_VECTEUR_3D vec2(xyz1,xyz3);
                OT_VECTEUR_3D n=vec1&vec2;
                double detj=n.get_longueur();
                n.norme();
                n=-1.*coface->get_orientation()*n;
                for (int k=0;k<3;k++)
                {
                    double xsi=1./6.;
                    double eta=1./6.;
                    double wi=1./6.;
                    if (k==1) xsi=2./3.;
                    if (k==2) eta=2./3.;
                    double x=(1-xsi-eta)*xyz1[0]+xsi*xyz2[0]+eta*xyz3[0];
                    double y=(1-xsi-eta)*xyz1[1]+xsi*xyz2[1]+eta*xyz3[1];
                    double z=(1-xsi-eta)*xyz1[2]+xsi*xyz2[2]+eta*xyz3[2];
                    OT_VECTEUR_3D psi1(x,0,0);
                    OT_VECTEUR_3D psi2(0.5*x*x,0,0);
                    OT_VECTEUR_3D psi3(0,0.5*y*y,0);
                    OT_VECTEUR_3D psi4(0,0,0.5*z*z);
                    OT_VECTEUR_3D psi5(0,y*y*y/3.,z*z*z/3.);
                    OT_VECTEUR_3D psi6(x*x*x/3.,0,z*z*z/3.);
                    OT_VECTEUR_3D psi7(x*x*x/3.,y*y*y/3.,0);
                    OT_VECTEUR_3D psi8(0.,y*y*z/2.,0);
                    OT_VECTEUR_3D psi9(x*x*z/2.,0.,0);
                    OT_VECTEUR_3D psi10(x*x*y/2.,0,0);
                    vol=vol+detj*wi*(psi1*n);
                    xg=xg+wi*detj*(psi2*n);
                    yg=yg+wi*detj*(psi3*n);
                    zg=zg+wi*detj*(psi4*n);
                    a=a+wi*detj*(psi5*n);
                    b=b+wi*detj*(psi6*n);
                    c=c+wi*detj*(psi7*n);
                    d=d+wi*detj*(psi8*n);
                    e=e+wi*detj*(psi9*n);
                    f=f+wi*detj*(psi10*n);
                }
            }
        }
    }
    volume=vol;
    cdm.change_x(xg/volume);
    cdm.change_y(yg/volume);
    cdm.change_z(zg/volume);
    inertieglobale(0,0)=a*dens;
    inertieglobale(1,1)=b*dens;
    inertieglobale(2,2)=c*dens;
    inertieglobale(1,0)=f*dens;
    inertieglobale(0,1)=f*dens;
    inertieglobale(2,0)=e*dens;
    inertieglobale(0,2)=e*dens;
    inertieglobale(2,1)=d*dens;
    inertieglobale(1,2)=d*dens;
    OT_VECTEUR_3D m1(volume*(cdm.get_y()*cdm.get_y()+cdm.get_z()*cdm.get_z()),volume*cdm.get_x()*cdm.get_y(),volume*cdm.get_x()*cdm.get_z());
    OT_VECTEUR_3D m2(volume*cdm.get_x()*cdm.get_y(),volume*(cdm.get_x()*cdm.get_x()+cdm.get_z()*cdm.get_z()),volume*cdm.get_y()*cdm.get_z());
    OT_VECTEUR_3D m3(volume*cdm.get_x()*cdm.get_z(),volume*cdm.get_y()*cdm.get_z(),volume*(cdm.get_x()*cdm.get_x()+cdm.get_y()*cdm.get_y()));
    m1=(-1.)*m1;
    m2=(-1.)*m2;
    m3=(-1.)*m3;
    OT_MATRICE_3D change(m1,m2,m3);
    inertiecdm=inertieglobale+change;
}

Revision:	763
Committed:	Wed Dec 2 19:55:53 2015 UTC (9 years, 5 months ago) by francois
File size:	8223 byte(s)
Log Message:	Le fichier MAGiC est maintenant versionné. LA version actuelle est 2.0. L'ancienne version est 1.0. Tout est transparent pour l'utilisateur. Les vieilles versions sont lisibles mais les nouveaux enregistrements sont dans la version la plus récente. Changement des conditions aux limites : ajout d'un parametre pour dire si la condition numerique est une valeur ou une formule ou un lien vers une autre entité magic. Les parametres pour saisir sont maintenant -ccf -ccfi -ccff -ccft -ccfit -ccfft
#	User	Rev	Content
1	francois	283	//------------------------------------------------------------
2			//------------------------------------------------------------
3			// MAGiC
4			// Jean Christophe Cuilli?re et Vincent FRANCOIS
5			// D?partement de G?nie M?canique - UQTR
6			//------------------------------------------------------------
7			// Le projet MAGIC est un projet de recherche du d?partement
8			// de g?nie m?canique de l'Universit? du Qu?bec ?
9			// Trois Rivi?res
10			// Les librairies ne peuvent ?tre utilis?es sans l'accord
11			// des auteurs (contact : francois@uqtr.ca)
12			//------------------------------------------------------------
13			//------------------------------------------------------------
14			//
15			// mg_volume.cpp
16			//
17			//------------------------------------------------------------
18			//------------------------------------------------------------
19			// COPYRIGHT 2000
20			// Version du 02/03/2006 ? 11H22
21			//------------------------------------------------------------
22			//------------------------------------------------------------
23
24
25			#include "gestionversion.h"
26			#include "mg_volume.h"
27			#include "ot_mathematique.h"
28			#include "mg_maillage.h"
29			#include "vct_volume.h"
30			//#include "message.h"
31
32			MG_VOLUME::MG_VOLUME(std::string idori,unsigned long num):MG_ELEMENT_TOPOLOGIQUE(num,idori),num_materiau(-1),vect(NULL)
33			{
34			}
35
36			MG_VOLUME::MG_VOLUME(std::string idori):MG_ELEMENT_TOPOLOGIQUE(idori),num_materiau(-1),vect(NULL)
37			{
38			}
39
40			MG_VOLUME::MG_VOLUME(MG_VOLUME& mdd):MG_ELEMENT_TOPOLOGIQUE(mdd),lst_coquille(mdd.lst_coquille),vect(mdd.vect)
41			{
42			}
43
44			MG_VOLUME::~MG_VOLUME()
45			{
46			}
47
48			void MG_VOLUME::ajouter_mg_coquille(class MG_COQUILLE* mgcoq)
49			{
50			lst_coquille.insert(lst_coquille.end(),mgcoq);
51			}
52
53			void MG_VOLUME::supprimer_mg_coquille(class MG_COQUILLE* mgcoq)
54			{
55			std::vector<MG_COQUILLE*>::iterator i;
56			for (i=lst_coquille.begin();i!=lst_coquille.end();i++)
57			{
58			if ((*i)==mgcoq)
59			{
60			lst_coquille.erase(i);
61			return;
62			}
63			}
64			}
65
66
67			int MG_VOLUME::get_nb_mg_coquille(void)
68			{
69			return lst_coquille.size();
70			}
71
72			MG_COQUILLE* MG_VOLUME::get_mg_coquille(int num)
73			{
74			return lst_coquille[num];
75			}
76
77			int MG_VOLUME::get_dimension(void)
78			{
79			return 3;
80			}
81			VCT& MG_VOLUME::get_vectorisation(void)
82			{
83			//VCT* p=NULL;
84			//return *p;
85			if (vect==NULL) vect=new VCT_VOLUME(this);
86			return *vect;
87			}
88			void MG_VOLUME::get_topologie_sousjacente(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE> lst)
89			{
90			int nbcoq=lst_coquille.size();
91			for (int i=0;i<nbcoq;i++)
92			{
93			MG_COQUILLE* coq=lst_coquille[i];
94			int nbface=coq->get_nb_mg_coface();
95			for (int j=0;j<nbface;j++)
96			{
97			MG_FACE* face=coq->get_mg_coface(j)->get_face();
98			lst->ajouter(face);
99			face->get_topologie_sousjacente(lst);
100			}
101			}
102			}
103
104
105	francois	632
106
107			BOITE_3D MG_VOLUME::get_boite_3D(int pas_echantillon)
108			{
109			BOITE_3D boitevol;
110			int nb_coquille=get_nb_mg_coquille();
111			for (int i=0;i<nb_coquille;i++)
112			{
113			MG_COQUILLE* coq=get_mg_coquille(i);
114			int nb_face=coq->get_nb_mg_coface();
115			for (int j=0;j<nb_face;j++)
116			{
117			MG_FACE* face=coq->get_mg_coface(j)->get_face();
118			BOITE_3D tmp=face->get_boite_3D(pas_echantillon);
119			if ((i==0) && (j==0)) boitevol=tmp;
120			else boitevol=tmp+boitevol;
121			}
122			}
123			return boitevol;
124			}
125
126	francois	763	void MG_VOLUME::enregistrer(std::ostream& o,double version)
127	francois	283	{
128			o << "%" << get_id() << "=VOLUME(" << get_idoriginal() << ",(";
129			for (unsigned int i=0;i<lst_coquille.size();i++)
130			{
131			o << "$" << lst_coquille[i]->get_id();
132			if (i!=lst_coquille.size()-1) o << ",";
133			else o << ")";
134			}
135			int nb=get_nb_ccf();
136	francois	763	if (version<2)
137			{
138	francois	283	o << "," << num_materiau << "," << nb;
139			if (nb!=0)
140			{
141			o << ",(";
142			for (int i=0;i<nb;i++)
143			{
144			char nom[3];
145			get_type_ccf(i,nom);
146			o << "(" << nom << "," << get_valeur_ccf(i) << ")";
147			if (i!=nb-1) o << "," ;
148			}
149			o << ")";
150			}
151			o << ");" << std::endl;
152	francois	763	}
153			else
154			{
155			o << "," ;
156			enregistrer_ccf(o,version);
157			o << ");" << std::endl;
158			}
159	francois	283	}
160
161			void MG_VOLUME::change_num_materiau(int num)
162			{
163			num_materiau=num;
164			}
165
166			int MG_VOLUME::get_num_materiau(void)
167			{
168			return num_materiau;
169			}
170
171	francois	576
172			bool MG_VOLUME::est_un_volume_element(void)
173			{
174			return false;
175			}
176	francois	283	void MG_VOLUME::get_propriete_massique(class MG_MAILLAGE* mai,double& volume,class OT_VECTEUR_3D& cdm,class OT_MATRICE_3D& inertieglobale,class OT_MATRICE_3D& inertiecdm,double dens)
177			{
178			double vol=0.;
179			double xg=0.;
180			double yg=0.;
181			double zg=0.;
182			double a=0.;
183			double b=0.;
184			double c=0.;
185			double d=0.;
186			double e=0.;
187			double f=0.;
188			int nbcoquille=get_nb_mg_coquille();
189			for (int i=0;i<nbcoquille;i++)
190			{
191			MG_COQUILLE* coq=get_mg_coquille(i);
192			int nbface=coq->get_nb_mg_coface();
193			for (int j=0;j<nbface;j++)
194			{
195			MG_COFACE* coface=coq->get_mg_coface(j);
196			MG_FACE* face=coface->get_face();
197			TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
198			for (MG_TRIANGLE* tri=(MG_TRIANGLE)face->get_lien_maillage()->get_premier(it);tri!=NULL;tri=(MG_TRIANGLE)face->get_lien_maillage()->get_suivant(it))
199			{
200			if (mai->get_mg_triangleid(tri->get_id())==NULL) continue;
201			MG_NOEUD* noeud1=tri->get_noeud1();
202			MG_NOEUD* noeud2=tri->get_noeud2();
203			MG_NOEUD* noeud3=tri->get_noeud3();
204			double *xyz1=noeud1->get_coord();
205			double *xyz2=noeud2->get_coord();
206			double *xyz3=noeud3->get_coord();
207			OT_VECTEUR_3D vec1(xyz1,xyz2);
208			OT_VECTEUR_3D vec2(xyz1,xyz3);
209			OT_VECTEUR_3D n=vec1&vec2;
210			double detj=n.get_longueur();
211			n.norme();
212			n=-1.coface->get_orientation()n;
213			for (int k=0;k<3;k++)
214			{
215			double xsi=1./6.;
216			double eta=1./6.;
217			double wi=1./6.;
218			if (k==1) xsi=2./3.;
219			if (k==2) eta=2./3.;
220			double x=(1-xsi-eta)xyz1[0]+xsixyz2[0]+eta*xyz3[0];
221			double y=(1-xsi-eta)xyz1[1]+xsixyz2[1]+eta*xyz3[1];
222			double z=(1-xsi-eta)xyz1[2]+xsixyz2[2]+eta*xyz3[2];
223			OT_VECTEUR_3D psi1(x,0,0);
224			OT_VECTEUR_3D psi2(0.5xx,0,0);
225			OT_VECTEUR_3D psi3(0,0.5yy,0);
226			OT_VECTEUR_3D psi4(0,0,0.5zz);
227			OT_VECTEUR_3D psi5(0,yyy/3.,zzz/3.);
228			OT_VECTEUR_3D psi6(xxx/3.,0,zzz/3.);
229			OT_VECTEUR_3D psi7(xxx/3.,yyy/3.,0);
230			OT_VECTEUR_3D psi8(0.,yyz/2.,0);
231			OT_VECTEUR_3D psi9(xxz/2.,0.,0);
232			OT_VECTEUR_3D psi10(xxy/2.,0,0);
233			vol=vol+detjwi(psi1*n);
234			xg=xg+widetj(psi2*n);
235			yg=yg+widetj(psi3*n);
236			zg=zg+widetj(psi4*n);
237			a=a+widetj(psi5*n);
238			b=b+widetj(psi6*n);
239			c=c+widetj(psi7*n);
240			d=d+widetj(psi8*n);
241			e=e+widetj(psi9*n);
242			f=f+widetj(psi10*n);
243			}
244			}
245			}
246			}
247			volume=vol;
248			cdm.change_x(xg/volume);
249			cdm.change_y(yg/volume);
250			cdm.change_z(zg/volume);
251			inertieglobale(0,0)=a*dens;
252			inertieglobale(1,1)=b*dens;
253			inertieglobale(2,2)=c*dens;
254			inertieglobale(1,0)=f*dens;
255			inertieglobale(0,1)=f*dens;
256			inertieglobale(2,0)=e*dens;
257			inertieglobale(0,2)=e*dens;
258			inertieglobale(2,1)=d*dens;
259			inertieglobale(1,2)=d*dens;
260			OT_VECTEUR_3D m1(volume(cdm.get_y()cdm.get_y()+cdm.get_z()cdm.get_z()),volumecdm.get_x()cdm.get_y(),volumecdm.get_x()*cdm.get_z());
261			OT_VECTEUR_3D m2(volumecdm.get_x()cdm.get_y(),volume(cdm.get_x()cdm.get_x()+cdm.get_z()cdm.get_z()),volumecdm.get_y()*cdm.get_z());
262			OT_VECTEUR_3D m3(volumecdm.get_x()cdm.get_z(),volumecdm.get_y()cdm.get_z(),volume(cdm.get_x()cdm.get_x()+cdm.get_y()*cdm.get_y()));
263			m1=(-1.)*m1;
264			m2=(-1.)*m2;
265			m3=(-1.)*m3;
266			OT_MATRICE_3D change(m1,m2,m3);
267			inertiecdm=inertieglobale+change;
268			}
269