aster/src/mgopt_simp.cpp

#include "gestionversion.h"
#include "mgopt_simp.h"
#include "mg_file.h"
#include "mg_export.h"
#include <string.h>
#include <math.h>
#include "fct_taille_fem_solution.h"
#include "fct_generateur_3d.h"
#include "ot_cpu.h"
#include "mglanceuraster.h"

class SIMP_TETRA
{
public:
FEM_ELEMENT3* tet;
double densite;
double new_densite;
int design;
double new_denergie;
double energie;
double denergie;
double volume;
int maille_niveau;
int indice;
double distance_ref;
double distance_ref2;
double centre[3];
BOITE_3D get_boite_3D(void) 
        {
        return tet->get_boite_3D();
        };
std::vector<SIMP_TETRA*> voisin;
std::vector<SIMP_TETRA*> voisin2;
unsigned long get_id(void) {return tet->get_id();};
double distance(SIMP_TETRA* tet2) 
        {
        double dx=centre[0]-tet2->centre[0];
        double dy=centre[1]-tet2->centre[1];
        double dz=centre[2]-tet2->centre[2];
        return sqrt(dx*dx+dy*dy+dz*dz);
        };
};


MGOPT_SIMP::MGOPT_SIMP(bool save):MGOPT(save)
{
params.ajouter("f",0.3,OT_PARAMETRES::DOUBLE,"Fraction volumique de la methode SIMP");
params.ajouter("rminc",0.,OT_PARAMETRES::DOUBLE,"0. Valeur de la carte de taille (si pas de maillage le parametre fichiercarte permet d'utiliser une carte de taille) sinon valeur de rmin pour le lissage de la compliance");
params.ajouter("rmind",0.,OT_PARAMETRES::DOUBLE,"0. Valeur de la carte de taille (si pas de maillage le parametre fichiercarte permet d'utiliser une carte de taille) sinon valeur de rmin pour le lissage de la densite");
params.ajouter("coefvoisinc",1.0,OT_PARAMETRES::DOUBLE,"Coefficient de multiplication de rminc");
params.ajouter("coefvoisind",1.0,OT_PARAMETRES::DOUBLE,"Coefficient de multiplication de rmind");
params.ajouter("nbrniveau",150.,OT_PARAMETRES::DOUBLE,"Nombre de matériaux utilisés");
params.ajouter("p",3.,OT_PARAMETRES::DOUBLE,"Coefficient de penalite du module d'Young");
params.ajouter("ro_void",0.001,OT_PARAMETRES::DOUBLE,"Densite minimale consideree comme nulle");
params.ajouter("m",0.2,OT_PARAMETRES::DOUBLE,"Limite d'evolution de densité entre deux itérations");
params.ajouter("eta",0.5,OT_PARAMETRES::DOUBLE,"Coefficent d'affectation du beta");
params.ajouter("type_lissage",1.,OT_PARAMETRES::DOUBLE,"0. Pas de lissage 1. Lissage de la compliance 2. Lissage de la densité 3. Lissage de la compliance et de la densite");
params.ajouter("lissage_densite",1.,OT_PARAMETRES::DOUBLE,"0. Complet 1. Derniere iteration");
params.ajouter("type_lissage_densite",1.,OT_PARAMETRES::DOUBLE,"0. Distance 1. Distance pondéré 2. Gaussien 3. Gaussien pondéré");
params.ajouter("type_lissage_compliance",1.,OT_PARAMETRES::DOUBLE,"0. Sigmund 1997 1. Sigmund 2007");
params.ajouter("kc",1.0,OT_PARAMETRES::DOUBLE,"Ponderation de la distance dans le lissage de la compliance");
params.ajouter("kd",1.0,OT_PARAMETRES::DOUBLE,"Ponderation de la distance dans le lissage de la densité");
params.ajouter("convergence_lagrange",0.1,OT_PARAMETRES::DOUBLE,"Critere de convergence de la recherche du multiplicateur de Lagrange en \%");
params.ajouter("iter_max",50.,OT_PARAMETRES::DOUBLE,"Nombre d'iteration maximale dans le processus SIMP");
params.ajouter("critere_convergence",0.5,OT_PARAMETRES::DOUBLE,"Critere de convergence de la méthode SIMP en \%");
params.ajouter("ro_min",0.8,OT_PARAMETRES::DOUBLE,"Seuil de conservation des éléments");
params.ajouter("iter_vue",0.,OT_PARAMETRES::DOUBLE,"generation d'un fichier gmsh tous les iter_vue. 0  pas de generation");
params.ajouter("type_filtrage",0.,OT_PARAMETRES::DOUBLE,"0. filtrage par distance 1. filtrage par couches");
params.ajouter("nb_couchesc",1.,OT_PARAMETRES::DOUBLE,"nombre de couches filtrage de la compliance");
params.ajouter("nb_couchesd",1.,OT_PARAMETRES::DOUBLE,"nombre de couches filtrage de la densite");
}

