Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 8403) +++ /issm/trunk/src/c/Makefile.am (revision 8404) @@ -630,10 +630,10 @@ ./modules/SystemMatricesx/SystemMatricesx.cpp\ ./modules/SystemMatricesx/SystemMatricesx.h\ - ./modules/PenaltyConstraintsx/PenaltyConstraintsx.cpp\ - ./modules/PenaltyConstraintsx/PenaltyConstraintsx.h\ - ./modules/PenaltyConstraintsx/PenaltyConstraintsLocal.h\ - ./modules/PenaltyConstraintsx/RiftConstraints.cpp\ - ./modules/PenaltyConstraintsx/MeltingConstraints.cpp\ - ./modules/PenaltyConstraintsx/MeltingIsPresent.cpp\ + ./modules/ConstraintsStatex/ConstraintsStatex.cpp\ + ./modules/ConstraintsStatex/ConstraintsStatex.h\ + ./modules/ConstraintsStatex/ConstraintsStateLocal.h\ + ./modules/ConstraintsStatex/RiftConstraintsState.cpp\ + ./modules/ConstraintsStatex/ThermalConstraintsState.cpp\ + ./modules/ConstraintsStatex/MeltingIsPresent.cpp\ ./modules/Responsex/Responsex.h\ ./modules/Responsex/Responsex.cpp\ @@ -1265,10 +1265,10 @@ ./modules/SystemMatricesx/SystemMatricesx.cpp\ ./modules/SystemMatricesx/SystemMatricesx.h\ - ./modules/PenaltyConstraintsx/PenaltyConstraintsx.cpp\ - ./modules/PenaltyConstraintsx/PenaltyConstraintsx.h\ - ./modules/PenaltyConstraintsx/PenaltyConstraintsLocal.h\ - ./modules/PenaltyConstraintsx/RiftConstraints.cpp\ - ./modules/PenaltyConstraintsx/MeltingConstraints.cpp\ - ./modules/PenaltyConstraintsx/MeltingIsPresent.cpp\ + ./modules/ConstraintsStatex/ConstraintsStatex.cpp\ + ./modules/ConstraintsStatex/ConstraintsStatex.h\ + ./modules/ConstraintsStatex/ConstraintsStateLocal.h\ + ./modules/ConstraintsStatex/RiftConstraintsState.cpp\ + ./modules/ConstraintsStatex/ThermalConstraintsState.cpp\ + ./modules/ConstraintsStatex/MeltingIsPresent.cpp\ ./modules/Responsex/Responsex.h\ ./modules/Responsex/Responsex.cpp\ Index: /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h =================================================================== --- /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h (revision 8404) +++ /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h (revision 8404) @@ -0,0 +1,31 @@ +/*!\file: ConstraintsStateLocal.h + * \brief local header files + */ + +#ifndef _CONSTRAINTSSTATELOCAL_H +#define _CONSTRAINTSSTATELOCAL_H + +#include "../../Container/Container.h" +#include "../../objects/objects.h" + +/*melting: */ +void ThermalConstraintsState(Loads* loads, int* pconverged, int* pnum_unstable_constraints,int analysis_type); +int MeltingIsPresent(Loads* loads,int analysis_type); + +/*rifts module: */ +int RiftIsPresent(Loads* loads,int analysis_type); +void RiftConstraintsState(int* pconverged, int* pnum_unstable_constraints,Loads* loads,int min_mechanical_constraints,int analysis_type); +void RiftConstrain(int* pnum_unstable_constraints,Loads* loads,int analysis_type); +int RiftIsFrozen(Loads* loads,int analysis_type); +void RiftFreezeConstraints(Loads* loads,int analysis_type); + +/*rifts, trial and errors: */ +int RiftIsPreStable(Loads* loads); +void RiftSetPreStable(Loads* loads); +void RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads); +void RiftMaxPenetrationInInputs(Loads* loads); +int RiftPotentialUnstableConstraints(Loads* loads); +int RiftIsMaterialStable(Loads* loads); + +#endif /* _CONSTRAINTSSTATEX_H */ + Index: /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp =================================================================== --- /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp (revision 8404) +++ /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp (revision 8404) @@ -0,0 +1,49 @@ +/*!\file ConstraintsStatex + * \brief: set up penalty constraints on loads + */ + +#include "./ConstraintsStatex.h" +#include "./ConstraintsStateLocal.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void ConstraintsStatex(int* pconverged, int* pnum_unstable_constraints, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads,Materials* materials, Parameters* parameters){ + + int i; + + extern int num_procs; + extern int my_rank; + + /*output: */ + int converged=0; + int num_unstable_constraints=0; + int min_mechanical_constraints=0; + int analysis_type; + + /*Display message*/ + _printf_(VerboseModule()," Constraining penalties\n"); + + /*recover parameters: */ + parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Do we have penalties linked to rifts? In this case, run our special rifts penalty + * management routine, otherwise, skip : */ + if (RiftIsPresent(loads,analysis_type)){ + RiftConstraintsState(&converged,&num_unstable_constraints,loads,min_mechanical_constraints,analysis_type); + } + else if(MeltingIsPresent(loads,analysis_type)){ + ThermalConstraintsState(loads,&converged,&num_unstable_constraints,analysis_type); + } + else{ + /*Do nothing, no constraints management!:*/ + num_unstable_constraints=0; + converged=1; + } + + /*Assign output pointers: */ + *pconverged=converged; + *pnum_unstable_constraints=num_unstable_constraints; +} Index: /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStatex.h =================================================================== --- /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStatex.h (revision 8404) +++ /issm/trunk/src/c/modules/ConstraintsStatex/ConstraintsStatex.h (revision 8404) @@ -0,0 +1,15 @@ +/*!\file: ConstraintsStatex.h + * \brief header file for penalty constraints module + */ + +#ifndef _CONSTRAINTSSTATEX_H +#define _CONSTRAINTSSTATEX_H + +#include "../../Container/Container.h" +#include "../../objects/objects.h" + +/* local prototypes: */ +void ConstraintsStatex(int* pconverged, int* pnum_unstable_constraints, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads,Materials* materials, Parameters* parameters); + +#endif /* _CONSTRAINTSSTATEX_H */ + Index: /issm/trunk/src/c/modules/ConstraintsStatex/MeltingIsPresent.cpp =================================================================== --- /issm/trunk/src/c/modules/ConstraintsStatex/MeltingIsPresent.cpp (revision 8404) +++ /issm/trunk/src/c/modules/ConstraintsStatex/MeltingIsPresent.cpp (revision 8404) @@ -0,0 +1,37 @@ +/*!\file: MeltingIsPresent.cpp + * \brief + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./ConstraintsStateLocal.h" + +int MeltingIsPresent(Loads* loads,int configuration_type){ + + int i; + int found=0; + int mpi_found=0; + + for(i=0;iSize();i++){ + Object* object=(Object*)loads->GetObjectByOffset(i); + Load* load=(Load*)object; + if(load->InAnalysis(configuration_type)){ + if (object->Enum()==PengridEnum){ + found=1; + break; + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; +} Index: /issm/trunk/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp =================================================================== --- /issm/trunk/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp (revision 8404) +++ /issm/trunk/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp (revision 8404) @@ -0,0 +1,377 @@ +/*!\file RiftConstraintsState.cpp + * \brief: manage penalties for rifts + */ + +#include "./ConstraintsStateLocal.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../include/include.h" +#include "../../shared/shared.h" + +#define _ZIGZAGCOUNTER_ + +/*current module: */ +/*RiftIsPresent(Loads* loads,int configuration_type){{{1*/ +int RiftIsPresent(Loads* loads,int configuration_type){ + + + int i; + + int found=0; + int mpi_found=0; + + /*go though loads, and figure out if one of the loads is a Riftfront: */ + for (i=0;iSize();i++){ + Load* load=(Load*)loads->GetObjectByOffset(i); + if(load->InAnalysis(configuration_type)){ + if(RiftfrontEnum==loads->GetEnum(i)){ + found=1; + break; + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; +} +/*}}}*/ +/*RiftConstraintsState(int* pconverged, int* pnum_unstable_constraints,Loads* loads,int min_mechanical_constraints,int configuration_type){{{1*/ +void RiftConstraintsState(int* pconverged, int* pnum_unstable_constraints,Loads* loads,int min_mechanical_constraints,int configuration_type){ + + int num_unstable_constraints=0; + int converged=0; + int potential; + + RiftConstrain(&num_unstable_constraints,loads,configuration_type); + if(num_unstable_constraints==0)converged=1; + + if(RiftIsFrozen(loads,configuration_type)){ + converged=1; + num_unstable_constraints=0; + } + else if(num_unstable_constraints<=min_mechanical_constraints){ + _printf_(VerboseModule()," freezing constraints\n"); + RiftFreezeConstraints(loads,configuration_type); + } + + /*Assign output pointers: */ + *pconverged=converged; + *pnum_unstable_constraints=num_unstable_constraints; +} +/*}}}*/ +/*RiftConstrain(int* pnum_unstable_constraints,Loads* loads,int configuration_type){{{1*/ +void RiftConstrain(int* pnum_unstable_constraints,Loads* loads,int configuration_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + Load* load=NULL; + + int unstable; + int sum_num_unstable_constraints; + int num_unstable_constraints=0; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + load=(Load*)loads->GetObjectByOffset(i); + if(load->InAnalysis(configuration_type)){ + + riftfront=(Riftfront*)load; + + riftfront->Constrain(&unstable); + + num_unstable_constraints+=unstable; + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + /*Assign output pointers: */ + *pnum_unstable_constraints=num_unstable_constraints; + +} +/*}}}*/ +/*RiftIsFrozen(Loads* loads,int configuration_type){{{1*/ +int RiftIsFrozen(Loads* loads,int configuration_type){ + + int i; + + /* generic object pointer: */ + Load* load=NULL; + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + load=(Load*)loads->GetObjectByOffset(i); + if(load->InAnalysis(configuration_type)){ + + riftfront=(Riftfront*)load; + if (riftfront->IsFrozen()){ + found=1; + break; + } + } + } + } + + /*Is there just one found? that would mean we have frozen! : */ + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_MAX,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; +} +/*}}}*/ +/*RiftFreezeConstraints(Loads* loads,int configuration_type){{{1*/ +void RiftFreezeConstraints(Loads* loads,int configuration_type){ + + int i; + + /* generic object pointer: */ + Load* load=NULL; + Riftfront* riftfront=NULL; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + load=(Load*)loads->GetObjectByOffset(i); + if(load->InAnalysis(configuration_type)){ + + riftfront=(Riftfront*)load; + riftfront->FreezeConstraints(); + } + + } + } + +} +/*}}}*/ + +/*diverse trials and errors: */ +/*RiftIsMaterialStable(Loads* loads){{{1*/ +int RiftIsMaterialStable(Loads* loads){ + + int i; + + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*go though loads, and if non-linearity of the material has converged, let all penalties know: */ + for (i=0;iSize();i++){ + + if(RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + if (riftfront->IsMaterialStable()){ + found=1; + /*do not break! all penalties should get informed the non-linearity converged!*/ + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; +} +/*}}}*/ +/*RiftIsPreStable(Loads* loads){{{1*/ +int RiftIsPreStable(Loads* loads){ + + + int i; + + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*go though loads, and figure out if one of the penpair loads is still not stable: */ + for (i=0;iSize();i++){ + + if(RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + if (riftfront->PreStable()==0){ + found=1; + break; + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + if (found){ + /*We found an unstable constraint. : */ + return 0; + } + else{ + return 1; + } +} +/*}}}*/ +/*RiftSetPreStable(Loads* loads){{{1*/ +void RiftSetPreStable(Loads* loads){ + + + int i; + + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*go though loads, and set loads to pre stable.:*/ + for (i=0;iSize();i++){ + + if(RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + riftfront->SetPreStable(); + } + } +} +/*}}}*/ +/*RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads){{{1*/ +void RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + int unstable; + int sum_num_unstable_constraints; + int num_unstable_constraints=0; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->PreConstrain(&unstable); + + num_unstable_constraints+=unstable; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + /*Assign output pointers: */ + *pnum_unstable_constraints=num_unstable_constraints; + +} +/*}}}*/ +/*RiftMaxPenetrationInInputs(Loads* loads){{{1*/ +void RiftMaxPenetrationInInputs(Loads* loads){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + /*rift penetration: */ + double max_penetration=0; + double mpi_max_penetration; + double penetration; + + /*Ok, we are going to find the node pairs which are not penetrating, even though they + * are penalised. We will release only the one with has least <0 penetration. : */ + + max_penetration=0; + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->MaxPenetration(&penetration); + + if (penetration>max_penetration)max_penetration=penetration; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&max_penetration,&mpi_max_penetration,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_max_penetration,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + max_penetration=mpi_max_penetration; + #endif + + /*feed max_penetration to inputs: */ + for(i=0;iSize();i++){ + Load* load=(Load*)loads->GetObjectByOffset(i); + load->InputUpdateFromVector(&max_penetration,MaxPenetrationEnum,ConstantEnum); + } +} +/*}}}*/ +/*RiftPotentialUnstableConstraints(Loads* loads){{{1*/ +int RiftPotentialUnstableConstraints(Loads* loads){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + /*Ok, we are going to find the node pairs which are not penetrating, even though they + * are penalised. We will release only the one with has least <0 penetration. : */ + int unstable=0; + int sum_num_unstable_constraints=0; + int num_unstable_constraints=0; + + for (i=0;iSize();i++){ + + if (RiftfrontEnum==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->PotentialUnstableConstraint(&unstable); + + num_unstable_constraints+=unstable; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + return num_unstable_constraints; +} +/*}}}*/ Index: /issm/trunk/src/c/modules/ConstraintsStatex/ThermalConstraintsState.cpp =================================================================== --- /issm/trunk/src/c/modules/ConstraintsStatex/ThermalConstraintsState.cpp (revision 8404) +++ /issm/trunk/src/c/modules/ConstraintsStatex/ThermalConstraintsState.cpp (revision 8404) @@ -0,0 +1,51 @@ +/*!\file: ThermalConstraintsState.cpp + * \brief melting rate constraints + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./ConstraintsStateLocal.h" + +void ThermalConstraintsState(Loads* loads, int* pconverged, int* pnum_unstable_constraints,int configuration_type){ + + int i; + + int unstable=0; + int num_unstable_constraints=0; + int converged=0; + int sum_num_unstable_constraints=0; + + num_unstable_constraints=0; + + /*Enforce constraints: */ + for(i=0;iSize();i++){ + Object* object=(Object*)loads->GetObjectByOffset(i); + Load* load=(Load*)object; + if(load->InAnalysis(configuration_type)){ + if(object->Enum()==PengridEnum){ + + Pengrid* pengrid=(Pengrid*)object; + pengrid->PenaltyConstrain(&unstable); + num_unstable_constraints+=unstable; + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + /*Have we converged? : */ + if (num_unstable_constraints==0) converged=1; + else converged=0; + + /*Assign output pointers: */ + *pconverged=converged; + *pnum_unstable_constraints=num_unstable_constraints; +} Index: /issm/trunk/src/c/modules/modules.h =================================================================== --- /issm/trunk/src/c/modules/modules.h (revision 8403) +++ /issm/trunk/src/c/modules/modules.h (revision 8404) @@ -77,5 +77,5 @@ #include "./OutputResultsx/OutputResultsx.h" #include "./OutputRiftsx/OutputRiftsx.h" -#include "./PenaltyConstraintsx/PenaltyConstraintsx.h" +#include "./ConstraintsStatex/ConstraintsStatex.h" #include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h" #include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h" Index: /issm/trunk/src/c/solvers/solver_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_nonlinear.cpp (revision 8403) +++ /issm/trunk/src/c/solvers/solver_nonlinear.cpp (revision 8404) @@ -69,5 +69,5 @@ InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); + ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); _printf_(VerboseConvergence()," number of unstable constraints: %i\n",num_unstable_constraints); Index: /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp (revision 8403) +++ /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp (revision 8404) @@ -70,5 +70,5 @@ InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,tg); - PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); if (!converged){ Index: /issm/trunk/src/m/classes/clusters/pfe.m =================================================================== --- /issm/trunk/src/m/classes/clusters/pfe.m (revision 8403) +++ /issm/trunk/src/m/classes/clusters/pfe.m (revision 8404) @@ -119,5 +119,5 @@ fprintf(fid,'#PBS -S /bin/bash\n'); - fprintf(fid,'#PBS -N %s\n',modelname); +% fprintf(fid,'#PBS -N %s\n',modelname); fprintf(fid,'#PBS -l select=%i:ncpus=%i:model=%s\n',cluster.numnodes,cluster.cpuspernode,cluster.processor); fprintf(fid,'#PBS -l walltime=%i\n',cluster.time*60); %walltime is in seconds. Index: /issm/trunk/src/m/solvers/solver_nonlinear.m =================================================================== --- /issm/trunk/src/m/solvers/solver_nonlinear.m (revision 8403) +++ /issm/trunk/src/m/solvers/solver_nonlinear.m (revision 8404) @@ -42,5 +42,5 @@ [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug); - [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( femmodel.elements,femmodel.nodes,femmodel.vertices,loads, femmodel.materials,femmodel.parameters); + [loads,constraints_converged,num_unstable_constraints] =ConstraintsState( femmodel.elements,femmodel.nodes,femmodel.vertices,loads, femmodel.materials,femmodel.parameters); issmprintf(VerboseConvergence(),'%s%i',' number of unstable constraints: ',num_unstable_constraints); Index: /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m =================================================================== --- /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m (revision 8403) +++ /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m (revision 8404) @@ -39,5 +39,5 @@ [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,t_g); - [femmodel.loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters); + [femmodel.loads,constraints_converged,num_unstable_constraints] =ConstraintsState(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters); if ~converged, Index: /issm/trunk/src/mex/ConstraintsState/ConstraintsState.cpp =================================================================== --- /issm/trunk/src/mex/ConstraintsState/ConstraintsState.cpp (revision 8404) +++ /issm/trunk/src/mex/ConstraintsState/ConstraintsState.cpp (revision 8404) @@ -0,0 +1,65 @@ +/*\file ConstraintsState.c + *\brief: set up penalty constraints + */ + +#include "./ConstraintsState.h" + +void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ + + /*input datasets: */ + Elements* elements=NULL; + Nodes* nodes=NULL; + Vertices* vertices=NULL; + Loads* loads=NULL; + Materials* materials=NULL; + Parameters* parameters=NULL; + + /*output: */ + int converged; + int num_unstable_constraints; + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&ConstraintsStateUsage); + + /*Input datasets: */ + FetchData((DataSet**)&elements,ELEMENTS); + FetchData((DataSet**)&nodes,NODES); + FetchData((DataSet**)&vertices,VERTICES); + FetchData((DataSet**)&loads,LOADSIN); + FetchData((DataSet**)&materials,MATERIALS); + FetchParams(¶meters,PARAMETERS); + + /*configure: */ + elements-> Configure(elements,loads, nodes,vertices, materials,parameters); + nodes-> Configure(elements,loads, nodes,vertices, materials,parameters); + loads-> Configure(elements, loads, nodes,vertices, materials,parameters); + + /*!Generate internal degree of freedom numbers: */ + ConstraintsStatex(&converged, &num_unstable_constraints, elements,nodes,vertices, loads,materials,parameters); + + /*write output datasets: */ + WriteData(LOADS,loads); + WriteData(CONVERGED,converged); + WriteData(NUMUNSTABLECONSTRAINTS,num_unstable_constraints); + + /*Free ressources: */ + delete elements; + delete nodes; + delete vertices; + delete loads; + delete materials; + delete parameters; + + /*end module: */ + MODULEEND(); +} + +void ConstraintsStateUsage(void) +{ + _printf_(true,"\n"); + _printf_(true," usage: [loads, constraints_converged, num_unstable_constraints] = %s(elements,nodes,vertices,loads,materials,params);\n",__FUNCT__); + _printf_(true,"\n"); +} Index: /issm/trunk/src/mex/ConstraintsState/ConstraintsState.h =================================================================== --- /issm/trunk/src/mex/ConstraintsState/ConstraintsState.h (revision 8404) +++ /issm/trunk/src/mex/ConstraintsState/ConstraintsState.h (revision 8404) @@ -0,0 +1,37 @@ +/* + ConstraintsState.h +*/ + +#ifndef _CONSTRAINTSSTATE_H +#define _CONSTRAINTSSTATE_H + +/* local prototypes: */ +void ConstraintsStateUsage(void); + +#include "../../c/modules/modules.h" +#include "../../c/Container/Container.h" +#include "../../c/shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "ConstraintsState" + +/* serial input macros: */ +#define ELEMENTS (mxArray*)prhs[0] +#define NODES (mxArray*)prhs[1] +#define VERTICES (mxArray*)prhs[2] +#define LOADSIN (mxArray*)prhs[3] +#define MATERIALS (mxArray*)prhs[4] +#define PARAMETERS (mxArray*)prhs[5] + +/* serial output macros: */ +#define LOADS (mxArray**)&plhs[0] +#define CONVERGED (mxArray**)&plhs[1] +#define NUMUNSTABLECONSTRAINTS (mxArray**)&plhs[2] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 3 +#undef NRHS +#define NRHS 6 + +#endif /* _CONSTRAINTSSTATE_H */ Index: /issm/trunk/src/mex/Makefile.am =================================================================== --- /issm/trunk/src/mex/Makefile.am (revision 8403) +++ /issm/trunk/src/mex/Makefile.am (revision 8404) @@ -53,5 +53,5 @@ OutputResults\ ParsePetscOptions\ - PenaltyConstraints\ + ConstraintsState\ PointCloudFindNeighbors\ PropagateFlagsFromConnectivity\ @@ -251,6 +251,6 @@ OutputResults/OutputResults.h -PenaltyConstraints_SOURCES = PenaltyConstraints/PenaltyConstraints.cpp\ - PenaltyConstraints/PenaltyConstraints.h +ConstraintsState_SOURCES = ConstraintsState/ConstraintsState.cpp\ + ConstraintsState/ConstraintsState.h PointCloudFindNeighbors_SOURCES = PointCloudFindNeighbors/PointCloudFindNeighbors.cpp\