Changeset 4699


Ignore:
Timestamp:
07/21/10 14:47:33 (15 years ago)
Author:
Mathieu Morlighem
Message:

oops... reverted to old version

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.cpp

    r4698 r4699  
    1 s:
     1/*!\file RiftConstraints.cpp
     2 * \brief: manage penalties for rifts
     3 */
     4
     5#include "./PenaltyConstraintsLocal.h"
     6#include "../../EnumDefinitions/EnumDefinitions.h"
     7#include "../../include/include.h"
     8#include "../../shared/shared.h"
     9
     10#define _ZIGZAGCOUNTER_
     11
     12int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,Loads* loads,int min_mechanical_constraints,int analysis_type){
     13
     14        int num_unstable_constraints=0;
     15        int converged=0;
     16        int potential;
     17        ISSMERROR(" check analysis_type across all routine!");
     18
     19        Constrain(&num_unstable_constraints,loads);
     20        if(num_unstable_constraints==0)converged=1;
     21       
     22       
     23        if(IsFrozen(loads)){
     24                converged=1;
     25                num_unstable_constraints=0;
     26        }
     27        else if(num_unstable_constraints<=min_mechanical_constraints){
     28                _printf_("   freezing constraints\n");
     29                FreezeConstraints(loads);
     30        }
     31
     32        /*Assign output pointers: */
     33        *pconverged=converged;
     34        *pnum_unstable_constraints=num_unstable_constraints;
     35}
     36
     37int IsMaterialStable(Loads* loads){
     38
     39        int i;
     40       
     41        Riftfront* riftfront=NULL;
     42        int found=0;
     43        int mpi_found=0;
     44
     45        /*go though loads, and if non-linearity of the material has converged, let all penalties know: */
     46        for (i=0;i<loads->Size();i++){
     47
     48                if(RiftfrontEnum==loads->GetEnum(i)){
     49
     50                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     51
     52                        if (riftfront->IsMaterialStable()){
     53                                found=1;
     54                                /*do not break! all penalties should get informed the non-linearity converged!*/
     55                        }
     56                }
     57        }
     58
     59        #ifdef _PARALLEL_
     60        MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
     61        MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);               
     62        found=mpi_found;
     63        #endif
     64
     65        return found;
     66}
     67
     68int RiftIsPresent(Loads* loads,int analysis_type){
     69
     70
     71        int i;
     72       
     73        int found=0;
     74        int mpi_found=0;
     75
     76        /*go though loads, and figure out if one of the loads is a Riftfront: */
     77        for (i=0;i<loads->Size();i++){
     78                Load* load=(Load*)loads->GetObjectByOffset(i);
     79                if(load->InAnalysis(analysis_type)){
     80                        if(RiftfrontEnum==loads->GetEnum(i)){
     81                                found=1;
     82                                break;
     83                        }
     84                }
     85        }
     86
     87        #ifdef _PARALLEL_
     88        MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
     89        MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);               
     90        found=mpi_found;
     91        #endif
     92
     93        return found;
     94}
     95
     96int IsPreStable(Loads* loads){
     97
     98
     99        int i;
     100       
     101        Riftfront* riftfront=NULL;
     102        int found=0;
     103        int mpi_found=0;
     104
     105        /*go though loads, and figure out if one of the penpair loads is still not stable: */
     106        for (i=0;i<loads->Size();i++){
     107
     108                if(RiftfrontEnum==loads->GetEnum(i)){
     109
     110                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     111
     112                        if (riftfront->PreStable()==0){
     113                                found=1;
     114                                break;
     115                        }
     116                }
     117        }
     118
     119        #ifdef _PARALLEL_
     120        MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
     121        MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);               
     122        found=mpi_found;
     123        #endif
     124
     125        if (found){
     126                /*We found an unstable constraint. : */
     127                return 0;
     128        }
     129        else{
     130                return 1;
     131        }
     132}
     133
     134int SetPreStable(Loads* loads){
     135
     136
     137        int i;
     138       
     139        Riftfront* riftfront=NULL;
     140        int found=0;
     141        int mpi_found=0;
     142
     143        /*go though loads, and set loads to pre stable.:*/
     144        for (i=0;i<loads->Size();i++){
     145
     146                if(RiftfrontEnum==loads->GetEnum(i)){
     147
     148                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     149                        riftfront->SetPreStable();
     150                }
     151        }
     152}
     153
     154int PreConstrain(int* pnum_unstable_constraints,Loads* loads){
     155
     156        int                     i;
     157       
     158        /* generic object pointer: */
     159        Riftfront* riftfront=NULL;
     160
     161        int unstable;
     162        int sum_num_unstable_constraints;
     163        int num_unstable_constraints=0;
     164               
     165        /*Enforce constraints: */
     166        for (i=0;i<loads->Size();i++){
     167
     168                if (RiftfrontEnum==loads->GetEnum(i)){
     169
     170                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     171
     172                        riftfront->PreConstrain(&unstable);
     173
     174                        num_unstable_constraints+=unstable;
     175                }
     176        }
     177
     178        #ifdef _PARALLEL_
     179        MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
     180        MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);               
     181        num_unstable_constraints=sum_num_unstable_constraints;
     182        #endif
     183       
     184        /*Assign output pointers: */
     185        *pnum_unstable_constraints=num_unstable_constraints;
     186
     187}
     188
     189int Constrain(int* pnum_unstable_constraints,Loads* loads){
     190
     191        int                     i;
     192       
     193        /* generic object pointer: */
     194        Riftfront* riftfront=NULL;
     195
     196        int unstable;
     197        int sum_num_unstable_constraints;
     198        int num_unstable_constraints=0;
     199               
     200        /*Enforce constraints: */
     201        for (i=0;i<loads->Size();i++){
     202
     203                if (RiftfrontEnum==loads->GetEnum(i)){
     204
     205                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     206
     207                        riftfront->Constrain(&unstable);
     208
     209                        num_unstable_constraints+=unstable;
     210                }
     211        }
     212
     213        #ifdef _PARALLEL_
     214        MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
     215        MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);               
     216        num_unstable_constraints=sum_num_unstable_constraints;
     217        #endif
     218       
     219        /*Assign output pointers: */
     220        *pnum_unstable_constraints=num_unstable_constraints;
     221
     222}
     223
     224void FreezeConstraints(Loads* loads){
     225
     226        int                     i;
     227       
     228        /* generic object pointer: */
     229        Riftfront* riftfront=NULL;
     230
     231        /*Enforce constraints: */
     232        for (i=0;i<loads->Size();i++){
     233
     234                if (RiftfrontEnum==loads->GetEnum(i)){
     235
     236                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     237
     238                        riftfront->FreezeConstraints();
     239
     240                }
     241        }
     242
     243}
     244
     245int IsFrozen(Loads* loads){
     246
     247        int                     i;
     248       
     249        /* generic object pointer: */
     250        Riftfront* riftfront=NULL;
     251        int found=0;
     252        int mpi_found=0;
     253
     254        /*Enforce constraints: */
     255        for (i=0;i<loads->Size();i++){
     256
     257                if (RiftfrontEnum==loads->GetEnum(i)){
     258
     259                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     260                        if (riftfront->IsFrozen()){
     261                                found=1;
     262                                break;
     263                        }
     264                }
     265        }
     266       
     267        /*Is there just one found? that would mean we have frozen! : */
     268        #ifdef _PARALLEL_
     269        MPI_Reduce (&found,&mpi_found,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD );
     270        MPI_Bcast(&mpi_found,1,MPI_DOUBLE,0,MPI_COMM_WORLD);               
     271        found=mpi_found;
     272        #endif
     273
     274        return found;
     275}
     276
     277int MaxPenetrationInInputs(Loads* loads){
     278
     279        int                     i;
     280       
     281        /* generic object pointer: */
     282        Riftfront* riftfront=NULL;
     283
     284        /*rift penetration: */
     285        double max_penetration=0;
     286        double mpi_max_penetration;
     287        double penetration;
     288
     289        /*Ok, we are going to find the grid pairs which are not penetrating, even though they
     290         * are penalised. We will release only the one with has least <0 penetration. : */
     291
     292        max_penetration=0;
     293        for (i=0;i<loads->Size();i++){
     294
     295                if (RiftfrontEnum==loads->GetEnum(i)){
     296
     297                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     298
     299                        riftfront->MaxPenetration(&penetration);
     300
     301                        if (penetration>max_penetration)max_penetration=penetration;
     302                }
     303        }
     304
     305        #ifdef _PARALLEL_
     306        MPI_Reduce (&max_penetration,&mpi_max_penetration,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD );
     307        MPI_Bcast(&mpi_max_penetration,1,MPI_DOUBLE,0,MPI_COMM_WORLD);               
     308        max_penetration=mpi_max_penetration;
     309        #endif
     310
     311        /*feed max_penetration to inputs: */
     312        for(i=0;i<loads->Size();i++){
     313                Load* load=(Load*)loads->GetObjectByOffset(i);
     314                load->InputUpdateFromVector(&max_penetration,MaxPenetrationEnum,ConstantEnum);
     315        }
     316}
     317
     318int PotentialUnstableConstraints(Loads* loads){
     319
     320        int                     i;
     321       
     322        /* generic object pointer: */
     323        Riftfront* riftfront=NULL;
     324
     325        /*Ok, we are going to find the grid pairs which are not penetrating, even though they
     326         * are penalised. We will release only the one with has least <0 penetration. : */
     327        int unstable=0;
     328        int sum_num_unstable_constraints=0;
     329        int num_unstable_constraints=0;
     330
     331        for (i=0;i<loads->Size();i++){
     332
     333                if (RiftfrontEnum==loads->GetEnum(i)){
     334
     335                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     336
     337                        riftfront->PotentialUnstableConstraint(&unstable);
     338
     339                        num_unstable_constraints+=unstable;
     340                }
     341        }
     342
     343        #ifdef _PARALLEL_
     344        MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
     345        MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);               
     346        num_unstable_constraints=sum_num_unstable_constraints;
     347        #endif
     348
     349        return num_unstable_constraints;
     350}
Note: See TracChangeset for help on using the changeset viewer.