MGOPT_SIMP::MGOPT_SIMP(MGOPT_SIMP &mdd):MGOPT(mdd)
{
}


MGOPT_SIMP::~MGOPT_SIMP()
{
}


void MGOPT_SIMP::optimisation(int num_adapt,FEM_MAILLAGE* fem,char *nomparamaster,int iter)
{
MG_GEOMETRIE* geo=fem->get_mg_geometrie();
affiche((char*)"     Initialisation");
std::vector<SIMP_TETRA*> lsttet;
double f=params.get_valeur("f");
double rminc=params.get_valeur("rminc");
double rmind=params.get_valeur("rmind");
if (((rminc<1e-16)||(rmind<1e-16)) && (carte==NULL))
        {
        affiche((char*)"     Lecture de la carte de taille");
        int typecarte=(int)params.get_valeur("typecarte");
  std::string fichiercarte=params.get_nom("fichiercarte");
  if ((typecarte==0) &&(typecarte==1)) 
  {
    carte=new FCT_GENERATEUR_3D<4>;
    carte->lire((char *)fichiercarte.c_str());
  }
  if ((typecarte==2) &&(typecarte==3)) carte=new FCT_TAILLE_FEM_SOLUTION((char *)fichiercarte.c_str());
  }
double volume_tot;
double volume_design=0;
double volume_non_design=0;
double unite=fem->get_mg_maillage()->get_mg_geometrie()->get_valeur_unite();
LISTE_FEM_ELEMENT3::iterator it;
int ntet=0;
for (FEM_ELEMENT3* tet=fem->get_premier_element3(it);tet!=NULL;tet=fem->get_suivant_element3(it))
        {
        SIMP_TETRA* tet2=new SIMP_TETRA;
        tet2->tet=tet;
        tet2->tet->change_numero(ntet);
        ntet++;
        tet2->centre[0]=0.;
        tet2->centre[1]=0.;
        tet2->centre[2]=0.;
        int nbnoeud=tet->get_nb_fem_noeud();
        for (int i=0;i<nbnoeud;i++)
                {
                tet2->centre[0]=tet2->centre[0]+tet->get_fem_noeud(i)->get_x()*unite;
                tet2->centre[1]=tet2->centre[1]+tet->get_fem_noeud(i)->get_y()*unite;
                tet2->centre[2]=tet2->centre[2]+tet->get_fem_noeud(i)->get_z()*unite;
                }
        tet2->centre[0]=tet2->centre[0]/nbnoeud;                
        tet2->centre[1]=tet2->centre[1]/nbnoeud;                
        tet2->centre[2]=tet2->centre[2]/nbnoeud;                
        int lig,col;
        double jac[9],uv[3]={1./4.,1./4.,1./4.};
        tet2->volume=1./6.*tet->get_jacobien(jac,uv,lig,col,unite);
        tet2->distance_ref=rminc;
        tet2->distance_ref2=rmind;
        double xyz[3]={tet2->centre[0]/unite,tet2->centre[1]/unite,tet2->centre[2]/unite};
        double val[9];
        if (tet2->distance_ref<1e-16) 
                {
                carte->evaluer(xyz,val);
                tet2->distance_ref=1./sqrt(val[0]);
                }
        if (tet2->distance_ref2<1e-16) 
                {
                carte->evaluer(xyz,val);
                tet2->distance_ref2=1./sqrt(val[0]);
                }
        tet2->distance_ref=params.get_valeur("coefvoisinc")*tet2->distance_ref*unite;
        tet2->distance_ref2=params.get_valeur("coefvoisind")*tet2->distance_ref2*unite;
        if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=MAGIC::ORIGINE::IMPOSE) 
                {
                tet2->design=1;
                tet2->densite=f;
                volume_design=volume_design+tet2->volume;
                }
         else 
                {
                tet2->design=0;
                tet2->densite=1.;
                volume_non_design=volume_non_design+tet2->volume;
                }
        lsttet.insert(lsttet.end(),tet2);
        }
