Changeset 16164
- Timestamp:
- 09/18/13 09:17:13 (12 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 32 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r16159 r16164 433 433 434 434 /*retrieve parameters: */ 435 ElementMatrix* Ke=NULL;436 435 int analysis_type; 437 436 parameters->FindParam(&analysis_type,AnalysisTypeEnum); … … 1413 1412 /*Recover input*/ 1414 1413 Input* input=inputs->GetInput(enumtype); 1415 if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); 1416 int numnodes = this->NumberofNodes(); 1414 if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); 1417 1415 1418 1416 /* Start looping on the number of vertices: */ 1419 1417 GaussPenta* gauss=new GaussPenta(); 1420 for 1418 for(int iv=0;iv<this->NumberofNodes();iv++){ 1421 1419 gauss->GaussNode(this->element_type,iv); 1422 1420 input->GetInputValue(&pvalue[iv],gauss); … … 4084 4082 /*Intermediaries */ 4085 4083 int stabilization; 4086 int i,j ,found=0;4084 int i,j; 4087 4085 IssmDouble Jdet,u,v,w,um,vm,wm; 4088 4086 IssmDouble h,hx,hy,hz,vx,vy,vz,vel; … … 4320 4318 /*Intermediaries */ 4321 4319 int stabilization; 4322 int i,j ,found=0;4320 int i,j; 4323 4321 IssmDouble Jdet,u,v,w,um,vm,wm,vel; 4324 4322 IssmDouble h,hx,hy,hz,vx,vy,vz; … … 4540 4538 4541 4539 /*Intermediaries*/ 4542 int i ,found=0;4540 int i; 4543 4541 int stabilization; 4544 4542 IssmDouble Jdet,phi,dt; … … 4550 4548 IssmDouble u,v,w; 4551 4549 IssmDouble scalar_def,scalar_transient; 4552 IssmDouble temperature_list[NUMVERTICES];4553 4550 IssmDouble xyz_list[NUMVERTICES][3]; 4554 4551 IssmDouble L[numdof]; … … 4809 4806 4810 4807 /*Intermediaries*/ 4811 int i ,found=0;4808 int i; 4812 4809 int stabilization; 4813 4810 IssmDouble Jdet,phi,dt; … … 4818 4815 IssmDouble u,v,w; 4819 4816 IssmDouble scalar_def,scalar_transient; 4820 IssmDouble temperature_list[NUMVERTICES];4821 4817 IssmDouble xyz_list[NUMVERTICES][3]; 4822 4818 IssmDouble L[numdof]; … … 5086 5082 IssmDouble B_average,s_average; 5087 5083 int *doflist = NULL; 5088 IssmDouble pressure[numdof];5084 //IssmDouble pressure[numdof]; 5089 5085 5090 5086 /*Get dof list: */ … … 7918 7914 int i,j; 7919 7915 IssmDouble Jdet,viscosity; 7920 IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/7921 7916 IssmDouble xyz_list[NUMVERTICES][3]; 7922 7917 IssmDouble B[3][numdof2d]; … … 8658 8653 IssmDouble xyz_list[NUMVERTICES][3]; 8659 8654 IssmDouble basis[6]; //for the six nodes of the penta 8660 Tria* tria=NULL;8661 8655 Friction* friction=NULL; 8662 8656 GaussPenta *gauss=NULL; … … 8821 8815 IssmDouble xyz_list[NUMVERTICES][3]; 8822 8816 IssmDouble basis[6]; //for the six nodes of the penta 8823 Tria* tria=NULL;8824 8817 Friction* friction=NULL; 8825 8818 GaussPenta *gauss=NULL; … … 9367 9360 IssmDouble xyz_list[NUMVERTICES][3]; 9368 9361 IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ 9369 IssmDouble l1l6[6]; //for the six nodes and the bubble9370 IssmDouble dh1dh6[3][NUMVERTICES];9371 9362 GaussPenta *gauss=NULL; 9372 9363 9373 9364 /*Stabilization*/ 9374 bool stabilization = true;9375 9365 IssmDouble dbasis[3][6]; 9376 9366 IssmDouble dmu[3]; … … 9412 9402 9413 9403 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); 9414 GetNodalFunctionsP1(&l1l6[0], gauss);9415 9404 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 9416 9405 material->GetViscosity3dFS(&viscosity,&epsilon[0]); … … 9790 9779 /*FUNCTION Penta::CreateJacobianStressbalanceHO{{{*/ 9791 9780 ElementMatrix* Penta::CreateJacobianStressbalanceHO(void){ 9792 9793 /*Constants*/9794 const int numdof=NDOF2*NUMVERTICES;9795 9781 9796 9782 /*Intermediaries */ … … 10670 10656 int* doflists = NULL; 10671 10657 int* doflistpressure = NULL; 10672 Penta *penta = NULL;10673 10674 /*OK, we have to add results of this element for HO10675 * and results from the penta at base for SSA. Now recover results*/10676 penta=GetBasalElement();10677 10658 10678 10659 /*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */ -
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r16162 r16164 350 350 351 351 /*retrive parameters: */ 352 ElementVector* pe=NULL;353 352 int analysis_type; 354 353 parameters->FindParam(&analysis_type,AnalysisTypeEnum); … … 1225 1224 /*Recover input*/ 1226 1225 Input* input=inputs->GetInput(enumtype); 1227 if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); 1228 int numnodes = this->NumberofNodes(); 1226 if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); 1229 1227 1230 1228 /* Start looping on the number of vertices: */ … … 2041 2039 2042 2040 /*Intermediaries*/ 2043 int i; 2044 int numberofresults = 0; 2045 int *resultsenums = NULL; 2046 int *resultssizes = NULL; 2047 IssmDouble *resultstimes = NULL; 2048 int *resultssteps = NULL; 2041 int *resultsenums = NULL; 2042 int *resultssizes = NULL; 2043 IssmDouble *resultstimes = NULL; 2044 int *resultssteps = NULL; 2049 2045 2050 2046 /*Checks*/ … … 2052 2048 2053 2049 /*Count number of results*/ 2054 for(i=0;i<this->results->Size();i++){ 2055 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); 2056 numberofresults++; 2057 } 2050 int numberofresults = this->results->Size(); 2058 2051 2059 2052 if(numberofresults){ … … 2066 2059 2067 2060 /*populate enums*/ 2068 for(i =0;i<this->results->Size();i++){2061 for(int i=0;i<this->results->Size();i++){ 2069 2062 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); 2070 2063 resultsenums[i]=elementresult->InstanceEnum(); … … 2420 2413 /*Intermediaries*/ 2421 2414 int i; 2422 int tria_vertex_ids[3];2423 2415 IssmDouble nodeinputs[3]; 2424 2416 IssmDouble yts; … … 2441 2433 /*Recover element type*/ 2442 2434 this->SetElementType(finiteelement_type,analysis_counter); 2443 2444 /*Recover vertices ids needed to initialize inputs*/2445 for(i=0;i<3;i++){2446 tria_vertex_ids[i]=reCast<int>(iomodel->elements[3*index+i]); //ids for vertices are in the elements array from Matlab2447 }2448 2435 2449 2436 /*Recover nodes ids needed to initialize the node hook.*/ … … 3181 3168 /*Fetch number of nodes and dof for this finite element*/ 3182 3169 int numnodes = this->NumberofNodes(); 3183 int numdof = numnodes*NDOF2;3184 3170 3185 3171 /*Initialize Element vector and vectors*/ … … 3314 3300 /*Fetch number of nodes and dof for this finite element*/ 3315 3301 int numnodes = this->NumberofNodes(); _assert_(numnodes==3); 3316 int numdof = numnodes*NDOF2;3317 3302 3318 3303 /*Initialize Element vector*/ … … 3371 3356 /*Fetch number of nodes and dof for this finite element*/ 3372 3357 int numnodes = this->NumberofNodes(); 3373 int numdof = numnodes*NDOF2;3374 3358 3375 3359 /*Initialize Element matrix, vectors and Gaussian points*/ … … 4033 4017 IssmDouble grade_g_gaussian[NUMVERTICES]; 4034 4018 IssmDouble basis[3]; 4035 IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/4036 4019 Friction* friction=NULL; 4037 4020 GaussTria *gauss=NULL; … … 4537 4520 IssmDouble Tria::SurfaceLogVxVyMisfit(int weight_index){ 4538 4521 4539 int fit=-1;4540 4522 IssmDouble Jelem=0, S=0; 4541 4523 IssmDouble epsvel=2.220446049250313e-16; … … 7258 7240 IssmDouble D_scalar,Jdet,thickness; 7259 7241 IssmDouble xyz_list[NUMVERTICES][3]; 7260 IssmDouble D[2][2];7261 7242 IssmDouble l=8.; 7262 7243 … … 7498 7479 ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); 7499 7480 IssmDouble* basis = xNew<IssmDouble>(numnodes); 7500 IssmDouble* Vx = xNew<IssmDouble>(numnodes);7501 IssmDouble* Vy = xNew<IssmDouble>(numnodes);7502 7481 7503 7482 /*Retrieve all inputs and parameters*/ -
issm/trunk-jpl/src/c/classes/FemModel.cpp
r16144 r16164 400 400 void FemModel::UpdateConstraintsx(void){ /*{{{*/ 401 401 402 Element *element = NULL;403 402 IssmDouble time; 404 403 int analysis_type; … … 1022 1021 1023 1022 /*Intermediary*/ 1024 int num_responses; 1025 Element *element = NULL; 1026 int *responses = NULL; 1023 int num_responses; 1024 int *responses = NULL; 1027 1025 1028 1026 /*output: */ -
issm/trunk-jpl/src/c/classes/IoModel.cpp
r16141 r16164 288 288 va_list ap; 289 289 int dataenum; 290 int i;291 DoubleMatParam* parameter=NULL;292 290 293 291 /*Go through the entire list of enums and delete the corresponding data from the iomodel-data dataset: */ 294 295 292 va_start(ap,num); 296 for(i = 0; i <num;i++){297 dataenum=va_arg(ap, 293 for(int i=0;i<num;i++){ 294 dataenum=va_arg(ap,int); 298 295 _assert_(dataenum<MaximumNumberOfDefinitionsEnum); 299 296 … … 1041 1038 int i; 1042 1039 bool defaulting = false; 1043 bool transient = false;1044 1045 FILE *fid = NULL;1046 1040 int code = 0; 1047 1041 int vector_layout = 0; … … 1061 1055 1062 1056 /*First of, find the record for the enum, and get code of data type: */ 1063 fid=this->SetFilePointerToData(&code, &vector_layout,vector_enum);1057 this->SetFilePointerToData(&code, &vector_layout,vector_enum); 1064 1058 1065 1059 switch(code){ -
issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp
r16042 r16164 39 39 int pos1,pos2,pos3,pos4; 40 40 int num_nodes; 41 int num_elems;42 41 43 42 /*numericalflux constructor data: */ 44 43 int numericalflux_elem_ids[2]; 45 int numericalflux_mparid;46 44 int numericalflux_vertex_ids[2]; 47 45 int numericalflux_node_ids[4]; 48 46 int numericalflux_type; 49 50 /* Get MatPar id */51 numericalflux_mparid=iomodel->numberofelements+1; //matlab indexing52 47 53 48 /*Get edge*/ … … 60 55 if(e2==-1){ 61 56 /* Boundary edge, only one element */ 62 num_ elems=1; num_nodes=2;57 num_nodes=2; 63 58 numericalflux_type=BoundaryEnum; 64 59 numericalflux_elem_ids[0]=e1; … … 66 61 else{ 67 62 /* internal edge: connected to 2 elements */ 68 num_elems=2;num_nodes=4;63 num_nodes=4; 69 64 numericalflux_type=InternalEnum; 70 65 numericalflux_elem_ids[0]=e1; -
issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
r16151 r16164 441 441 442 442 // The penalty is stable if it doesn't change during to successive iterations. 443 444 int found=0;445 443 const int numnodes=1; 446 444 IssmDouble pressure; … … 449 447 int new_active; 450 448 int unstable=0; 451 int reset_penalties=0;452 449 int penalty_lock; 453 450 … … 508 505 ElementMatrix* Pengrid::PenaltyCreateKMatrixMelting(IssmDouble kmax){ 509 506 510 const int numdof=NUMVERTICES*NDOF1;511 507 IssmDouble pressure,temperature,t_pmp; 512 508 IssmDouble penalty_factor; … … 538 534 ElementMatrix* Pengrid::PenaltyCreateKMatrixThermal(IssmDouble kmax){ 539 535 540 const int numdof=NUMVERTICES*NDOF1;541 536 IssmDouble penalty_factor; 542 537 … … 557 552 ElementVector* Pengrid::PenaltyCreatePVectorMelting(IssmDouble kmax){ 558 553 559 const int numdof=NUMVERTICES*NDOF1;560 554 IssmDouble pressure; 561 555 IssmDouble temperature; … … 600 594 ElementVector* Pengrid::PenaltyCreatePVectorThermal(IssmDouble kmax){ 601 595 602 const int numdof=NUMVERTICES*NDOF1;603 596 IssmDouble pressure; 604 597 IssmDouble t_pmp; … … 632 625 const int numnodes = 1; 633 626 int unstable = 0; 634 int reset_penalties = 0;635 int found = 0;636 627 int new_active; 637 628 IssmDouble pressure; … … 669 660 ElementMatrix* Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax){ 670 661 671 const int numdof=NUMVERTICES*NDOF1;672 662 IssmDouble penalty_factor; 673 663 … … 688 678 ElementVector* Pengrid::PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax){ 689 679 690 const int numdof=NUMVERTICES*NDOF1;691 680 IssmDouble h_max; 692 681 IssmDouble penalty_factor; -
issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp
r16125 r16164 454 454 455 455 const int numdof = NDOF2*NUMVERTICES; 456 int dofs[1] = {0};457 456 IssmDouble thickness; 458 457 IssmDouble h[2]; … … 534 533 ElementVector* Riftfront::PenaltyCreatePVectorStressbalanceHoriz(IssmDouble kmax){ 535 534 536 const int numdof = NDOF2*NUMVERTICES;537 535 int j; 538 536 IssmDouble rho_ice; … … 636 634 int Riftfront::Constrain(int* punstable){ 637 635 638 const int numnodes = 2; 639 IssmDouble penetration; 636 IssmDouble penetration; 640 637 int activate; 641 638 int unstable; 642 IssmDouble 643 IssmDouble 644 IssmDouble 645 IssmDouble 646 IssmDouble 639 IssmDouble vx1; 640 IssmDouble vy1; 641 IssmDouble vx2; 642 IssmDouble vy2; 643 IssmDouble fractionincrement; 647 644 648 645 /*Objects: */ 649 Tria *tria1= NULL;650 Tria *tria2= NULL;646 Tria *tria1 = NULL; 647 Tria *tria2 = NULL; 651 648 652 649 /*enum of element? */ … … 756 753 int Riftfront::MaxPenetration(IssmDouble* ppenetration){ 757 754 758 const int numnodes=2; 759 IssmDouble penetration=0; 760 IssmDouble vx1; 761 IssmDouble vy1; 762 IssmDouble vx2; 763 IssmDouble vy2; 755 IssmDouble penetration; 756 IssmDouble vx1; 757 IssmDouble vy1; 758 IssmDouble vx2; 759 IssmDouble vy2; 764 760 765 761 /*Objects: */ 766 Tria *tria1= NULL;767 Tria *tria2= NULL;762 Tria *tria1 = NULL; 763 Tria *tria2 = NULL; 768 764 769 765 /*enum of element? */ … … 773 769 tria1=(Tria*)elements[0]; 774 770 tria2=(Tria*)elements[1]; 775 776 //initialize:777 penetration=-1;778 771 779 772 /*recover velocity: */ … … 787 780 788 781 /*Now, we return penetration only if we are active!: */ 789 if(this->active==0) penetration=-1;782 if(this->active==0) penetration=-1.; 790 783 791 784 /*If we are zigzag locked, same thing: */ 792 if(this->counter>this->penalty_lock) penetration=-1;785 if(this->counter>this->penalty_lock) penetration=-1.; 793 786 794 787 /*assign output pointer: */ … … 838 831 int Riftfront::PotentialUnstableConstraint(int* punstable){ 839 832 840 const int numnodes = 2;841 833 IssmDouble penetration; 842 834 int unstable; … … 887 879 int Riftfront::PreConstrain(int* punstable){ 888 880 889 const int numnodes = 2; 890 IssmDouble penetration; 881 IssmDouble penetration; 891 882 int unstable; 892 IssmDouble 893 IssmDouble 894 IssmDouble 895 IssmDouble 883 IssmDouble vx1; 884 IssmDouble vy1; 885 IssmDouble vx2; 886 IssmDouble vy2; 896 887 897 888 /*Objects: */ 898 Tria *tria1= NULL;899 Tria *tria2= NULL;889 Tria *tria1 = NULL; 890 Tria *tria2 = NULL; 900 891 901 892 /*enum of element? */ -
issm/trunk-jpl/src/c/classes/Node.cpp
r16125 r16164 26 26 /*Intermediary*/ 27 27 int k,l; 28 int gsize;29 28 30 29 /*id: */ … … 41 40 this->indexingupdate = true; 42 41 DistributeNumDofs(&this->indexing,analysis_type,in_approximation); //number of dofs per node 43 gsize=this->indexing.gsize;44 42 45 43 if(analysis_type==StressbalanceAnalysisEnum) … … 1128 1126 void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array){/*{{{*/ 1129 1127 1130 int i,counter;1131 int numdofs= 0;1128 int i,counter; 1129 int numdofs = 0; 1132 1130 IssmDouble norm; 1133 IssmDouble *transform = NULL; 1134 IssmDouble *values = NULL; 1131 IssmDouble *transform = NULL; 1135 1132 IssmDouble coord_system[3][3]; 1136 1133 -
issm/trunk-jpl/src/c/classes/kriging/Observations.cpp
r15741 r16164 125 125 /*Output and Intermediaries*/ 126 126 int nobs,i,index; 127 IssmPDouble hmin,h2,hmin2 ,radius2;127 IssmPDouble hmin,h2,hmin2; 128 128 int *indices = NULL; 129 129 Observation *observation = NULL; … … 139 139 if(hmin<radius) radius=hmin; 140 140 } 141 142 /*Compute radius square*/143 radius2 = radius*radius;144 141 145 142 /*Find all observations that are in radius*/ -
issm/trunk-jpl/src/c/classes/kriging/Quadtree.cpp
r15741 r16164 257 257 QuadtreeBox *box = NULL; 258 258 int xi,yi; 259 int level ,levelbin;259 int levelbin; 260 260 int index = -1; 261 261 double length,length2; … … 264 264 this->IntergerCoordinates(&xi,&yi,x,y); 265 265 266 /*Initialize levels*/ 267 level = 0; 266 /*Initialize level*/ 268 267 levelbin = (1L<<this->MaxDepth);// = 2^30 269 268 … … 273 272 /*Find the smallest box where this point is located*/ 274 273 while((box=*pbox) && (box->nbitems<0)){ 275 276 levelbin>>=1; level+=1; 277 274 levelbin>>=1; 278 275 pbox = &box->box[IJ(xi,yi,levelbin)]; 279 276 } -
issm/trunk-jpl/src/c/main/kriging.cpp
r15839 r16164 14 14 FILE *output_fid = NULL; 15 15 FILE *input_fid = NULL; 16 bool waitonlock = false;17 16 18 17 /*File names*/ -
issm/trunk-jpl/src/c/modules/AverageOntoPartitionx/AverageOntoPartitionx.cpp
r14999 r16164 19 19 int npart; 20 20 double *qmu_part = NULL; 21 int numberofvertices;22 21 23 22 /*output: */ … … 31 30 32 31 /*Some parameters: */ 33 numberofvertices=vertices->NumberOfVertices();34 32 parameters->FindParam(&npart,QmuNumberofpartitionsEnum); 35 33 -
issm/trunk-jpl/src/c/modules/ComputeBasalStressx/ComputeBasalStressx.cpp
r14999 r16164 11 11 12 12 /*Intermediary*/ 13 int i; 14 int found=0; 15 IssmDouble numberofelements; 16 Element* element=NULL; 13 IssmDouble numberofelements; 14 Element* element=NULL; 17 15 18 16 /*output: */ … … 26 24 27 25 /*Compute basal stress for each element: */ 28 for (i=0;i<elements->Size();i++){26 for(int i=0;i<elements->Size();i++){ 29 27 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i)); 30 28 element->ComputeBasalStress(sigma); -
issm/trunk-jpl/src/c/modules/ComputeStrainRatex/ComputeStrainRatex.cpp
r14999 r16164 8 8 #include "../../toolkits/toolkits.h" 9 9 10 void ComputeStrainRatex( 10 void ComputeStrainRatex(Vector<IssmDouble>** peps,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,Parameters* parameters){ 11 11 12 12 /*Intermediary*/ 13 int i;14 int found=0;15 13 int numberofelements; 16 14 Element* element=NULL; … … 26 24 27 25 /*Compute basal stress for each element: */ 28 for (i=0;i<elements->Size();i++){26 for(int i=0;i<elements->Size();i++){ 29 27 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i)); 30 28 element->ComputeStrainRate(eps); -
issm/trunk-jpl/src/c/modules/ConfigureObjectsx/ConfigureObjectsx.cpp
r15439 r16164 12 12 13 13 /*Intermediary*/ 14 int i; 15 int noerr=1; 16 int configuration_type; 17 Element *element = NULL; 18 Load *load = NULL; 19 Node *node = NULL; 20 Material *material = NULL; 14 int i; 15 int noerr = 1; 16 int configuration_type; 17 Element *element = NULL; 18 Load *load = NULL; 19 Material *material = NULL; 21 20 22 21 /*Get analysis type: */ … … 24 23 25 24 if(VerboseMProcessor()) _printf0_(" Configuring elements...\n"); 26 for 25 for(i=0;i<elements->Size();i++){ 27 26 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i)); 28 27 element->Configure(elements,loads,nodes,vertices,materials,parameters); 29 28 } 30 29 if(VerboseMProcessor()) _printf0_(" Configuring loads...\n"); 31 for 30 for(i=0;i<loads->Size();i++){ 32 31 load=(Load*)loads->GetObjectByOffset(i); 33 32 if (load->InAnalysis(configuration_type)){ … … 36 35 } 37 36 if(VerboseMProcessor()) _printf0_(" Configuring materials...\n"); 38 for 37 for(i=0;i<materials->Size();i++){ 39 38 material=(Material*)materials->GetObjectByOffset(i); 40 39 material->Configure(elements); -
issm/trunk-jpl/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp
r15838 r16164 230 230 void RiftSetPreStable(Loads* loads){ 231 231 232 int i;233 234 Riftfront* riftfront=NULL;235 int found=0;236 int mpi_found=0;237 238 232 /*go though loads, and set loads to pre stable.:*/ 239 for (i=0;i<loads->Size();i++){ 240 233 for(int i=0;i<loads->Size();i++){ 241 234 if(RiftfrontEnum==loads->GetEnum(i)){ 242 243 riftfront=(Riftfront*)loads->GetObjectByOffset(i); 235 Riftfront* riftfront=(Riftfront*)loads->GetObjectByOffset(i); 244 236 riftfront->SetPreStable(); 245 237 } -
issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp
r15335 r16164 13 13 void* ContourToMeshxt(void* vpthread_handle){ 14 14 15 int noerr=1;16 17 15 /*gate variables :*/ 18 ContourToMeshxThreadStruct* gate=NULL; 19 pthread_handle* handle=NULL; 20 int my_thread; 21 int num_threads; 22 int i0; 23 int i1; 24 25 int i; 26 27 /*Contour:*/ 28 Contours* contours=NULL; 29 30 /*parameters: */ 31 int nods; 32 int edgevalue; 33 double* x=NULL; 34 double* y=NULL; 35 double* in_nod=NULL; 16 ContourToMeshxThreadStruct *gate = NULL; 17 pthread_handle *handle = NULL; 18 int i,i1,i0; 36 19 37 20 /*recover handle and gate: */ 38 handle =(pthread_handle*)vpthread_handle;39 gate =(ContourToMeshxThreadStruct*)handle->gate;40 my_thread=handle->id;41 num_threads=handle->num;21 handle = (pthread_handle*)vpthread_handle; 22 gate = (ContourToMeshxThreadStruct*)handle->gate; 23 int my_thread = handle->id; 24 int num_threads = handle->num; 42 25 43 26 /*recover parameters :*/ 44 contours=gate->contours;45 nods=gate->nods;46 edgevalue=gate->edgevalue;47 in_nod=gate->in_nod;48 x=gate->x;49 y=gate->y;27 Contours* contours = gate->contours; 28 int nods = gate->nods; 29 int edgevalue = gate->edgevalue; 30 double *in_nod = gate->in_nod; 31 double *x = gate->x; 32 double *y = gate->y; 50 33 51 34 /*distribute indices across threads :*/ -
issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.cpp
r15335 r16164 34 34 int ContourToNodesx(IssmPDouble** pflags,double* x, double* y, int nods, Contours* contours, int edgevalue){ 35 35 36 /*Contour:*/37 Contour<IssmPDouble>* contouri=NULL;38 double* xc=NULL;39 double* yc=NULL;40 41 36 /*output: */ 42 37 IssmPDouble* flags=NULL; -
issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp
r16126 r16164 11 11 12 12 int i,connectivity; 13 int numberofdofspernode; 14 int fsize,configuration_type; 13 int configuration_type; 15 14 Element *element = NULL; 16 15 Load *load = NULL; … … 23 22 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); 24 23 femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum); 25 fsize=femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum);26 numberofdofspernode=femmodel->nodes->MaxNumDofs(configuration_type,GsetEnum);27 24 28 25 /*Initialize Jacobian Matrix*/ -
issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp
r15140 r16164 92 92 93 93 /*Recover material parameters and loading history: see GiaDeflectionCoreArgs for more details {{{*/ 94 ri =arguments->ri;95 re =arguments->re;96 hes =arguments->hes;97 times =arguments->times;98 numtimes =arguments->numtimes;99 currenttime =arguments->currenttime;100 lithosphere_shear_modulus =arguments->lithosphere_shear_modulus;101 lithosphere_density =arguments->lithosphere_density;102 mantle_shear_modulus =arguments->mantle_shear_modulus;103 mantle_viscosity =arguments->mantle_viscosity;104 mantle_density =arguments->mantle_density;105 lithosphere_thickness =arguments->lithosphere_thickness;106 rho_ice =arguments->rho_ice;107 disk_id =arguments->idisk;108 iedge =arguments->iedge;109 yts =arguments->yts;94 ri = arguments->ri; 95 re = arguments->re; 96 hes = arguments->hes; 97 times = arguments->times; 98 numtimes = arguments->numtimes; 99 currenttime = arguments->currenttime; 100 lithosphere_shear_modulus = arguments->lithosphere_shear_modulus; 101 lithosphere_density = arguments->lithosphere_density; 102 mantle_shear_modulus = arguments->mantle_shear_modulus; 103 mantle_viscosity = arguments->mantle_viscosity; 104 mantle_density = arguments->mantle_density; 105 lithosphere_thickness = arguments->lithosphere_thickness; 106 rho_ice = arguments->rho_ice; 107 disk_id = arguments->idisk; 108 iedge = arguments->iedge; 109 yts = arguments->yts; 110 110 111 111 /*}}}*/ -
issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp
r15969 r16164 80 80 /*FUNCTION PropagateFloatingiceToGroundedNeighbors {{{*/ 81 81 IssmDouble* PropagateFloatingiceToGroundedNeighbors(Elements* elements,Nodes* nodes,Vertices* vertices,Parameters* parameters,IssmDouble* vertices_potentially_ungrounding){ 82 int i,analysis_type ,numberofvertices;82 int i,analysis_type; 83 83 int nflipped,local_nflipped; 84 84 IssmDouble* phi = NULL; … … 86 86 Vector<IssmDouble>* vec_elements_neighboring_floatingice = NULL; 87 87 Vector<IssmDouble>* vec_phi = NULL; 88 Node* node = NULL;89 88 Element* element = NULL; 90 89 91 90 /*recover parameters: */ 92 91 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 93 numberofvertices=vertices->NumberOfVertices();94 92 95 93 /*recover vec_phi*/ -
issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp
r15557 r16164 134 134 135 135 bool debug = M*N>1? true:false; 136 debug = true;137 136 138 137 /*partition loop across threads: */ -
issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp
r15694 r16164 11 11 12 12 /*intermediary: */ 13 FILE *fid = NULL;14 13 int code,vector_layout; 15 14 IssmDouble *spcdata = NULL; … … 17 16 18 17 /*First of, find the record for the enum, and get code of data type: */ 19 fid=iomodel->SetFilePointerToData(&code, &vector_layout,vector_enum);18 iomodel->SetFilePointerToData(&code, &vector_layout,vector_enum); 20 19 if(code!=7)_error_("expecting a IssmDouble vector for constraints with enum " << EnumToStringx(vector_enum)); 21 20 if(vector_layout!=1)_error_("expecting a nodal vector for constraints with enum " << EnumToStringx(vector_enum)); -
issm/trunk-jpl/src/c/modules/Kml2Expx/Kml2Expx.cpp
r15104 r16164 8 8 #include "../KMLFileReadx/KMLFileReadx.h" 9 9 10 int Kml2Expx(char* filkml,char* filexp, 11 int sgn){ 10 int Kml2Expx(char* filkml,char* filexp,int sgn){ 12 11 13 12 double cm,sp; 14 15 13 Ll2xydef(&cm,&sp,sgn); 16 14 17 return(Kml2Expx(filkml,filexp, 18 sgn,cm,sp)); 15 return(Kml2Expx(filkml,filexp,sgn,cm,sp)); 19 16 } 20 17 21 int Kml2Expx(char* filkml,char* filexp, 22 int sgn,double cm,double sp){ 18 int Kml2Expx(char* filkml,char* filexp,int sgn,double cm,double sp){ 23 19 24 int iret=0; 25 double *lat=NULL,*lon=NULL; 26 27 KML_Object* kobj=NULL; 28 29 FILE* fidi=NULL; 30 FILE* fido=NULL; 31 32 clock_t clock0,clock1; 33 time_t time0, time1; 20 int iret = 0; 21 KML_Object *kobj = NULL; 22 FILE *fidi = NULL; 23 FILE *fido = NULL; 24 clock_t clock0,clock1; 25 time_t time0 ,time1; 34 26 35 27 clock0=clock(); … … 37 29 _printf0_("\nKml2Expx Module -- " << ctime(&time0)); 38 30 39 /* read kml file */ 40 31 /*read kml file*/ 41 32 fidi=fopen(filkml,"r"); 42 33 if (!(kobj=KMLFileReadx(fidi))) 43 34 _error_("Error reading kml file."); 44 35 fclose(fidi); 45 36 46 /* open exp file */ 47 37 /*open exp file*/ 48 38 _printf0_("Writing exp profiles to file.\n"); 49 39 fido=fopen(filexp,"w"); 50 40 51 /* write the polygons and linestrings */ 52 41 /*write the polygons and linestrings */ 53 42 kobj->WriteExp(fido,"",sgn,cm,sp); 54 43 55 /* close exp file */ 56 44 /*close exp file */ 57 45 fclose(fido); 58 59 46 delete kobj; 60 47 -
issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp
r14219 r16164 75 75 76 76 /*output: */ 77 double* segments=NULL; 78 Segment<double>* segment=NULL; 77 double* segments=NULL; 79 78 int numsegs; 80 79 -
issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
r16144 r16164 188 188 temp_m=mdims_array[i]; 189 189 temp_n=ndims_array[i]; 190 _assert_(temp_n==5); 190 191 191 192 m=0; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/CreateConstraintsStressbalance.cpp
r15986 r16164 85 85 } 86 86 } 87 else{ 88 _error_("model not supported yet"); 89 } 87 90 88 91 IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceAnalysisEnum,finiteelement,1); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/CreateLoadsStressbalance.cpp
r15986 r16164 9 9 void CreateLoadsStressbalance(Loads** ploads, IoModel* iomodel){ 10 10 11 /*DataSets*/12 Pengrid *pengrid = NULL;13 14 11 /*Intermediary*/ 15 12 int segment_width; … … 19 16 int penpair_ids[2]; 20 17 bool isSSA,isL1L2,isHO,isFS; 21 int numpenalties,numberofpressureloads,numrifts,numriftsegments; 22 IssmDouble *pressureload = NULL; 23 IssmDouble *elements_type = NULL; 24 IssmDouble *nodeoniceshelf = NULL; 18 int numpenalties,numrifts,numriftsegments; 25 19 IssmDouble *riftinfo = NULL; 26 IssmDouble *nodeonbed = NULL;27 IssmDouble *nodeonFS = NULL;28 IssmDouble *nodeonicesheet = NULL;29 IssmDouble *vertices_type = NULL;30 20 IssmDouble *penalties = NULL; 31 21 int assert_int; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/StressbalanceSIA/CreateConstraintsStressbalanceSIA.cpp
r15771 r16164 17 17 /*Output*/ 18 18 Constraints* constraints = NULL; 19 SpcStatic* spcstatic = NULL;20 19 21 20 /*Recover pointer: */ -
issm/trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp
r15838 r16164 10 10 void NodalValuex( IssmDouble* pnodalvalue, int natureofdataenum,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){ 11 11 12 int my_rank;13 int i;14 int index;15 Element* element=NULL;16 12 IssmDouble value; 17 int found; 18 int sumfound; 19 int cpu_found; 20 21 /*Get my_rank:*/ 22 my_rank=IssmComm::GetRank(); 13 int index; 14 int found,sumfound,cpu_found; 23 15 24 16 /*retrieve element we are interested in: */ … … 27 19 /*This is the vertex id for which we want to collect the data. Go through elements, and for each 28 20 *element, figure out if they hold the vertex, and the data. If so, return it: */ 29 for(i =0;i<elements->Size();i++){21 for(int i=0;i<elements->Size();i++){ 30 22 Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i)); 31 23 found=element->NodalValue(&value,index,natureofdataenum); 32 if 33 cpu_found= my_rank;24 if(found){ 25 cpu_found=IssmComm::GetRank(); 34 26 break; 35 27 } … … 37 29 38 30 /*Broadcast whether we found the element: */ 39 ISSM_MPI_Allreduce (&found,&sumfound,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());31 ISSM_MPI_Allreduce(&found,&sumfound,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm()); 40 32 if(!sumfound)_error_("could not find element with vertex with id" << index << " to compute nodal value " << EnumToStringx(natureofdataenum)); 41 33 -
issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp
r14656 r16164 12 12 /*threading: */ 13 13 PointCloudFindNeighborsThreadStruct gate; 14 int num=1;15 14 16 15 #ifdef _MULTITHREADING_ 17 num=_NUMTHREADS_; 16 int num=_NUMTHREADS_; 17 #else 18 int num=1; 18 19 #endif 19 20 -
issm/trunk-jpl/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp
r14951 r16164 16 16 IssmDouble delta18oLapseRate=-6.2*pow(10.,-3); 17 17 IssmDouble glacialindex; // used to vary present day temperature 18 int imonth=0;19 18 20 glacialindex = 0 ;//(Delta18oTime-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceTime-Delta18oSurfacePresent))19 glacialindex = 0.;//(Delta18oTime-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceTime-Delta18oSurfacePresent)) 21 20 // /(Delta18oLgm-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceLgm-Delta18oSurfacePresent)); 22 21 23 22 for (int imonth = 0; imonth<12; imonth++){ 24 monthlytemperaturestmp[imonth] = glacialindex*TemperaturesLgm[imonth] + (1 -glacialindex)*TemperaturesPresentday[imonth];23 monthlytemperaturestmp[imonth] = glacialindex*TemperaturesLgm[imonth] + (1.-glacialindex)*TemperaturesPresentday[imonth]; 25 24 //monthlyprectmp[imonth] = 1.5*pow(2,((monthlytemperaturestmp[imonth]-273.15-0)/10)); //equation from rob's paper 26 25 monthlyprectmp[imonth] = PrecipitationsPresentday[imonth];
Note:
See TracChangeset
for help on using the changeset viewer.