Changeset 9002


Ignore:
Timestamp:
07/15/11 08:49:34 (14 years ago)
Author:
Eric.Larour
Message:

Implemented transient boundary conditions. Applied to thermal core only for now

Location:
issm/trunk/src
Files:
13 added
32 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Container/DataSet.cpp

    r8798 r9002  
    288288                                dataset->AddObject(spc);}
    289289                                break;
     290                        case SpctEnum:{
     291                                Spct* spct=NULL;
     292                                spct=new Spct();
     293                                spct->Demarshall(&marshalled_dataset);
     294                                dataset->AddObject(spct);}
     295                                break;
     296
    290297                        case PengridEnum:{
    291298                                Pengrid* pengrid=NULL;

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r8942 r9002  
    154154        /*Spc: */
    155155        SpcEnum,
     156        SpctEnum,
    156157        /*}}}*/
    157158        /*Geography {{{1*/

  • issm/trunk/src/c/Makefile.am

    r8967 r9002  
    230230                                        ./objects/Materials/Matpar.h\
    231231                                        ./objects/Materials/Matpar.cpp\
     232                                        ./objects/Constraints/Constraint.h\
    232233                                        ./objects/Constraints/Spc.cpp\
    233234                                        ./objects/Constraints/Spc.h\
     235                                        ./objects/Constraints/Spct.cpp\
     236                                        ./objects/Constraints/Spct.h\
    234237                                        ./objects/Loads/Penpair.cpp\
    235238                                        ./objects/Loads/Penpair.h\
    … …  
    566569                                        ./modules/IoModelToInputsx/IoModelToInputsx.h\
    567570                                        ./modules/IoModelToInputsx/IoModelToInputsx.cpp\
     571                                        ./modules/IoModelToConstraintsx/IoModelToConstraintsx.h\
     572                                        ./modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp\
    568573                                        ./modules/InputUpdateFromConstantx/InputUpdateFromConstantx.h\
    569574                                        ./modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp\
    … …  
    596601                                        ./modules/SpcNodesx/SpcNodesx.h\
    597602                                        ./modules/SpcNodesx/SpcNodesx.cpp\
     603                                        ./modules/UpdateConstraintsx/UpdateConstraintsx.h\
     604                                        ./modules/UpdateConstraintsx/UpdateConstraintsx.cpp\
    598605                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
    599606                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\
    … …  
    909916                                        ./objects/Materials/Matpar.h\
    910917                                        ./objects/Materials/Matpar.cpp\
     918                                        ./objects/Constraints/Constraint.h\
    911919                                        ./objects/Constraints/Spc.cpp\
    912920                                        ./objects/Constraints/Spc.h\
     921                                        ./objects/Constraints/Spct.cpp\
     922                                        ./objects/Constraints/Spct.h\
    913923                                        ./objects/Loads/Penpair.cpp\
    914924                                        ./objects/Loads/Penpair.h\
    … …  
    12381248                                        ./modules/IoModelToInputsx/IoModelToInputsx.h\
    12391249                                        ./modules/IoModelToInputsx/IoModelToInputsx.cpp\
     1250                                        ./modules/IoModelToConstraintsx/IoModelToConstraintsx.h\
     1251                                        ./modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp\
    12401252                                        ./modules/InputUpdateFromConstantx/InputUpdateFromConstantx.h\
    12411253                                        ./modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp\
    … …  
    12681280                                        ./modules/SpcNodesx/SpcNodesx.h\
    12691281                                        ./modules/SpcNodesx/SpcNodesx.cpp\
     1282                                        ./modules/UpdateConstraintsx/UpdateConstraintsx.h\
     1283                                        ./modules/UpdateConstraintsx/UpdateConstraintsx.cpp\
    12701284                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
    12711285                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\

  • issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r8942 r9002  
    130130                case ExternalResultEnum : return "ExternalResult";
    131131                case SpcEnum : return "Spc";
     132                case SpctEnum : return "Spct";
    132133                case GeographyEnum : return "Geography";
    133134                case IceSheetEnum : return "IceSheet";

  • issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r8926 r9002  
    7373        parameters->AddObject(new BoolParam(IsprognosticEnum,iomodel->isprognostic));
    7474        parameters->AddObject(new BoolParam(IsthermalEnum,iomodel->isthermal));
     75        parameters->AddObject(new DoubleParam(TimeEnum,0.0));  //start at time 0 by default for all solutions.
    7576
    7677        /*Deal with more complex parameters*/

  • issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp

    r8926 r9002  
    44
    55#include "../../../Container/Container.h"
     6#include "../../../modules/modules.h"
    67#include "../../../io/io.h"
    78#include "../../../toolkits/toolkits.h"
    … …  
    1920        /*Output*/
    2021        Constraints* constraints = NULL;
    21         Spc*    spc  = NULL;
    2222
    2323        /*Recover pointer: */
    … …  
    3030        if (iomodel->dim==2) goto cleanup_and_return;
    3131
    32         /*Fetch data: */
    33         IoModelFetchData(&iomodel->spctemperature,NULL,NULL,iomodel_handle,SpctemperatureEnum);
    34 
    35         /*Initialize counter*/
    36         count=0;
    37 
    38         /*Create spcs from x,y,z, as well as the spc values on those spcs: */
    39         for (i=0;i<iomodel->numberofvertices;i++){
    40                 /*keep only this partition's nodes:*/
    41                 if((iomodel->my_vertices[i])){
    42 
    43                         if ((int)iomodel->spctemperature[2*i]){
    44 
    45                                 constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->spctemperature[2*i+1],ThermalAnalysisEnum));
    46                                 count++;
    47 
    48                         }
    49 
    50                 } //if((my_nodes[i]==1))
    51         }
    52 
    53         /*Free data: */
    54         xfree((void**)&iomodel->spctemperature);
     32        /*Create constraints: */
     33        IoModelToConstraintsx(constraints,iomodel,iomodel_handle,SpctemperatureEnum,ThermalAnalysisEnum);
    5534
    5635        cleanup_and_return:
    57        
     36
    5837        /*Assign output pointer: */
    5938        *pconstraints=constraints;

  • issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp

    r8926 r9002  
    3939                if((iomodel->my_vertices[i]==1)){
    4040
    41                         if (!iomodel->spctemperature[2*i]){ //No penalty applied on spc nodes!
     41                        if (isnan(iomodel->spctemperature[i])){ //No penalty applied on spc nodes!
    4242
    4343                                loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,ThermalAnalysisEnum));

  • issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.cpp

    r8303 r9002  
    1010#include "../../EnumDefinitions/EnumDefinitions.h"
    1111
    12 void SpcNodesx(Nodes* nodes,Constraints* constraints,int analysis_type){
     12void SpcNodesx(Nodes* nodes,Constraints* constraints,Parameters* parameters, int analysis_type){
    1313
    1414        int i;
    15         Node* node=NULL;
    16         int nodeid;
    17         int dof;
    18         double value;
     15                       
     16        for(i=0;i<constraints->Size();i++){
     17       
     18                Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i);
    1919
    20         for(i=0;i<constraints->Size();i++){
     20                /*Check this constraint belongs to this analysis: */
     21                if(constraint->InAnalysis(analysis_type)){
    2122
    22                 Object* object=(Object*)constraints->GetObjectByOffset(i);
    23 
    24                 /*Check constraints is a single point constraint (spc): */
    25                 if(object->Enum()==SpcEnum){
    26 
    27                         Spc* spc=(Spc*)object;
    28 
    29                         if(spc->InAnalysis(analysis_type)){
    30 
    31                                 /*Ok, constraints object is a constraint. Get the nodeid from the node it applies to: */
    32                                 nodeid=spc->GetNodeId();
    33                                 dof=spc->GetDof();
    34                                 value=spc->GetValue();
    35 
    36                                 /*Now, chase through nodes and find the corect node: */
    37                                 node=(Node*)nodes->GetObjectById(NULL,nodeid);
    38 
    39                                 /*Apply constraint: */
    40                                 if(node){ //in case the spc is dealing with a node on another cpu
    41                                         node->ApplyConstraint(dof,value);
    42                                 }
    43                         }
    44 
     23                        /*Ok, apply constraint onto corresponding node: */
     24                        constraint->ConstrainNode(nodes,parameters);
    4525                }
    4626        }

  • issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.h

    r5772 r9002  
    1111
    1212/* local prototypes: */
    13 void SpcNodesx(Nodes* nodesin,Constraints* constraints,int analysis_type);
     13void SpcNodesx(Nodes* nodes, Constraints* constraints,Parameters* parameters,int analysis_type);
    1414
    1515#endif  /* _SPCNODESX_H */

  • issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r8942 r9002  
    128128        else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
    129129        else if (strcmp(name,"Spc")==0) return SpcEnum;
     130        else if (strcmp(name,"Spct")==0) return SpctEnum;
    130131        else if (strcmp(name,"Geography")==0) return GeographyEnum;
    131132        else if (strcmp(name,"IceSheet")==0) return IceSheetEnum;

  • issm/trunk/src/c/modules/modules.h

    r8967 r9002  
    5050#include "./InputArtificialNoisex/InputArtificialNoisex.h"
    5151#include "./IoModelToInputsx/IoModelToInputsx.h"
     52#include "./IoModelToConstraintsx/IoModelToConstraintsx.h"
    5253#include "./KMLMeshWritex/KMLFileReadx.h"
    5354#include "./KMLMeshWritex/KMLMeshWritex.h"
    … …  
    9596#include "./Solverx/Solverx.h"
    9697#include "./SpcNodesx/SpcNodesx.h"
     98#include "./UpdateConstraintsx/UpdateConstraintsx.h"
    9799#include "./SurfaceAreax/SurfaceAreax.h"
    98100#include "./SystemMatricesx/SystemMatricesx.h"

  • issm/trunk/src/c/objects/Constraints/Spc.cpp

    r8224 r9002  
    147147/*}}}1*/
    148148
     149/*Constraint virtual functions definitions: */
     150/*FUNCTION Spc::InAnalysis{{{1*/
     151bool Spc::InAnalysis(int in_analysis_type){
     152        if (in_analysis_type==this->analysis_type) return true;
     153        else return false;
     154}
     155/*}}}*/
     156/*FUNCTION Spc::ConstrainNode{{{1*/
     157void Spc::ConstrainNode(Nodes* nodes,Parameters* parameters){
     158
     159        Node* node=NULL;
     160
     161        /*Chase through nodes and find the node to which this Spc applys: */
     162        node=(Node*)nodes->GetObjectById(NULL,nodeid);
     163
     164        /*Apply constraint: */
     165        if(node){ //in case the spc is dealing with a node on another cpu
     166                node->ApplyConstraint(dof,value);
     167        }
     168}
     169/*}}}*/
     170
    149171/*Spc functions*/
    150172/*FUNCTION Spc::GetDof {{{1*/
    … …  
    164186}
    165187/*}}}1*/
    166 /*FUNCTION Spc::InAnalysis{{{1*/
    167 bool Spc::InAnalysis(int in_analysis_type){
    168         if (in_analysis_type==this->analysis_type) return true;
    169         else return false;
    170 }
    171 /*}}}*/

  • issm/trunk/src/c/objects/Constraints/Spc.h

    r4248 r9002  
    1212/*}}}*/
    1313
    14 class Spc: public Object{
     14class Spc: public Constraint{
    1515
    1616        private:
    … …  
    3939                Object* copy();
    4040                /*}}}*/
     41                /*Constraint virtual functions definitions: {{{1*/
     42                void   ConstrainNode(Nodes* nodes,Parameters* parameters);
     43                bool   InAnalysis(int analysis_type);
     44                /*}}}*/
    4145                /*Spc management:{{{1 */
    4246                int    GetNodeId();
    4347                int    GetDof();
    4448                double GetValue();
    45                 bool    InAnalysis(int analysis_type);
    4649                /*}}}*/
    4750

  • issm/trunk/src/c/objects/FemModel.cpp

    r8926 r9002  
    5959
    6060                _printf_(VerboseMProcessor(),"      resolve node constraints\n");
    61                 SpcNodesx(nodes,constraints,analysis_type);
    62        
     61                SpcNodesx(nodes,constraints,parameters,analysis_type);
     62
    6363                _printf_(VerboseMProcessor(),"      create nodal degrees of freedom\n");
    6464                NodesDofx(nodes,parameters,analysis_type);

  • issm/trunk/src/c/objects/Loads/Load.h

    r8800 r9002  
    2323
    2424                virtual       ~Load(){};
     25               
    2526                /*Virtual functions: {{{1*/
    2627                virtual void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;

  • issm/trunk/src/c/objects/Node.cpp

    r8827 r9002  
    529529}
    530530/*}}}*/
     531/*FUNCTION Node::RelaxConstraint{{{1*/
     532void  Node::RelaxConstraint(int dof){
     533
     534        /*Dof should be added to the f-set, and taken out of the s-set:*/
     535        DofInFSet(dof-1);
     536        this->indexing.svalues[dof-1]=NAN;
     537}
     538/*}}}*/
    531539/*FUNCTION Node::CreateVecSets {{{1*/
    532540void  Node::CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s){
    … …  
    593601}
    594602/*}}}*/
     603/*FUNCTION Node::DofInFSet {{{1*/
     604void  Node::DofInFSet(int dof){
     605
     606        /*Put dof for this node into the f set (ie, this dof will NOT be constrained
     607         * to a fixed value during computations. */
     608
     609        this->indexing.f_set[dof]=1;
     610        this->indexing.s_set[dof]=0;
     611}
     612/*}}}*/
    595613/*FUNCTION Node::FreezeDof{{{1*/
    596614void  Node::FreezeDof(int dof){

  • issm/trunk/src/c/objects/Node.h

    r8800 r9002  
    7575                int   IsClone();
    7676                void  ApplyConstraint(int dof,double value);
     77                void  RelaxConstraint(int dof);
    7778                void  DofInSSet(int dof);
     79                void  DofInFSet(int dof);
    7880                int   GetDof(int dofindex,int setenum);
    7981                void  CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s);

  • issm/trunk/src/c/objects/objects.h

    r8800 r9002  
    2020
    2121/*Constraints: */
     22#include "./Constraints/Constraint.h"
    2223#include "./Constraints/Spc.h"
     24#include "./Constraints/Spct.h"
    2325
    2426/*Gauss*/

  • issm/trunk/src/c/shared/Numerics/UnitConversion.cpp

    r8967 r9002  
    4747        double scale;
    4848        switch(type_enum){
     49                case TimeEnum:        scale=1.0/yts;break; //yr
    4950                case VxEnum:          scale=yts;break; //m/yr
    5051                case VxObsEnum:       scale=yts;break; //m/yr

  • issm/trunk/src/c/solutions/thermal_core.cpp

    r8392 r9002  
    4343                if(nsteps)_printf_(VerboseSolution(),"time step: %i/%i\n",i+1,nsteps);
    4444                time=(i+1)*dt;
     45                femmodel->parameters->SetParam(time,TimeEnum);
    4546
    4647                /*call thermal_core_step: */

  • issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp

    r8815 r9002  
    4848        InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,reset_penalties,ResetPenaltiesEnum);
    4949        InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,false,ConvergedEnum);
     50
     51        UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
     52
    5053        for(;;){
    5154

  • issm/trunk/src/m/model/extrude.m

    r8823 r9002  
    136136
    137137%drag_coefficient is limited to nodes that are on the bedrock.
    138 md.drag_coefficient=project3d(md,md.drag_coefficient,'node',1);
     138md.drag_coefficient=project3d(md,'vector',md.drag_coefficient,'type','node','layer',1);
    139139
    140140%p and q (same deal, except for element that are on the bedrock: )
    141 md.drag_p=project3d(md,md.drag_p,'element');
    142 md.drag_q=project3d(md,md.drag_q,'element');
     141md.drag_p=project3d(md,'vector',md.drag_p,'type','element');
     142md.drag_q=project3d(md,'vector',md.drag_q,'type','element');
    143143
    144144%observations
    145 md.vx_obs=project3d(md,md.vx_obs,'node');
    146 md.vy_obs=project3d(md,md.vy_obs,'node');
    147 md.vel_obs=project3d(md,md.vel_obs,'node');
    148 md.vx_bal=project3d(md,md.vx_bal,'node');
    149 md.vy_bal=project3d(md,md.vy_bal,'node');
    150 md.vel_bal=project3d(md,md.vel_bal,'node');
    151 md.vel_obs_raw=project3d(md,md.vel_obs_raw,'node');
    152 md.surface_mass_balance=project3d(md,md.surface_mass_balance,'node');
    153 md.surface_accumulation_rate=project3d(md,md.surface_accumulation_rate,'node');
    154 md.surface_ablation_rate=project3d(md,md.surface_ablation_rate,'node');
    155 md.dhdt=project3d(md,md.dhdt,'node');
    156 md.firn_layer=project3d(md,md.firn_layer,'node',md.numlayers);
    157 if ~isempty(md.forcings),
    158         forcingnames=fieldnames(md.forcings);
    159         numforcings=length(fieldnames(md.forcings));
    160         for i=1:numforcings
    161                 forcing=md.forcings.(forcingnames{i});
    162                 md.forcings.(forcingnames{i})=[project3d(md,forcing(1:end-1,:),'node');forcing(end,:)];
    163         end
    164 end
     145md.vx_obs=project3d(md,'vector',md.vx_obs,'type','node');
     146md.vy_obs=project3d(md,'vector',md.vy_obs,'type','node');
     147md.vel_obs=project3d(md,'vector',md.vel_obs,'type','node');
     148md.vx_bal=project3d(md,'vector',md.vx_bal,'type','node');
     149md.vy_bal=project3d(md,'vector',md.vy_bal,'type','node');
     150md.vel_bal=project3d(md,'vector',md.vel_bal,'type','node');
     151md.vel_obs_raw=project3d(md,'vector',md.vel_obs_raw,'type','node');
     152md.surface_mass_balance=project3d(md,'vector',md.surface_mass_balance,'type','node');
     153md.surface_accumulation_rate=project3d(md,'vector',md.surface_accumulation_rate,'type','node');
     154md.surface_ablation_rate=project3d(md,'vector',md.surface_ablation_rate,'type','node');
     155md.dhdt=project3d(md,'vector',md.dhdt,'type','node');
     156md.firn_layer=project3d(md,'vector',md.firn_layer,'type','node','layer',md.numlayers);
    165157
    166158%results
    167 if ~isnan(md.vx),md.vx=project3d(md,md.vx,'node');end;
    168 if ~isnan(md.vy),md.vy=project3d(md,md.vy,'node');end;
    169 if ~isnan(md.vz),md.vz=project3d(md,md.vz,'node');end;
    170 if ~isnan(md.vel),md.vel=project3d(md,md.vel,'node');end;
    171 if ~isnan(md.temperature),md.temperature=project3d(md,md.temperature,'node');end;
    172 if ~isnan(md.waterfraction),md.waterfraction=project3d(md,md.waterfraction,'node');end;
    173 if ~isnan(md.surface_slopex),md.surface_slopex=project3d(md,md.surface_slopex,'node');end;
    174 if ~isnan(md.surface_slopey),md.surface_slopey=project3d(md,md.surface_slopey,'node');end;
    175 if ~isnan(md.bed_slopex),md.bed_slopex=project3d(md,md.bed_slopex,'node');end;
    176 if ~isnan(md.bed_slopey),md.bed_slopey=project3d(md,md.bed_slopey,'node');end;
     159if ~isnan(md.vx),md.vx=project3d(md,'vector',md.vx,'type','node');end;
     160if ~isnan(md.vy),md.vy=project3d(md,'vector',md.vy,'type','node');end;
     161if ~isnan(md.vz),md.vz=project3d(md,'vector',md.vz,'type','node');end;
     162if ~isnan(md.vel),md.vel=project3d(md,'vector',md.vel,'type','node');end;
     163if ~isnan(md.temperature),md.temperature=project3d(md,'vector',md.temperature,'type','node');end;
     164if ~isnan(md.waterfraction),md.waterfraction=project3d(md,'vector',md.waterfraction,'type','node');end;
     165if ~isnan(md.surface_slopex),md.surface_slopex=project3d(md,'vector',md.surface_slopex,'type','node');end;
     166if ~isnan(md.surface_slopey),md.surface_slopey=project3d(md,'vector',md.surface_slopey,'type','node');end;
     167if ~isnan(md.bed_slopex),md.bed_slopex=project3d(md,'vector',md.bed_slopex,'type','node');end;
     168if ~isnan(md.bed_slopey),md.bed_slopey=project3d(md,'vector',md.bed_slopey,'type','node');end;
    177169
    178170%bedinfo and surface info
    179 md.elementonbed=project3d(md,ones(md.numberofelements2d,1),'element',1);
    180 md.elementonsurface=project3d(md,ones(md.numberofelements2d,1),'element',md.numlayers-1);
    181 md.nodeonbed=project3d(md,ones(md.numberofnodes2d,1),'node',1);
    182 md.nodeonsurface=project3d(md,ones(md.numberofnodes2d,1),'node',md.numlayers);
     171md.elementonbed=project3d(md,'vector',ones(md.numberofelements2d,1),'type','element','layer',1);
     172md.elementonsurface=project3d(md,'vector',ones(md.numberofelements2d,1),'type','element','layer',md.numlayers-1);
     173md.nodeonbed=project3d(md,'vector',ones(md.numberofnodes2d,1),'type','node','layer',1);
     174md.nodeonsurface=project3d(md,'vector',ones(md.numberofnodes2d,1),'type','node','layer',md.numlayers);
    183175
    184176%elementstype
    … …  
    186178        oldelements_type=md.elements_type2d;
    187179        md.elements_type=zeros(number_el3d,1);
    188         md.elements_type=project3d(md,oldelements_type,'element');
    189         md.nodeonhutter=project3d(md,md.nodeonhutter,'node');
    190         md.nodeonmacayeal=project3d(md,md.nodeonmacayeal,'node');
    191         md.nodeonpattyn=project3d(md,md.nodeonpattyn,'node');
    192         md.nodeonstokes=project3d(md,md.nodeonstokes,'node');
     180        md.elements_type=project3d(md,'vector',oldelements_type,'type','element');
     181        md.nodeonhutter=project3d(md,'vector',md.nodeonhutter,'type','node');
     182        md.nodeonmacayeal=project3d(md,'vector',md.nodeonmacayeal,'type','node');
     183        md.nodeonpattyn=project3d(md,'vector',md.nodeonpattyn,'type','node');
     184        md.nodeonstokes=project3d(md,'vector',md.nodeonstokes,'type','node');
    193185end
    194186
    … …  
    197189        oldvertices_type=md.vertices_type2d;
    198190        md.vertices_type=zeros(number_nodes3d,1);
    199         md.vertices_type=project3d(md,oldvertices_type,'node');
     191        md.vertices_type=project3d(md,'vector',oldvertices_type,'type','node');
    200192end
    201193
    202194%boundary conditions
    203 md.spcvx=project3d(md,md.spcvx,'node');
    204 md.spcvy=project3d(md,md.spcvy,'node');
    205 md.spcvz=project3d(md,md.spcvz,'node');
    206 md.spctemperature=project3d(md,md.spctemperature,'node',md.numlayers);
    207 md.spcthickness=project3d(md,md.spcthickness,'node');
    208 md.diagnostic_ref=project3d(md,md.diagnostic_ref,'node');
     195md.spcvx=project3d(md,'vector',md.spcvx,'type','node');
     196md.spcvy=project3d(md,'vector',md.spcvy,'type','node');
     197md.spcvz=project3d(md,'vector',md.spcvz,'type','node');
     198md.spctemperature=project3d(md,'vector',md.spctemperature,'type','node','layer',md.numlayers,'padding',NaN);
     199md.spcthickness=project3d(md,'vector',md.spcthickness,'type','node');
     200md.diagnostic_ref=project3d(md,'vector',md.diagnostic_ref,'type','node');
    209201
    210202%in 3d, pressureload: [node1 node2 node3 node4 element]
    … …  
    226218
    227219%materials
    228 md.rheology_B=project3d(md,md.rheology_B,'node');
    229 md.rheology_n=project3d(md,md.rheology_n,'element');
     220md.rheology_B=project3d(md,'vector',md.rheology_B,'type','node');
     221md.rheology_n=project3d(md,'vector',md.rheology_n,'type','element');
    230222
    231223%parameters
    232 md.surface=project3d(md,md.surface,'node');
    233 md.thickness=project3d(md,md.thickness,'node');
    234 md.bed=project3d(md,md.bed,'node');
    235 md.nodeonboundary=project3d(md,md.nodeonboundary,'node');
    236 md.elementoniceshelf=project3d(md,md.elementoniceshelf,'element');
    237 md.nodeoniceshelf=project3d(md,md.nodeoniceshelf,'node');
    238 md.elementonicesheet=project3d(md,md.elementonicesheet,'element');
    239 md.nodeonicesheet=project3d(md,md.nodeonicesheet,'node');
    240 md.elementonwater=project3d(md,md.elementonwater,'element');
    241 md.nodeonwater=project3d(md,md.nodeonwater,'node');
    242 if ~isnan(md.weights),md.weights=project3d(md,md.weights,'node');end;
     224md.surface=project3d(md,'vector',md.surface,'type','node');
     225md.thickness=project3d(md,'vector',md.thickness,'type','node');
     226md.bed=project3d(md,'vector',md.bed,'type','node');
     227md.nodeonboundary=project3d(md,'vector',md.nodeonboundary,'type','node');
     228md.elementoniceshelf=project3d(md,'vector',md.elementoniceshelf,'type','element');
     229md.nodeoniceshelf=project3d(md,'vector',md.nodeoniceshelf,'type','node');
     230md.elementonicesheet=project3d(md,'vector',md.elementonicesheet,'type','element');
     231md.nodeonicesheet=project3d(md,'vector',md.nodeonicesheet,'type','node');
     232md.elementonwater=project3d(md,'vector',md.elementonwater,'type','element');
     233md.nodeonwater=project3d(md,'vector',md.nodeonwater,'type','node');
     234if ~isnan(md.weights),md.weights=project3d(md,'vector',md.weights,'type','node');end;
    243235
    244236%Put lithostatic pressure is there is an existing pressure
    … …  
    248240
    249241%special for thermal modeling:
    250 md.basal_melting_rate=project3d(md,md.basal_melting_rate,'node',1);
    251 md.observed_temperature=project3d(md,md.observed_temperature,'node');
     242md.basal_melting_rate=project3d(md,'vector',md.basal_melting_rate,'type','node','layer',1);
     243md.observed_temperature=project3d(md,'vector',md.observed_temperature,'type','node');
    252244if ~isnan(md.geothermalflux)
    253         md.geothermalflux=project3d(md,md.geothermalflux,'node',1); %bedrock only gets geothermal flux
     245        md.geothermalflux=project3d(md,'vector',md.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux
    254246end
    255247
    256248%mask:
    257 md.mask=project3d(md,md.mask,'node');
     249md.mask=project3d(md,'vector',md.mask,'type','node');
    258250
    259251%increase connectivity if less than 25:

  • issm/trunk/src/m/model/ismodelselfconsistent.m

    r9001 r9002  
    485485
    486486                        %CHECK THAT WE ARE NOT FULLY CONSTRAINED
    487                         if isnan(find(~md.spctemperature(:,1))),
     487                        if ~any(~isnan(md.spctemperature))
    488488                                message(['model not consistent: model ' md.name ' is totally constrained for temperature, no need to solve!']);
    489489                        end
    … …  
    507507
    508508                                %INITIAL TEMPERATURE
    509                                 fields={'temperature','spctemperature(:,2)','observed_temperature'};
     509                                fields={'temperature','observed_temperature'};
    510510                                checkgreater(md,fields,0)
     511
     512                                %CHECK SPCTEMPERATURE that are not NaN are >0.
     513                                if find(any(md.spctemperature(find(~isnan(md.spctemperature)))<=0)),
     514                                        message(['model not consistent: model ' md.name ' is constrained with negative or nil temperatures!']);
     515                                end
     516
    511517                        end
    512518                        %}}}

  • issm/trunk/src/m/model/project3d.m

    r8298 r9002  
    1 function projected_vector=project3d(md3d,vector2d,type,varargin);
     1function projected_vector=project3d(md3d,varargin);
    22%PROJECT3D - vertically project a vector from 2d mesh
    33%
    … …  
    55%   This vector can be a node vector of size (md3d.numberofnodes2d,N/A) or an
    66%   element vector of size (md3d.numberofelements2d,N/A).
    7 %   type is 'element' or 'node'. layer an optional layer number where vector
    8 %   should keep its values. If not specified, all layers adopt the value of the
    9 %   2d vector.
     7%   arguments:
     8%      'vector': 2d vector
     9%      'type': 'element' or 'node'.
     10%   options:
     11%      'layer' a layer number where vector should keep its values. If not specified, all layers adopt the
     12%             value of the 2d vector.
     13%      'padding': default to 0 (value adopted by other 3d layers not being projected0
    1014%
    11 %   Usage:
    12 %      extruded_vector=project3d(md3d,vector2d,type,layer);
     15%   Egs:
     16%      extruded_vector=project3d(md3d,'vector',vector2d,'type','node','layer',1,'padding',NaN);
     17%      extruded_vector=project3d(md3d,'vector',vector2d,'type','element','padding',0);
     18%      extruded_vector=project3d(md3d,'vector',vector2d,'type','node');
    1319
    14 if nargin==4,
    15         layer=varargin{1};
    16         if ((layer<1) || (layer>md3d.numlayers)),
    17                 error(['project3d error message: layer shoud be between 1 and ' num2str(md3d.numlayers)]);
    18         end
    19 else
    20         layer=0;
     20
     21if nargin==0,
     22        Project3dUsage;
     23        error();
    2124end
     25
     26options=pairoptions(varargin{:});
     27
     28%retrieve parameters from options.
     29vector2d=getfieldvalue(options,'vector'); %mandatory
     30type=getfieldvalue(options,'type'); %mandatory
     31layer=getfieldvalue(options,'layer',0); %optional (do all layers otherwise)
     32paddingvalue=getfieldvalue(options,'padding',0); %0 by default
    2233
    2334if strcmpi(type,'node'),
    2435
    25         projected_vector=zeros(md3d.numberofnodes,size(vector2d,2));
     36        projected_vector=paddingvalue*ones(md3d.numberofnodes,size(vector2d,2));
    2637       
    2738        if layer==0,
    … …  
    3243                projected_vector(((layer-1)*md3d.numberofnodes2d+1):(layer*md3d.numberofnodes2d),:)=vector2d;
    3344        end
    34 else
     45elseif strcmpi(type,'element'),
    3546
    36         projected_vector=zeros(md3d.numberofelements,size(vector2d,2));
     47        projected_vector=paddingvalue*ones(md3d.numberofelements,size(vector2d,2));
    3748
    3849        if layer==0,
    … …  
    4455                projected_vector( ((layer-1)*md3d.numberofelements2d+1):(layer*md3d.numberofelements2d),:)=vector2d;
    4556        end
     57else
     58        error('project3d error message: unknown projection type');
    4659end
     60
     61function Project3dUsage;
     62
     63disp('function projected_vector=project3d(md3d,varargin);');
     64disp('PROJECT3D - vertically project a vector from 2d mesh...');
     65disp('');
     66disp('   vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh.');
     67disp('   This vector can be a node vector of size (md3d.numberofnodes2d,N/A) or an ');
     68disp('   element vector of size (md3d.numberofelements2d,N/A). ');
     69disp('   arguments: ');
     70disp('      ''vector'': 2d vector');
     71disp('      ''type'': ''element'' or ''node''. ');
     72disp('   options: ');
     73disp('      ''layer'' a layer number where vector should keep its values. If not specified, all layers adopt the ');
     74disp('             value of the 2d vector.');
     75disp('      ''padding'': default to 0 (value adopted by other 3d layers not being projected0');
     76disp('');
     77disp('   Egs:');
     78disp('      extruded_vector=project3d(md3d,''vector'',vector2d,''type'',''node'',''layer'',1,''padding'',NaN);');
     79disp('      extruded_vector=project3d(md3d,''vector'',vector2d,''type'',''element'',''padding'',0);');
     80disp('      extruded_vector=project3d(md3d,''vector'',vector2d,''type'',''node'');');

  • issm/trunk/src/m/solutions/NewFemModel.m

    r8803 r9002  
    3737
    3838                issmprintf(VerboseMProcessor(),'%s','      resolve node constraints');
    39                 [femmodel.nodes]=SpcNodes(femmodel.nodes,femmodel.constraints,analysis_type);
     39                [femmodel.nodes]=SpcNodes(femmodel.nodes,femmodel.constraints,femmodel.parameters,analysis_type);
    4040
    4141                issmprintf(VerboseMProcessor(),'%s','      create nodal degrees of freedom');

  • issm/trunk/src/m/solvers/solver_thermal_nonlinear.m

    r8814 r9002  
    2020        [femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,0,ConvergedEnum);
    2121        [femmodel.elements femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,1,ResetPenaltiesEnum);
     22
     23        %update constraints
     24        [femmodel.nodes]=UpdateConstraints(femmodel.nodes,femmodel.constraints,femmodel.parameters);
    2225
    2326        while(~converged),

  • issm/trunk/src/m/utils/BC/SetIceSheetBC.m

    r8823 r9002  
    6161
    6262if (length(md.observed_temperature)==md.numberofnodes),
    63         md.spctemperature=[md.nodeonsurface md.observed_temperature]; %impose observed temperature on surface
     63        md.spctemperature=NaN*ones(md.numberofnodes,1);
     64        pos=find(md.nodeonsurface); md.spctemperature(pos)=md.observed_temperature(pos); %impose observed temperature on surface
    6465        if (length(md.geothermalflux)~=md.numberofnodes),
    6566                md.geothermalflux=50*10^-3*ones(md.numberofnodes,1); %50 mW/m^2

  • issm/trunk/src/m/utils/BC/SetIceShelfBC.m

    r8823 r9002  
    9393
    9494if (length(md.observed_temperature)==md.numberofnodes),
    95         md.spctemperature=[md.nodeonsurface md.observed_temperature]; %impose observed temperature on surface
     95        md.spctemperature=NaN*ones(md.numberofnodes,1);
     96        pos=find(md.nodeonsurface); md.spctemperature(pos)=md.observed_temperature(pos); %impose observed temperature on surface
    9697        if (length(md.geothermalflux)~=md.numberofnodes),
    9798                md.geothermalflux=zeros(md.numberofnodes,1);

  • issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m

    r8823 r9002  
    104104
    105105if (length(md.observed_temperature)==md.numberofnodes),
    106         md.spctemperature=[md.nodeonsurface md.observed_temperature]; %impose observed temperature on surface
     106        md.spctemperature=NaN*ones(md.numberofnodes,1);
     107        pos=find(md.nodeonsurface); md.spctemperature(pos)=md.observed_temperature(pos); %impose observed temperature on surface
    107108        if (length(md.geothermalflux)~=md.numberofnodes),
    108109                md.geothermalflux=zeros(md.numberofnodes,1);

  • issm/trunk/src/mex/Makefile.am

    r8801 r9002  
    6969                                SparseToVector\
    7070                                SpcNodes\
     71                                UpdateConstraints\
    7172                                SurfaceArea\
    7273                                SystemMatrices\
    … …  
    301302                          SpcNodes/SpcNodes.h
    302303
     304UpdateConstraints_SOURCES = UpdateConstraints/UpdateConstraints.cpp\
     305                          UpdateConstraints/UpdateConstraints.h
     306
    303307SystemMatrices_SOURCES = SystemMatrices/SystemMatrices.cpp\
    304308                          SystemMatrices/SystemMatrices.h

  • issm/trunk/src/mex/SpcNodes/SpcNodes.cpp

    r8910 r9002  
    1010        /*input datasets: */
    1111        Nodes* nodes=NULL;
     12        Parameters* parameters=NULL;
    1213        Constraints* constraints=NULL;
    1314        int      analysis_type;
    … …  
    2223        FetchMatlabData((DataSet**)&nodes,NODESIN);
    2324        FetchMatlabData((DataSet**)&constraints,CONSTRAINTS);
     25        FetchMatlabData(&parameters,PARAMETERS);
    2426        FetchMatlabData(&analysis_type,ANALYSISTYPE);
    2527       
    2628        /*!Generate internal degree of freedom numbers: */
    27         SpcNodesx(nodes,constraints,analysis_type);
     29        SpcNodesx(nodes,constraints,parameters,analysis_type);
    2830
    2931        /*write output datasets: */
    … …  
    3335        delete nodes;
    3436        delete constraints;
     37        delete parameters;
    3538
    3639        /*end module: */
    … …  
    4144{
    4245        _printf_(true,"\n");
    43         _printf_(true,"   usage: [m.node]=%s(m.nodes,m.constraints);\n",__FUNCT__);
     46        _printf_(true,"   usage: [m.node]=%s(m.nodes,m.constraints,m.parameters);\n",__FUNCT__);
    4447        _printf_(true,"\n");
    4548}

  • issm/trunk/src/mex/SpcNodes/SpcNodes.h

    r5778 r9002  
    2121#define NODESIN (mxArray*)prhs[0]
    2222#define CONSTRAINTS (mxArray*)prhs[1]
    23 #define ANALYSISTYPE (mxArray*)prhs[2]
     23#define PARAMETERS (mxArray*)prhs[2]
     24#define ANALYSISTYPE (mxArray*)prhs[3]
    2425
    2526/* serial output macros: */
    … …  
    3031#define NLHS  1
    3132#undef NRHS
    32 #define NRHS  3
     33#define NRHS  4
    3334
    3435
Note: See TracChangeset for help on using the changeset viewer.