volume_tot=volume_design+volume_non_design;
double volumemoyen=volume_tot/fem->get_nb_fem_element3();
affiche((char*)"     Recherche de voisins");
int nbsimptet=lsttet.size();
if (params.get_valeur("type_filtrage")==0.)
        {
        for (int i=0;i<nbsimptet;i++)
                {
                SIMP_TETRA* tet=lsttet[i];
                if (tet->design==0) continue;
                ajouter_voisin_distance(i,tet,lsttet);
                }
        }
if (params.get_valeur("type_filtrage")==1.)
        {
        for (int i=0;i<nbsimptet;i++)
                {
                SIMP_TETRA* tet=lsttet[i];
                if (tet->design==0) continue;
                ajouter_voisin_couche(i,tet,lsttet,params.get_valeur("nb_couchesc"),params.get_valeur("nb_couchesd"));
                }
        }

int nbiteration=1;
int ok=0;
int niveaumax=(int)params.get_valeur("nbrniveau");
int itervue=(int)params.get_valeur("iter_vue");
double p=params.get_valeur("p");
std::vector<double> Ctotiter;
while (ok==0)
        {
        if (nbiteration>params.get_valeur("iter_max")) break;
        char message[255];
        if (num_adapt==0) sprintf(message,"     Iteration \033[1;32m%d\033[1;33m",nbiteration);
        else sprintf(message,"     Iteration \033[1;32m%d.%d\033[1;33m",num_adapt,nbiteration);
        affiche(message);
        std::vector<FEM_ELEMENT3*>  *lstniveau=new std::vector<FEM_ELEMENT3*>[niveaumax+1];
        double densmin = params.get_valeur("ro_void");
        for (int i=0;i<nbsimptet;i++)
                {
                SIMP_TETRA* tet=lsttet[i];
                tet->maille_niveau=(int)((tet->densite-densmin)*niveaumax/(1.-densmin))+1;
                if (tet->maille_niveau>niveaumax) tet->maille_niveau=niveaumax;
                lstniveau[tet->maille_niveau].insert(lstniveau[tet->maille_niveau].end(),tet->tet);
                }
        MG_EXPORT exp;
        //int version_aster=params.get_valeur("version_aster");
        exp.lire_params_aster(nomparamaster);
        double version=exp.aster(geo,fem,nometude,MAGIC::CALCUL_ASTER::OPTIMISATIONTOPOLOGIQUE,(char*)"00000001",p,niveaumax,lstniveau);
        delete [] lstniveau;
        char messagever[255];
        sprintf(messagever,"        Calcul aster \033[1;32m%.1lf\033[1;33m",version);
        affiche(messagever);
        MGLANCEURASTER aster;
        int code=aster.execute(param,nometude);
        if (code!=0)
                {
                char message[500];
                sprintf(message,"          Code de sortie aster : %d",code);
                affiche(message);
                }
        affiche((char*)"        Analyse resultat");
        recupere_energie(lsttet);
        double Ctot=0;
        for (int i=0;i<nbsimptet;i++)
                {
                SIMP_TETRA* tet=lsttet[i];
                tet->denergie= -p*2.*tet->energie/tet->densite;
                Ctot=Ctot+2.*tet->energie;
                }
        Ctotiter.insert(Ctotiter.end(),Ctot);
        sprintf(message,"             Compliance \033[1;31m%le\033[1;33m",Ctot);
        affiche(message);
        if (nbiteration!=1)
                {
                int nb=Ctotiter.size();
                double c1=Ctotiter[nb-2];
                double c2=Ctotiter[nb-1];
                double ecart=fabs((c2-c1)/c1)*100.;
                if (ecart<params.get_valeur("critere_convergence")) ok=1;
                sprintf(message,"             Ecart \033[1;31m%lf%%\033[1;33m",ecart);
                affiche(message);
                }
        // lissage compliance
        double kc=params.get_valeur("kc");
        int type_lissage =(int)params.get_valeur("type_lissage");
        int lissage_densite =(int)params.get_valeur("lissage_densite");
        int type_lissage_densite =(int)params.get_valeur("type_lissage_densite");
        int type_lissage_compliance=(int)params.get_valeur("type_lissage_compliance");
        for (int i=0;i<nbsimptet;i++)
                {
                SIMP_TETRA* tet=lsttet[i];
                if ((type_lissage==0) || (type_lissage==2))
                        tet->new_denergie=tet->denergie;
                if ((type_lissage==1) || (type_lissage==3))
                        {
                        if (tet->design == 1)
                                {
                                double hi= pow(tet->distance_ref,kc);
                                double hisensi = hi*tet->densite*tet->denergie;
                                if (type_lissage_compliance==1)
                                        {
                                        hi=hi/tet->volume;
                                        hisensi=hisensi/std::max(volumemoyen*0.00001,tet->volume);
                                        }
                                int nbvoisin=tet->voisin.size();
                                for (int j = 0 ; j<nbvoisin ; j++)
                                        {
                                        SIMP_TETRA* tet2=tet->voisin[j];
                                        if (tet2->design == 1)
                                                {
                                                double dist=tet->distance(tet2);
                                                double hj =fabs(pow(tet->distance_ref - dist,kc));
                                                if (type_lissage_compliance==0)
                                                        {
                                                        hisensi=hisensi+tet2->densite*hj*tet2->denergie;
                                                        hi=hi+hj;
                                                        }
                                                if (type_lissage_compliance==1)
                                                        {
                                                        hisensi=hisensi+tet2->densite*hj*tet2->denergie/std::max(volumemoyen,tet2->volume);
                                                        hi=hi+hj/tet2->volume;
                                                        }
                                                hisensi=hisensi+tet2->densite*hj*tet2->denergie;
                                                hi=hi+hj;
                                                }
                                        }
                                tet->new_denergie = hisensi/hi/std::max(tet->densite,densmin);
                                }
                        else tet->new_denergie=tet->denergie;
                        }
                }
        // application formule pas encore connue
        double l1= 1e-8 ;
        double l2= 1e8;
        double m=params.get_valeur("m");
        double eta=params.get_valeur("eta");
        double convergence_lagrange=params.get_valeur("convergence_lagrange");
        double critere_densite = 0.0;
        int oklagrange=0;
        int compteurlagrange=0;
        while (oklagrange==0)
                {
                compteurlagrange++;
                double lmid = 0.5*(l2+l1);
                critere_densite = 0.0;
                        for (int i=0;i<nbsimptet;i++)
                        {
                        SIMP_TETRA* tet=lsttet[i];
                        if (tet->design==1)
                                {
                                double x1 = tet->densite+m;
                                double beta=-tet->new_denergie/lmid/tet->volume;
                                double x2 = tet->densite*pow(beta,eta);
                                double x3 = std::min(x1,x2);
                                double x4 = 1.;
                                double x5 = std::min(x3,x4);
                                double x6 = tet->densite-m;
                                double x7 = std::max(x5,x6);
                                double x8 = densmin;
                                tet->new_densite = std::max(x7,x8);
                                critere_densite = critere_densite + tet->new_densite*tet->volume;
                                }
                        else
                                {
                                tet->new_densite=1.;
                                }
                        }
                if (critere_densite - f*volume_design  > 0.)
                        l1=lmid;
                else
                        l2=lmid;
                if (100.*fabs(critere_densite- f*volume_design )/(f*volume_design)<convergence_lagrange) oklagrange=1;
                if (compteurlagrange>500) oklagrange=2;
                }
        if (oklagrange==1) sprintf(message,"             Convergence multiplicateur lagrange \033[1;31m%le%%\033[1;33m",100.*(critere_densite- f*volume_design )/(f*volume_design));
        if (oklagrange==2) sprintf(message,"             Divergence multiplicateur lagrange \033[1;31m%le %le\033[1;33m",l1,l2);
        affiche(message);
        // lissage de densite.
        double kd=params.get_valeur("kd");
        for (int i=0;i<nbsimptet;i++)
                {
                SIMP_TETRA* tet=lsttet[i];
                if ((type_lissage==2) || (type_lissage==3))
                        {
                        if (((lissage_densite==1)&&(ok==1)) || (lissage_densite==0))
                                {
                                if (tet->design == 1)
                                        {
                                        double widensite=0.;
                                        double wi= 0.;
                                        wi=poid_lissage(0.,tet->distance_ref2,kd,tet->volume,type_lissage_densite);
                                        widensite = wi*tet->new_densite;
                                        int nbvoisin=tet->voisin2.size();
                                        for (int j = 0 ; j<nbvoisin ; j++)
                                                {
                                                SIMP_TETRA* tet2=tet->voisin2[j];
                                                if (tet2->design == 1)
                                                        {
                                                        double dist=tet->distance(tet2);
                                                        double wj=poid_lissage(dist,tet->distance_ref2,kd,tet2->volume,type_lissage_densite);
                                                        wi = wi+wj;
                                                        widensite = widensite + tet2->new_densite*wj;
                                                        }
                                                }
                                        tet->densite = widensite/wi;
                                        }
                                else tet->densite=tet->new_densite;
                                }
                        else tet->densite=tet->new_densite;
                        }
                else tet->densite=tet->new_densite;
                }
        if (itervue!=0)
          if (nbiteration%itervue==0)
              {
              affiche((char*)"        Sauvegarde résultat iteration");
              sprintf(message,"%s_densite_iter%d.soltmp",nometudesortie,nbiteration);
              //int nbsolution=gestd->get_nb_fem_solution();
              //for (int i=nbsolution;i>0;i--version 2 du for)
                    //gestd->supprimer_fem_solution_du_gestionnaire(i-1);
              char message2[50];
              static int step=0;
              sprintf(message2,"Iteration%%%d",step);
              step++;
              FEM_SOLUTION* solution=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_element3(),message2,MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT3);
              gestd->ajouter_fem_solution(solution);
              sprintf(message,"Densite");
              solution->change_legende(0,message); 
              for (int i=0;i<nbsimptet;i++)
                  {
                  SIMP_TETRA* tet=lsttet[i];
                  solution->ecrire(tet->densite,i,0);
                  }
              }
        nbiteration++;
        }
if (itervue!=0)
  {
  affiche((char*)"     Sauvegarde GMSH résultats iterations");
  MG_EXPORT exp;
  char nomfichier[500];
  sprintf(nomfichier,"%s_iterations",nometudesortie);
  exp.gmsh(fem,nomfichier);
  }
double seuil=params.get_valeur("ro_min");
affiche((char*)"     Ecriture des donnees finales");
char message[255];
if (iter==0) sprintf(message,"%s.compliance",nometudesortie);
else sprintf(message,"%s_iter%d.compliance",nometudesortie,iter);
FILE *out=fopen(message,"wt");
time_t heurefin=time(NULL);
struct tm  tfin = *localtime(&heurefin);
fprintf(out,"*****************************************************\n");
fprintf(out,"           Optimisateur de topologie\n");
fprintf(out,"                   ERICCA - UQTR\n");
fprintf(out,"*****************************************************\n");
fprintf(out,"    Etude        : %s\n",nometude);
fprintf(out,"    Date         : %02u-%02u-%04u\n", tdebut.tm_mday, tdebut.tm_mon+1, 1900 + tdebut.tm_year);
fprintf(out,"    Heure debut  : %02u:%02u:%02u\n", tdebut.tm_hour, tdebut.tm_min, tdebut.tm_sec);
fprintf(out,"    Heure fin    : %02u:%02u:%02u\n", tfin.tm_hour, tfin.tm_min, tfin.tm_sec);
fprintf(out,"    Parametres   : \n");
int nbparam=params.get_nb();
for (int i=0;i<nbparam;i++)
  {
  if (params.get_type(i)==OT_PARAMETRES::DOUBLE) fprintf(out,"                   %s = %lf\n",params.get_nom(i).c_str(),params.get_valeur(i));
  if (params.get_type(i)==OT_PARAMETRES::STRING) fprintf(out,"                   %s = %s\n",params.get_nom(i).c_str(),params.get_nom(params.get_nom(i).c_str()).c_str());
  }
fprintf(out,"*****************************************************\n");
fprintf(out,"         Compliance après chaque itération\n");
fprintf(out,"*****************************************************\n");
for (int i=0;i<Ctotiter.size();i++)
        fprintf(out,"%le\n",Ctotiter[i]);
fprintf(out,"*****************************************************\n");
fprintf(out,"  Variation de la compliance après chaque itération\n");
fprintf(out,"*****************************************************\n");
for (int i=1;i<Ctotiter.size();i++)
        fprintf(out,"%lf%%\n",(Ctotiter[i]-Ctotiter[i-1])*100./Ctotiter[i-1]);
double crit1=0.;
double crit2=0.;
double crit3=0.;
double critv1[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
double critv2[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
for (int i=0;i<nbsimptet;i++)
        {
        SIMP_TETRA* tet=lsttet[i];
        if (tet->design==1)
                {
                crit1=crit1+tet->densite*tet->volume;
                if (tet->densite>seuil)
                        {
                        crit2=crit2+tet->densite*tet->volume;
                        crit3=crit3+tet->volume;
                        }
                for (int j=1;j<10;j++)
                        if (tet->densite>0.1*j)
                                {
                                critv1[j]=critv1[j]+tet->densite*tet->volume;
                                critv2[j]=critv2[j]+tet->volume;
                                }
                }
        }
fprintf(out,"*****************************************************\n");
fprintf(out,"         Résultat\n");
fprintf(out,"*****************************************************\n");
fprintf(out,"Volume du design : %le \n",volume_design);
fprintf(out,"Objectif du volume de design : %le \n",volume_design*f);
fprintf(out,"Volume de design obtenu : %le \n",crit1);
fprintf(out,"Volume de design obtenu avec le seuil: %le \n",crit2);
fprintf(out,"Volume reel de design obtenu avec le seuil: %le \n\n",crit3);
fprintf(out," : Volume : Volume reel : Objectif\n");
fprintf(out," : %le : %le : %le\n",volume_design,volume_design,volume_design);
for (int j=1;j<10;j++)
        fprintf(out,"Volume de design obtenu avec le seuil de %.1f : %le : %le : %le \n",j*0.1,critv1[j],critv2[j],volume_design*f);
fprintf(out,"*************************************************************\n");
fclose(out);
LISTE_FEM_NOEUD::iterator itnoeud;
int nbfemnoeud=fem->get_nb_fem_noeud();
/*double *nume=new double[nbfemnoeud];
double *deno=new double[nbfemnoeud];
int cpt=0;
for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
        {
        nd->change_numero(cpt);
        nume[cpt]=0.;
        deno[cpt]=0.;
        cpt++;
        }*/
if (iter==0) sprintf(message,"%s_densite.sol",nometudesortie);
else sprintf(message,"%s_densite_iter%d.sol",nometudesortie,iter);
int nbsolution=gestd->get_nb_fem_solution();
for (int i=nbsolution;i>0;i--)
        {
          FEM_SOLUTION* sol=gestd->get_fem_solution(i-1);
          std::string nom=sol->get_nom_fichier();
          char message[500],extension[500];
          sprintf(message,"%s",nom.c_str());
          char *p=strrchr(message,'.');
          strncpy(extension,p,strlen(message)+message-p);
          if (strcmp(extension,".soltmp")==0)
               gestd->supprimer_fem_solution(i-1);
            else gestd->supprimer_fem_solution_du_gestionnaire(i-1);
        }


FEM_SOLUTION* solution=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_element3(),"Optimisation",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT3);
gestd->ajouter_fem_solution(solution);
solution->change_legende(0,"Densite"); 
for (int i=0;i<nbsimptet;i++)
        {
        SIMP_TETRA* tet=lsttet[i];
        if (tet->design==1) 
                if (tet->densite>seuil) 
                        ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(MAGIC::ORIGINE::OPTIMISE); 
                else
                        ((MG_TETRA*)tet->tet->get_mg_element_maillage())->change_origine(MAGIC::ORIGINE::MAILLEUR_AUTO); 
        tet->tet->change_solution(tet->densite);
        solution->ecrire(tet->densite,i,0);
        /*for (int j=0;j<tet->tet->get_nb_fem_noeud();j++)
                {
                FEM_NOEUD* noeud=tet->tet->get_fem_noeud(j);
                nume[noeud->get_numero()]=nume[noeud->get_numero()]+tet->volume*tet->densite;
                deno[noeud->get_numero()]=deno[noeud->get_numero()]+tet->volume;
                }       */
        }
/*sprintf(message,"%s_densite2.sol",nometudesortie);
FEM_SOLUTION* solution2=new FEM_SOLUTION(fem,1,message,fem->get_nb_fem_noeud(),"Optimisation",ENTITE_NOEUD);
gestd->ajouter_fem_solution(solution2);
solution2->change_legende(0,"Densite");
cpt=0;
for (FEM_NOEUD *nd=fem->get_premier_noeud(itnoeud);nd!=NULL;nd=fem->get_suivant_noeud(itnoeud))
        {
        solution2->ecrire(nume[cpt]/deno[cpt],cpt,0,);
        cpt++;
        }
delete [] deno;
delete [] nume;*/
int nb=lsttet.size();
for (int i=0;i<nb;i++) delete lsttet[i];
}

double MGOPT_SIMP::poid_lissage(double dist,double distref,double k,double volume,int type)
{
double wi;
if (type==0) wi=pow(distref-dist,k);
if (type==1) wi=pow(distref-dist,k)*volume;
if (type==2) wi=exp(-dist*dist/2./distref/distref/9.)/2./M_PI/(distref/3.);
if (type==3) wi=volume*exp(-dist*dist/2./distref/distref/9.)/2./M_PI/(distref/3.);
return fabs(wi);
        
}

void MGOPT_SIMP::adapte_resultat(char *nomgestd,char *nomparam)
{
if (nomparam!=NULL) lire_params(nomparam);
affiche((char*)"");
affiche((char*)"*************************");
affiche((char*)"Optimisation de topologie");
affiche((char*)"*************************");
affiche((char*)"");
affiche((char*)"");
affiche((char*)"Changement du seuil dans les resultats");
double seuil=params.get_valeur("seuil");
gestd=new MG_FILE(nomgestd);
FEM_MAILLAGE* fem=gestd->get_fem_maillage(0);
FEM_SOLUTION* solution=gestd->get_fem_solution(0);
solution->active_solution(0);
LISTE_FEM_ELEMENT3::iterator it;
for (FEM_ELEMENT3 *tet=fem->get_premier_element3(it);tet!=NULL;tet=fem->get_suivant_element3(it))
        {
        if (((MG_TETRA*)tet->get_mg_element_maillage())->get_origine()!=MAGIC::ORIGINE::IMPOSE) 
                if (tet->get_solution()>seuil) 
                        ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAGIC::ORIGINE::OPTIMISE); 
                else
                        ((MG_TETRA*)tet->get_mg_element_maillage())->change_origine(MAGIC::ORIGINE::MAILLEUR_AUTO); 
        
        }
affiche((char*)"Enregistrement");
gestd->enregistrer(nomgestd);
affiche((char*)"Enregistrement sous GMSH");
char *p=strchr(nomgestd,'.');
strncpy(nometude,nomgestd,p-nomgestd);
nometude[p-nomgestd]=0;
MG_EXPORT exp;
char nomfichier[500];
sprintf(nomfichier,"%s_mg",nometude);
exp.gmsh(fem->get_mg_maillage(),nomfichier);
sprintf(nomfichier,"%s_fem",nometude);
exp.gmsh(fem,nomfichier);
affiche((char*)"Fin");
}


void MGOPT_SIMP::recupere_energie(std::vector<class SIMP_TETRA*> lsttet)
{
char message[750];
sprintf(message,"%s.resu",nometude);
FILE* in=fopen(message,"rt");
int fin=0;
do
        {
        char* res=fgets(message,750,in);
        if (feof(in)) fin=1;
        char mot1[100]; 
        char mot2[100]; 
        char mot3[100]; 
        char mot4[100]; 
        char mot5[100]; 
        char mot6[100]; 
        char mot7[100]; 
        char mot8[100]; 
        char mot9[100]; 
        char mot10[100]; 
        int numlu=sscanf(message,"%s %s %s %s %s %s %s %s %s %s",mot1,mot2,mot3,mot4,mot5,mot6,mot7,mot8,mot9,mot10);
        
        
        if (numlu>2)
        if (strcmp(mot2,"TOTAL_JOB")==0)
          {
          double val;
          sscanf(mot8,"%lf",&val);  
          CODE_ASTER_CPU=CODE_ASTER_CPU+val;
          sscanf(mot10,"%lf",&val);  
          CODE_ASTER_ECOULE=CODE_ASTER_ECOULE+val;
          }
        
        
        if (numlu>9)
        if (strcmp(mot1,"CHAMP")==0)
         if (strcmp(mot2,"PAR")==0) 
          if (strcmp(mot3,"ELEMENT")==0) 
           if (strcmp(mot4,"CONSTANT")==0) 
            if (strcmp(mot5,"SUR")==0) 
             if (strcmp(mot6,"ELEMENT")==0) 
              if (strcmp(mot7,"DE")==0) 
               if (strcmp(mot8,"NOM")==0) 
                if (strcmp(mot9,"SYMBOLIQUE")==0) 
                 if (strcmp(mot10,"ENERGIE")==0) 
                        {
                        int fin2=0;
                        int passe=0;
                        int nbelement=0;
                        do 
                                {
                                char message[750];
                                char* res=fgets(message,750,in);
                                char mot1[500];
                                char mot2[500];
                                int numlu=sscanf(message,"%s %s",mot1,mot2);
                                int decalage;
                                if ((numlu==2) && (strcmp(mot2,"TOTALE")==0))
                                        {
                                        char *p=strchr(mot1,'M')+1;
                                        int num=atoi(p);
                                        if (passe==0) {passe=1;decalage=num;}
                                        char* res=fgets(message,750,in);
                                        double val;
                                        sscanf(message,"%lf",&val);
                                        lsttet[num-decalage]->energie=val;
                                        nbelement++;
                                        }
                                if (nbelement == lsttet.size()) {fin2=1;fin=0;}
                                }
                                
                        while (fin2==0);
                        }
        }
while (fin==0);
fclose(in);
}


void MGOPT_SIMP::ajouter_voisin_distance(int i,SIMP_TETRA* tet,std::vector<SIMP_TETRA*> &lst)
{
tet->indice=i;
int nbnoeud=tet->tet->get_nb_fem_noeud();
int correspondance[4];
if (nbnoeud==4) 
        {
        correspondance[0]=0;
        correspondance[1]=1;
        correspondance[2]=2;
        correspondance[3]=3;
        }
if (nbnoeud==10) 
        {
        correspondance[0]=0;
        correspondance[1]=2;
        correspondance[2]=4;
        correspondance[3]=9;
        }
SIMP_TETRA* tetcour=tet;
int ok=0;
int compteur=0;
std::vector<SIMP_TETRA*>& lstvoisin=tet->voisin;
if (tet->distance_ref<tet->distance_ref2)
        lstvoisin=tet->voisin2;
while (ok==0)
        {
        for (int j=0;j<4;j++)
                {
                int num=correspondance[j];
                FEM_NOEUD *noeud=tetcour->tet->get_fem_noeud(num);
                int nbtetra=noeud->get_lien_element3()->get_nb();
                for (int k=0;k<nbtetra;k++)
                        {
                        FEM_ELEMENT3* ftet=noeud->get_lien_element3()->get(k);
                        SIMP_TETRA* stet=lst[ftet->get_numero()];
                        if (stet->indice!=i)
                                {
                                stet->indice=i;
                                double dist=tet->distance(stet);
                                if (dist<tet->distance_ref) 
                                                tet->voisin.insert(tet->voisin.end(),stet);
                                if (dist<tet->distance_ref2) 
                                                tet->voisin2.insert(tet->voisin2.end(),stet);
                                }
                        }
                }
        if (compteur>=lstvoisin.size()) ok=1;
        else 
                {
                tetcour=lstvoisin[compteur];
                compteur++;
                }
                
        }

}


void MGOPT_SIMP::ajouter_voisin_couche(int i,SIMP_TETRA* tet,std::vector<SIMP_TETRA*> &lst,int nb_couches,int nb_couches2)
{
tet->indice=i;
int nbnoeud=tet->tet->get_nb_fem_noeud();
int correspondance[4];
if (nbnoeud==4) 
        {
        correspondance[0]=0;
        correspondance[1]=1;
        correspondance[2]=2;
        correspondance[3]=3;
        }
if (nbnoeud==10) 
        {
        correspondance[0]=0;
        correspondance[1]=2;
        correspondance[2]=4;
        correspondance[3]=9;
        }
std::vector<SIMP_TETRA*>& lstvoisin=tet->voisin;
if (tet->distance_ref<tet->distance_ref2)
        lstvoisin=tet->voisin2;
tet->voisin.insert(tet->voisin.end(),tet);
tet->voisin2.insert(tet->voisin2.end(),tet);
tet->indice=i;
int debut=0,fin=1;
for (int ok=0;ok<std::max(nb_couches,nb_couches2);ok++)
        {
        for (int iliste=debut;iliste<fin;iliste++)
          {
          SIMP_TETRA* tetcour;
         if (nb_couches>nb_couches2 ) tetcour=tet->voisin[iliste]; else tetcour=tet->voisin2[iliste];
          for (int j=0;j<4;j++)
                {
                int num=correspondance[j];
                FEM_NOEUD *noeud=tetcour->tet->get_fem_noeud(num);
                int nbtetra=noeud->get_lien_element3()->get_nb();
                for (int k=0;k<nbtetra;k++)
                        {
                        FEM_ELEMENT3* ftet=noeud->get_lien_element3()->get(k);
                        SIMP_TETRA* stet=lst[ftet->get_numero()];
                        if (stet->indice!=i)
                                {
                                stet->indice=i;
                                if (ok<nb_couches) tet->voisin.insert(tet->voisin.end(),stet);
                                if (ok<nb_couches2) tet->voisin2.insert(tet->voisin2.end(),stet);
                                }
                        }
                }
            }
         debut=fin;
         if (nb_couches>nb_couches2 ) fin=tet->voisin.size(); else fin=tet->voisin2.size();
        }
}
Revision:	951
Committed:	Fri Aug 10 15:17:17 2018 UTC (6 years, 9 months ago) by couturad
File size:	26704 byte(s)
Log Message:	-> Ajout de Project Chrono (voir CMakeLists.txt). -> Ajout d'un générateur de microstructure basé sur la dynamique des corps rigides (MSTRUCT_GENERATEUR_DCR). -> Ajout d'un opérateur de décallage de la topologie (MG_CG_OP_TRANSF_DECALLAGE). -> Retrait de «using namespace std» (conflit avec namespace chrono) et modification des fichiers affectés. -> Modification de mailleur2d.cpp afin d'enregistrer un fichier MAGiC (void.magic) lorsque le nombre d'itération dépasse la valeur maximale.