Changeset 13622
- Timestamp:
- 10/11/12 11:23:47 (12 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 288 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/Container/Constraints.cpp
r13604 r13622 54 54 numberofconstraints=localconstraints; 55 55 #endif 56 57 56 58 57 return numberofconstraints; -
issm/trunk-jpl/src/c/Container/DataSet.cpp
r13056 r13622 30 30 /*FUNCTION DataSet::DataSet(){{{*/ 31 31 DataSet::DataSet(){ 32 32 33 33 sorted=0; 34 34 sorted_ids=NULL; … … 40 40 DataSet::DataSet(int dataset_enum){ 41 41 enum_type=dataset_enum; 42 42 43 43 sorted=0; 44 44 sorted_ids=NULL; … … 130 130 void DataSet::DeepEcho(){ 131 131 132 133 132 vector<Object*>::iterator object; 134 133 -
issm/trunk-jpl/src/c/Container/Elements.cpp
r13612 r13622 65 65 /*FUNCTION Elements::DeleteResults{{{*/ 66 66 void Elements::DeleteResults(void){ 67 67 68 68 for (int i=0;i<this->Size();i++){ 69 69 Element* element=(Element*)this->GetObjectByOffset(i); … … 108 108 * We will use the Patch object, which will store all of the information needed, and will be able 109 109 * to output itself to disk on its own. See object/Patch.h for format of this object.*/ 110 110 111 111 /*First, determine maximum number of vertices, nodes, and number of results: */ 112 112 numrows=0; … … 167 167 168 168 int i; 169 169 170 170 int my_rank; 171 171 int num_procs; … … 184 184 int rank; 185 185 int minrank; 186 186 187 187 /*recover my_rank:*/ 188 188 my_rank=IssmComm::GetRank(); -
issm/trunk-jpl/src/c/Container/Inputs.cpp
r13056 r13622 207 207 /*FUNCTION Inputs::ConstrainMin{{{*/ 208 208 void Inputs::ConstrainMin(int constrain_enum, IssmDouble minimum){ 209 209 210 210 Input* constrain_input=NULL; 211 211 /*Find x and y inputs: */ … … 407 407 /*FUNCTION Inputs::AXPY{{{*/ 408 408 void Inputs::AXPY(int MeshYEnum, IssmDouble scalar, int MeshXEnum){ 409 409 410 410 Input* xinput=NULL; 411 411 Input* yinput=NULL; -
issm/trunk-jpl/src/c/Container/Loads.cpp
r13604 r13622 70 70 #endif 71 71 72 73 72 return numberofloads; 74 73 } -
issm/trunk-jpl/src/c/Container/Nodes.cpp
r13612 r13622 66 66 int* alltruedofs=NULL; 67 67 int numnodes=0; 68 68 69 69 /*recover my_rank:*/ 70 70 my_rank=IssmComm::GetRank(); … … 235 235 236 236 int i; 237 237 238 238 int numdofs=0; 239 239 int allnumdofs; … … 271 271 272 272 int i; 273 273 274 274 int numnodes=0; 275 275 int allnumnodes=0; … … 326 326 /*sid starts at 0*/ 327 327 max_sid++; 328 328 329 329 /*return*/ 330 330 return max_sid; … … 337 337 int my_rank; 338 338 int sid; 339 339 340 340 /*recover my_rank:*/ 341 341 my_rank=IssmComm::GetRank(); -
issm/trunk-jpl/src/c/Container/Options.cpp
r13216 r13622 42 42 43 43 char* name=NULL; 44 44 45 45 vector<Object*>::iterator object; 46 46 Option* option=NULL; -
issm/trunk-jpl/src/c/Container/Parameters.cpp
r13425 r13622 53 53 /*FUNCTION Parameters::FindParam(bool* pbool,int enum_type){{{*/ 54 54 void Parameters::FindParam(bool* pbool,int enum_type){ _assert_(this); 55 55 56 56 vector<Object*>::iterator object; 57 57 Param* param=NULL; … … 70 70 /*FUNCTION Parameters::FindParam(int* pinteger,int enum_type){{{*/ 71 71 void Parameters::FindParam(int* pinteger,int enum_type){ _assert_(this); 72 72 73 73 vector<Object*>::iterator object; 74 74 Param* param=NULL; … … 87 87 /*FUNCTION Parameters::FindParam(IssmDouble* pscalar, int enum_type){{{*/ 88 88 void Parameters::FindParam(IssmDouble* pscalar, int enum_type){ _assert_(this); 89 89 90 90 vector<Object*>::iterator object; 91 91 Param* param=NULL; … … 121 121 /*FUNCTION Parameters::FindParam(char** pstring,int enum_type){{{*/ 122 122 void Parameters::FindParam(char** pstring,int enum_type){ _assert_(this); 123 123 124 124 vector<Object*>::iterator object; 125 125 Param* param=NULL; … … 139 139 /*FUNCTION Parameters::FindParam(char*** pstringarray,int* pM,int enum_type){{{*/ 140 140 void Parameters::FindParam(char*** pstringarray,int* pM,int enum_type){ _assert_(this); 141 141 142 142 vector<Object*>::iterator object; 143 143 Param* param=NULL; … … 229 229 /*FUNCTION Parameters::FindParam(IssmDouble*** parray,int* pM,int** pmdims_array,int** pndims_array,int enum_type){{{*/ 230 230 void Parameters::FindParam(IssmDouble*** parray,int* pM,int** pmdims_array,int** pndims_array,int enum_type){ _assert_(this); 231 231 232 232 vector<Object*>::iterator object; 233 233 Param* param=NULL; … … 246 246 /*FUNCTION Parameters::FindParam(Vector<IssmDouble>** pvec,int enum_type){{{*/ 247 247 void Parameters::FindParam(Vector<IssmDouble>** pvec,int enum_type){ _assert_(this); 248 248 249 249 vector<Object*>::iterator object; 250 250 Param* param=NULL; … … 264 264 /*FUNCTION Parameters::FindParam(Matrix<IssmDouble>** pmat,int enum_type){{{*/ 265 265 void Parameters::FindParam(Matrix<IssmDouble>** pmat,int enum_type){ _assert_(this); 266 266 267 267 vector<Object*>::iterator object; 268 268 Param* param=NULL; … … 320 320 321 321 Param* param=NULL; 322 322 323 323 /*first, figure out if the param has already been created: */ 324 324 param=(Param*)this->FindParamObject(enum_type); … … 332 332 333 333 Param* param=NULL; 334 334 335 335 /*first, figure out if the param has already been created: */ 336 336 param=(Param*)this->FindParamObject(enum_type); … … 344 344 345 345 Param* param=NULL; 346 346 347 347 /*first, figure out if the param has already been created: */ 348 348 param=(Param*)this->FindParamObject(enum_type); … … 356 356 357 357 Param* param=NULL; 358 358 359 359 /*first, figure out if the param has already been created: */ 360 360 param=(Param*)this->FindParamObject(enum_type); … … 368 368 369 369 Param* param=NULL; 370 370 371 371 /*first, figure out if the param has already been created: */ 372 372 param=(Param*)this->FindParamObject(enum_type); … … 380 380 381 381 Param* param=NULL; 382 382 383 383 /*first, figure out if the param has already been created: */ 384 384 param=(Param*)this->FindParamObject(enum_type); … … 392 392 393 393 Param* param=NULL; 394 394 395 395 /*first, figure out if the param has already been created: */ 396 396 param=(Param*)this->FindParamObject(enum_type); … … 428 428 429 429 Param* param=NULL; 430 430 431 431 /*first, figure out if the param has already been created: */ 432 432 param=(Param*)this->FindParamObject(enum_type); … … 440 440 441 441 Param* param=NULL; 442 442 443 443 /*first, figure out if the param has already been created: */ 444 444 param=(Param*)this->FindParamObject(enum_type); -
issm/trunk-jpl/src/c/Container/Results.cpp
r12365 r13622 66 66 /*FUNCTION Results::Write{{{*/ 67 67 void Results::Write(Parameters* parameters){ 68 68 69 69 int i; 70 70 FILE *fid = NULL; -
issm/trunk-jpl/src/c/Container/Vertices.cpp
r13612 r13622 49 49 int *truepids = NULL; 50 50 int *alltruepids = NULL; 51 51 52 52 /*recover my_rank:*/ 53 53 my_rank=IssmComm::GetRank(); … … 119 119 int* ranks=NULL; 120 120 int* minranks=NULL; 121 121 122 122 /*recover num_procs:*/ 123 123 num_procs=IssmComm::GetSize(); … … 186 186 int my_rank; 187 187 int sid; 188 188 189 189 /*recover my_rank:*/ 190 190 my_rank=IssmComm::GetRank(); -
issm/trunk-jpl/src/c/classes/DofIndexing.cpp
r13413 r13622 164 164 _printLine_(" ssize: " << ssize); 165 165 _printLine_(" clone: " << clone); 166 166 167 167 _printLine_(" set membership: f,s sets "); 168 168 for(i=0;i<gsize;i++){ -
issm/trunk-jpl/src/c/classes/FemModel.cpp
r13616 r13622 88 88 this->parameters->FindParam(&outbinfilename,OutputFileNameEnum); 89 89 pfclose(output_fid,outbinfilename); 90 90 91 91 /*Write lock file if requested: */ 92 92 this->parameters->FindParam(&waitonlock,SettingsWaitonlockEnum); … … 214 214 /*FUNCTION FemModel::OutputResults {{{*/ 215 215 void FemModel::OutputResults(void){ 216 216 217 217 _pprintLine_("write results to disk:"); 218 218 … … 253 253 int solution_type; 254 254 void (*solutioncore)(FemModel*)=NULL; //core solution function pointer 255 255 256 256 _pprintLine_("call computational core:"); 257 257 -
issm/trunk-jpl/src/c/classes/Hook.cpp
r13413 r13622 126 126 output->offsets = xNew<int>(output->num); 127 127 } 128 128 129 129 for(int i=0;i<output->num;i++){ 130 130 output->objects[i] = this->objects[i]; … … 182 182 /*FUNCTION Hook::delivers{{{*/ 183 183 Object* Hook::delivers(void){ 184 184 185 185 /*first, check that we only have one T object in our object list: */ 186 186 if (this->num!=1) _error_("trying to delivery a single hook object when hook holds " << this->num << " objects" << "\n"); -
issm/trunk-jpl/src/c/classes/IoModel.cpp
r13612 r13622 28 28 this->independents=NULL; 29 29 this->constants=NULL; 30 30 31 31 this->my_elements=NULL; 32 32 this->my_nodes=NULL; … … 34 34 this->singlenodetoelementconnectivity=NULL; 35 35 this->numbernodetoelementconnectivity=NULL; 36 36 37 37 this->nodecounter=0; 38 38 this->loadcounter=0; … … 68 68 this->singlenodetoelementconnectivity=NULL; 69 69 this->numbernodetoelementconnectivity=NULL; 70 70 71 71 this->nodecounter=0; 72 72 this->loadcounter=0; … … 205 205 void IoModel::DeclareIndependents(void){ 206 206 207 208 207 int i; 209 208 bool autodiff = false; … … 217 216 bool keep=false; 218 217 219 220 218 /*Initialize array detecting whether data[i] is an independent AD mode variable: */ 221 219 this->independents=xNew<bool>(MaximumNumberOfEnums); … … 226 224 227 225 this->FetchData(&numberofvertices,MeshNumberofverticesEnum); 228 226 229 227 #ifdef _HAVE_ADOLC_ 230 228 … … 266 264 else this->independent_objects=NULL; 267 265 268 269 266 } 270 267 /*}}}*/ … … 283 280 dataenum=va_arg(ap, int); 284 281 _assert_(dataenum<MaximumNumberOfEnums); 285 282 286 283 /*do not erase independent variables for the AD mode computations!: */ 287 284 if (!this->independents[dataenum]) xDelete<IssmDouble>(this->data[dataenum]); … … 300 297 301 298 int my_rank; 302 299 303 300 /*record descriptions; */ 304 301 int record_enum; … … 323 320 /*Go find in the binary file, the position of the data we want to fetch: */ 324 321 if(my_rank==0){ //cpu 0{{{ 325 322 326 323 /*First set FILE* position to the beginning of the file: */ 327 324 fseek(this->fid,0,SEEK_SET); … … 339 336 } 340 337 else{ 341 338 342 339 /* Read the record length and the data type code: */ 343 340 fread(&record_length,sizeof(int),1,this->fid); 344 341 fread(&record_code,sizeof(int),1,this->fid); 345 342 346 343 #ifdef _HAVE_MPI_ 347 344 /*Tell other cpus what we are doing: */ … … 352 349 MPI_Bcast(&record_length,1,MPI_INT,0,IssmComm::GetComm()); 353 350 #endif 354 351 355 352 switch(record_code){ 356 353 case 1: … … 416 413 string[0]='\0'; 417 414 } 418 415 419 416 /*Add string to parameters: */ 420 417 this->constants->AddObject(new StringParam(record_enum,string)); … … 478 475 /*boolean. get it from cpu 0 */ 479 476 MPI_Bcast(&booleanint,1,MPI_INT,0,IssmComm::GetComm()); 480 477 481 478 /*create BoolParam: */ 482 479 this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to a boolean … … 486 483 /*integer. get it from cpu 0 */ 487 484 MPI_Bcast(&integer,1,MPI_INT,0,IssmComm::GetComm()); 488 485 489 486 /*create IntParam: */ 490 487 this->constants->AddObject(new IntParam(record_enum,integer)); … … 494 491 /*scalar. get it from cpu 0 */ 495 492 MPI_Bcast(&scalar,1,MPI_DOUBLE,0,IssmComm::GetComm()); 496 493 497 494 /*create DoubleParam: */ 498 495 this->constants->AddObject(new DoubleParam(record_enum,scalar)); … … 530 527 } 531 528 532 533 529 } 534 530 } … … 541 537 542 538 int my_rank; 543 544 539 545 540 /*output: */ 546 541 int booleanint; 547 542 int code; 548 543 549 544 /*recover my_rank:*/ 550 545 my_rank=IssmComm::GetRank(); 551 546 552 553 547 /*Set file pointer to beginning of the data: */ 554 548 fid=this->SetFilePointerToData(&code,NULL,data_enum); 555 549 556 550 if(code!=1)_error_("expecting a boolean for enum " << EnumToStringx(data_enum)); 557 551 558 552 /*We have to read a boolean from disk. */ 559 553 if(my_rank==0){ … … 578 572 int integer; 579 573 int code; 580 574 581 575 /*recover my_rank:*/ 582 576 my_rank=IssmComm::GetRank(); 583 577 584 578 /*Set file pointer to beginning of the data: */ 585 579 fid=this->SetFilePointerToData(&code,NULL,data_enum); 586 580 587 581 if(code!=2)_error_("expecting an integer for enum " << EnumToStringx(data_enum)); 588 582 589 583 /*We have to read a integer from disk. First read the dimensions of the integer, then the integer: */ 590 584 if(my_rank==0){ … … 603 597 void IoModel::FetchData(IssmDouble* pscalar,int data_enum){ 604 598 605 606 599 int my_rank; 607 608 600 609 601 /*output: */ 610 602 IssmPDouble scalar; 611 603 int code; 612 604 613 605 /*recover my_rank:*/ 614 606 my_rank=IssmComm::GetRank(); … … 616 608 /*Set file pointer to beginning of the data: */ 617 609 fid=this->SetFilePointerToData(&code,NULL,data_enum); 618 610 619 611 if(code!=3)_error_("expecting a IssmDouble for enum " << EnumToStringx(data_enum)); 620 612 621 613 /*We have to read a scalar from disk. First read the dimensions of the scalar, then the scalar: */ 622 614 if(my_rank==0){ … … 629 621 /*Assign output pointers: */ 630 622 *pscalar=scalar; 631 623 632 624 } 633 625 /*}}}*/ … … 636 628 637 629 int my_rank; 638 639 630 640 631 /*output: */ … … 642 633 int string_size; 643 634 int code=0; 644 635 645 636 /*recover my_rank:*/ 646 637 my_rank=IssmComm::GetRank(); … … 648 639 /*Set file pointer to beginning of the data: */ 649 640 fid=this->SetFilePointerToData(&code,NULL,data_enum); 650 641 651 642 if(code!=4)_error_("expecting a string for enum " << EnumToStringx(data_enum)); 652 643 653 644 /*Now fetch: */ 654 645 655 646 /*We have to read a string from disk. First read the dimensions of the string, then the string: */ 656 647 if(my_rank==0){ … … 679 670 string[0]='\0'; 680 671 } 681 682 672 683 673 /*Assign output pointers: */ … … 697 687 int code=0; 698 688 int vector_type=0; 699 689 700 690 /*recover my_rank:*/ 701 691 my_rank=IssmComm::GetRank(); … … 705 695 706 696 if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum)); 707 697 708 698 /*Now fetch: */ 709 699 … … 733 723 if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix "); 734 724 } 735 725 736 726 #ifdef _HAVE_MPI_ 737 727 MPI_Bcast(matrix,M*N,MPI_DOUBLE,0,IssmComm::GetComm()); … … 771 761 int code=0; 772 762 int vector_type=0; 773 763 774 764 /*recover my_rank:*/ 775 765 my_rank=IssmComm::GetRank(); … … 778 768 fid=this->SetFilePointerToData(&code,&vector_type,data_enum); 779 769 if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum)); 780 770 781 771 /*Now fetch: */ 782 772 … … 832 822 833 823 int my_rank; 834 824 835 825 int i; 836 826 … … 838 828 int numstrings=0; 839 829 char** strings=NULL; 840 830 841 831 /*intermediary: */ 842 832 char* string=NULL; 843 833 int string_size; 844 834 int code; 845 835 846 836 /*recover my_rank:*/ 847 837 my_rank=IssmComm::GetRank(); … … 849 839 /*Set file pointer to beginning of the data: */ 850 840 fid=this->SetFilePointerToData(&code,NULL,data_enum); 851 841 852 842 if(code!=9)_error_("expecting a string array for enum " << EnumToStringx(data_enum)); 853 843 854 844 /*We have to read a bunch of strings from disk. First read the number of strings, and allocate: */ 855 845 if(my_rank==0){ … … 867 857 /*Go through strings, and read: */ 868 858 for(i=0;i<numstrings;i++){ 869 859 870 860 if(my_rank==0){ 871 861 if(fread(&string_size,sizeof(int),1,fid)!=1) _error_("could not read length of string "); … … 917 907 IssmPDouble *matrix = NULL; 918 908 int code; 919 909 920 910 /*recover my_rank:*/ 921 911 my_rank=IssmComm::GetRank(); … … 924 914 fid=this->SetFilePointerToData(&code,NULL,data_enum); 925 915 if(code!=8)_error_("expecting a IssmDouble mat array for enum " << EnumToStringx(data_enum)); 926 916 927 917 /*Now fetch: */ 928 918 if(my_rank==0){ … … 1000 990 int code; 1001 991 char *name = NULL; 1002 992 1003 993 /*First get option name*/ 1004 994 this->FetchData(&name,index); … … 1051 1041 /*Go through the entire list of enums and fetch the corresponding data. Add it to the iomodel->data dataset. Everything 1052 1042 *we fetch is a IssmDouble* : */ 1053 1043 1054 1044 va_start(ap,num); 1055 1045 for(i=0; i<num; i++){ 1056 1046 1057 1047 dataenum=va_arg(ap, int); 1058 1048 1059 1049 if (this->independents[dataenum]){ 1060 1050 /*this data has already been checked out! Continue: */ … … 1095 1085 int nel; 1096 1086 int numberofelements; 1097 1098 1087 1099 1088 /*variables being fetched: */ … … 1282 1271 int lastindex,index; 1283 1272 int record_length; 1284 1273 1285 1274 /*recover my_rank:*/ 1286 1275 my_rank=IssmComm::GetRank(); -
issm/trunk-jpl/src/c/classes/Patch.cpp
r13612 r13622 119 119 MPI_Status status; 120 120 #endif 121 121 122 122 /*recover my_rank:*/ 123 123 my_rank=IssmComm::GetRank(); -
issm/trunk-jpl/src/c/classes/bamg/AdjacentTriangle.cpp
r12821 r13622 8 8 9 9 namespace bamg { 10 11 10 12 11 /*Constructors/Destructors*/ -
issm/trunk-jpl/src/c/classes/bamg/BamgQuadtree.cpp
r13530 r13622 428 428 } 429 429 430 431 430 if ( n0 > 0) 432 431 { -
issm/trunk-jpl/src/c/classes/bamg/EigenMetric.cpp
r12821 r13622 29 29 b=-a11-a22; 30 30 delta=b*b - 4*(a11*a22-a21*a21); 31 32 31 33 32 /*Compute norm of M to avoid round off errors*/ -
issm/trunk-jpl/src/c/classes/bamg/Geometry.cpp
r13530 r13622 173 173 } 174 174 delete [] verticeslength; 175 175 176 176 } 177 177 else{ … … 905 905 } 906 906 907 908 907 if ((*eg0)(direction0)==(GeomVertex*)vg0) 909 908 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,direction0); //vg0 = absisce -
issm/trunk-jpl/src/c/classes/bamg/ListofIntersectionTriangles.cpp
r13530 r13622 319 319 ilast=NewItem(t,ba[0],ba[1],ba[2]); } 320 320 } // outside departure 321 322 323 321 324 322 // recherche the intersection of [a,b] with Bh Mesh. … … 392 390 Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b); 393 391 394 395 392 if (detbij >= 0) { //we find the triangle contening b 396 393 dt[0]=bamg::det((*t)[1],(*t)[2],b); -
issm/trunk-jpl/src/c/classes/bamg/Mesh.cpp
r13530 r13622 200 200 VertexOnBThEdge=0; 201 201 } 202 203 202 204 203 if(nbe) … … 1245 1244 tt[2]->SetSingleVertexToTriangleConnectivity(); 1246 1245 1247 1248 1246 // swap if the point s is on a edge 1249 1247 if(izerodet>=0) { … … 1574 1572 //initialize subdomains[isd].head as 0 1575 1573 for (isd=0;isd<nbsubdomains;isd++) subdomains[isd].head =0; 1576 1574 1577 1575 k=0; 1578 1576 for (it=0;it<nbt;it++){ … … 2539 2537 delete [] HeapArete; 2540 2538 delete [] HeapTriangle; 2541 2542 2539 2543 2540 if (OutSide|| !Gh.subdomains || !Gh.nbsubdomains ) … … 2695 2692 nbsubdomains=inew;} 2696 2693 2697 2698 2694 for (it=0;it<nbt;it++) 2699 2695 if ( mark[it] ==-1 ) … … 3075 3071 if (verbose>5) _printLine_(" do nothing: costheta > 1"); 3076 3072 } 3077 3078 3073 3079 3074 long nbqq = (nbt*3)/2; … … 4273 4268 newedges[ie].adj[1]=newedges + ie +1; 4274 4269 R2 A = edges[i][0],B = edges[i][1]; 4275 4276 4270 4277 4271 kk += (i == edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1]))); … … 4343 4337 nbv = k; 4344 4338 4345 4346 4339 kedge = new long[3*nbt+1]; 4347 4340 ksplitarray = new long[nbt+1]; … … 4412 4405 const Triangle & tt = ta; 4413 4406 4414 4415 4407 const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]]; 4416 4408 const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]]; … … 4648 4640 BTh.vertices[i].m = BTh.vertices[i].m*2.; 4649 4641 4650 4651 4642 ret = 2; 4652 4643 if (nbt>= maxnbt) goto Error; // bug … … 5001 4992 _assert_(nbv<maxnbv); 5002 4993 vertices[nbv]=Gh[i]; 5003 4994 5004 4995 //Add pointer from geometry (Gh) to vertex from mesh (Th) 5005 4996 Gh[i].MeshVertexHook=vertices+nbv; … … 5196 5187 va = vb; 5197 5188 } 5198 5189 5199 5190 /*We just added one edge to the curve: Go to the next one*/ 5200 5191 lcurve = lcurveb; … … 5371 5362 /*Get curve number*/ 5372 5363 int nc=ei.GeomEdgeHook->CurveNumber; 5373 5364 5374 5365 //_printLine_("Dealing with curve number " << nc); 5375 5366 //_printLine_("edge on geometry is same as GhCurve? " << (ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no"); … … 5419 5410 /*Get edge of Bth with index iedge*/ 5420 5411 Edge &ei = BTh.edges[iedge]; 5421 5412 5422 5413 /*Initialize variables*/ 5423 5414 double Lstep=0,Lcurve=0; // step between two points (phase==1) … … 5536 5527 }// for(;;) end of the curve 5537 5528 } 5538 5539 5529 5540 5530 if (phase){ // construction of the last edge … … 5787 5777 BamgVertex & s2= (*t2)[OppositeVertex[a2]]; 5788 5778 5789 5790 5779 Icoor2 dets2 = det(*pva,*pvb,s2); 5791 5780 Icoor2 det1=t1->det , det2=t2->det ; -
issm/trunk-jpl/src/c/classes/bamg/Metric.cpp
r13036 r13622 16 16 /*FUNCTION Metric::Metric(double a){{{*/ 17 17 Metric::Metric(double a): a11(1/(a*a)),a21(0),a22(1/(a*a)){ 18 18 19 19 }/*}}}*/ 20 20 /*FUNCTION Metric::Metric(double a,double b,double c){{{*/ 21 21 Metric::Metric(double a,double b,double c) :a11(a),a21(b),a22(c){ 22 22 23 23 }/*}}}*/ 24 24 /*FUNCTION Metric::Metric(const double a[3],const Metric& m0, const Metric& m1,const Metric& m2 ){{{*/ -
issm/trunk-jpl/src/c/classes/bamg/Triangle.cpp
r13530 r13622 424 424 sinb12 = double(det2), 425 425 sinba2 = double(t2->det); 426 427 426 428 427 // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) -
issm/trunk-jpl/src/c/classes/bamg/VertexOnVertex.cpp
r12821 r13622 17 17 /*FUNCTION VertexOnVertex::VertexOnVertex(BamgVertex * w,BamgVertex *bw){{{*/ 18 18 VertexOnVertex::VertexOnVertex(BamgVertex * w,BamgVertex *bw) :v(w),bv(bw){ 19 19 20 20 }/*}}}*/ 21 21 -
issm/trunk-jpl/src/c/classes/dakota/DakotaPlugin.cpp
r13548 r13622 14 14 * the entire computations. 15 15 */ 16 17 16 18 17 #ifdef HAVE_CONFIG_H … … 77 76 memcpy(variable_descriptor,label.c_str(),(strlen(label.c_str())+1)*sizeof(char)); 78 77 79 80 78 variable_descriptors[i]=variable_descriptor; 81 79 } … … 109 107 } // namespace SIM 110 108 111 112 109 #endif //only works if dakota library has been compiled in. -
issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp
r13216 r13622 230 230 /*FUNCTION ElementMatrix::~ElementMatrix(){{{*/ 231 231 ElementMatrix::~ElementMatrix(){ 232 232 233 233 xDelete<IssmDouble>(this->values); 234 234 xDelete<int>(this->gglobaldoflist); … … 278 278 } 279 279 280 281 280 if((this->row_ssize!=0) && (this->row_fsize!=0)){ 282 281 /*first, retrieve values that are in the f and s-set from the g-set values matrix: */ -
issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp
r13216 r13622 139 139 /*fill values with 0: */ 140 140 this->values=xNewZeroInit<IssmDouble>(this->nrows); 141 141 142 142 /*g list*/ 143 143 this->gglobaldoflist=GetGlobalDofList(nodes,numnodes,GsetEnum,approximation); … … 151 151 /*FUNCTION ElementVector::~ElementVector(){{{*/ 152 152 ElementVector::~ElementVector(){ 153 153 154 154 xDelete<IssmDouble>(this->values); 155 155 xDelete<int>(this->gglobaldoflist); … … 181 181 xDelete<IssmDouble>(localvalues); 182 182 } 183 183 184 184 } 185 185 /*}}}*/ -
issm/trunk-jpl/src/c/classes/objects/Constraints/SpcDynamic.cpp
r13414 r13622 41 41 } 42 42 /*}}}*/ 43 43 44 44 /*Object virtual functions definitions:*/ 45 45 /*FUNCTION SpcDynamic::Echo {{{*/ … … 111 111 /*FUNCTION SpcDynamic::GetNodeId {{{*/ 112 112 int SpcDynamic::GetNodeId(){ 113 113 114 114 return nodeid; 115 115 } -
issm/trunk-jpl/src/c/classes/objects/Constraints/SpcStatic.cpp
r13414 r13622 40 40 } 41 41 /*}}}*/ 42 42 43 43 /*Object virtual functions definitions:*/ 44 44 /*FUNCTION SpcStatic::Echo {{{*/ … … 112 112 /*FUNCTION SpcStatic::GetNodeId {{{*/ 113 113 int SpcStatic::GetNodeId(){ 114 114 115 115 return nodeid; 116 116 } -
issm/trunk-jpl/src/c/classes/objects/Constraints/SpcTransient.cpp
r13414 r13622 54 54 } 55 55 /*}}}*/ 56 56 57 57 /*Object virtual functions definitions:*/ 58 58 /*FUNCTION SpcTransient::Echo {{{*/ … … 97 97 /*FUNCTION SpcTransient::InAnalysis{{{*/ 98 98 bool SpcTransient::InAnalysis(int in_analysis_type){ 99 99 100 100 if (in_analysis_type==this->analysis_type) return true; 101 101 else return false; … … 114 114 /*Chase through nodes and find the node to which this SpcTransient applys: */ 115 115 node=(Node*)nodes->GetObjectById(NULL,nodeid); 116 116 117 117 if(node){ //in case the spc is dealing with a node on another cpu 118 118 … … 159 159 /*FUNCTION SpcTransient::GetNodeId {{{*/ 160 160 int SpcTransient::GetNodeId(){ 161 161 162 162 return nodeid; 163 163 } … … 168 168 } 169 169 /*}}}*/ 170 -
issm/trunk-jpl/src/c/classes/objects/DependentObject.cpp
r13608 r13622 41 41 } 42 42 /*}}}*/ 43 44 43 45 44 /*Object virtual functions definitions:*/ … … 100 99 parameters->SetParam(this->index,IndexEnum); 101 100 } 102 101 103 102 ::Responsex(poutput_value,elements,nodes,vertices,loads,materials,parameters,this->name,false,0); 104 103 -
issm/trunk-jpl/src/c/classes/objects/ElementResults/BoolElementResult.cpp
r13414 r13622 110 110 /*FUNCTION BoolElementResult::PatchFill{{{*/ 111 111 void BoolElementResult::PatchFill(int row, Patch* patch){ 112 112 113 113 /*Here, we fill the result information into the patch object. First, let's remember what is in a row 114 114 * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values -
issm/trunk-jpl/src/c/classes/objects/ElementResults/DoubleElementResult.cpp
r13414 r13622 100 100 /*FUNCTION DoubleElementResult::ProcessUnits{{{*/ 101 101 void DoubleElementResult::ProcessUnits(Parameters* parameters){ 102 102 103 103 this->value=UnitConversion(this->value,IuToExtEnum,this->enum_type); 104 104 … … 112 112 /*FUNCTION DoubleElementResult::PatchFill{{{*/ 113 113 void DoubleElementResult::PatchFill(int row, Patch* patch){ 114 114 115 115 /*Here, we fill the result information into the patch object. First, let's remember what is in a row 116 116 * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values -
issm/trunk-jpl/src/c/classes/objects/ElementResults/PentaP1ElementResult.cpp
r13414 r13622 70 70 /*FUNCTION PentaP1ElementResult::copy{{{*/ 71 71 Object* PentaP1ElementResult::copy() { 72 72 73 73 return new PentaP1ElementResult(this->enum_type,this->values,this->step,this->time); 74 74 … … 111 111 /*FUNCTION PentaP1ElementResult::ProcessUnits{{{*/ 112 112 void PentaP1ElementResult::ProcessUnits(Parameters* parameters){ 113 113 114 114 UnitConversion(this->values,6,IuToExtEnum,this->enum_type); 115 115 … … 123 123 /*FUNCTION PentaP1ElementResult::PatchFill{{{*/ 124 124 void PentaP1ElementResult::PatchFill(int row, Patch* patch){ 125 125 126 126 /*Here, we fill the result information into the patch object. First, let's remember what is in a row 127 127 * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values -
issm/trunk-jpl/src/c/classes/objects/ElementResults/TriaP1ElementResult.cpp
r13414 r13622 49 49 /*FUNCTION TriaP1ElementResult::DeepEcho{{{*/ 50 50 void TriaP1ElementResult::DeepEcho(void){ 51 51 52 52 _printLine_("TriaP1ElementResult:"); 53 53 _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")"); … … 69 69 /*FUNCTION TriaP1ElementResult::copy{{{*/ 70 70 Object* TriaP1ElementResult::copy() { 71 71 72 72 return new TriaP1ElementResult(this->enum_type,this->values,this->step,this->time); 73 73 … … 99 99 /*FUNCTION TriaP1ElementResult::ProcessUnits{{{*/ 100 100 void TriaP1ElementResult::ProcessUnits(Parameters* parameters){ 101 101 102 102 UnitConversion(this->values,3,IuToExtEnum,this->enum_type); 103 103 … … 111 111 /*FUNCTION TriaP1ElementResult::PatchFill{{{*/ 112 112 void TriaP1ElementResult::PatchFill(int row, Patch* patch){ 113 113 114 114 /*Here, we fill the result information into the patch object. First, let's remember what is in a row 115 115 * of the patch object: enum_type step time element_id interpolation vertices_ids nodal_values -
issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp
r13528 r13622 82 82 this->inputs=new Inputs(); 83 83 this->results=new Results(); 84 84 85 85 /*initialize pointers:*/ 86 86 this->nodes=NULL; … … 186 186 GaussPenta* gauss=NULL; 187 187 188 189 188 /* Basal friction can only be found at the base of an ice sheet: */ 190 189 if (!IsOnBed() || IsFloating()){ … … 199 198 Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); 200 199 201 202 200 /*Build friction element, needed later: */ 203 201 friction=new Friction("3d",inputs,matpar,DiagnosticHorizAnalysisEnum); … … 216 214 count++; 217 215 } 218 216 219 217 /*Create PentaVertex input, which will hold the basal friction:*/ 220 218 this->inputs->AddInput(new PentaP1Input(BasalFrictionEnum,&basalfriction[0])); … … 257 255 /*retrieve some parameters: */ 258 256 this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum); 259 257 260 258 if(!IsOnBed()){ 261 259 //put zero … … 364 362 sigma_yz[iv]=2*viscosity*epsilon[5]; 365 363 } 366 364 367 365 /*Add Stress tensor components into inputs*/ 368 366 this->inputs->AddInput(new PentaP1Input(StressTensorxxEnum,&sigma_xx[0])); … … 381 379 382 380 int analysis_counter; 383 381 384 382 /*go into parameters and get the analysis_counter: */ 385 383 parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); … … 421 419 _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); 422 420 /*}}}*/ 423 421 424 422 /*Skip if water element*/ 425 423 if(IsOnWater()) return; … … 652 650 653 651 int i; 654 652 655 653 _printLine_("Penta:"); 656 654 _printLine_(" id: " << id); … … 1191 1189 /*FUNCTION Penta::Sid {{{*/ 1192 1190 int Penta::Sid(){ 1193 1191 1194 1192 return sid; 1195 1193 … … 1251 1249 /*Check that name is an element input*/ 1252 1250 if (!IsInput(name)) return; 1253 1251 1254 1252 if ((code==5) || (code==1)){ //boolean 1255 1253 this->inputs->AddInput(new BoolInput(name,reCast<bool,IssmDouble>(scalar))); … … 1873 1871 penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id); 1874 1872 } 1875 1873 1876 1874 /*Free ressources:*/ 1877 1875 xDelete<int>(doflist); … … 1897 1895 /*Add input to the element: */ 1898 1896 this->inputs->AddInput(new PentaP1Input(enum_type,values)); 1899 1897 1900 1898 /*Free ressources:*/ 1901 1899 xDelete<int>(doflist); … … 1937 1935 penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id); 1938 1936 } 1939 1937 1940 1938 /*Free ressources:*/ 1941 1939 xDelete<int>(doflist); … … 2170 2168 rho_ice=matpar->GetRhoIce(); 2171 2169 density=rho_ice/rho_water; 2172 2170 2173 2171 /*go through vertices, and update inputs, considering them to be PentaVertex type: */ 2174 2172 for(i=0;i<NUMVERTICES;i++){ … … 2211 2209 } 2212 2210 } 2213 2211 2214 2212 /*Add basal melting rate if element just ungrounded*/ 2215 2213 if(!this->IsFloating() && elementonshelf==true){ … … 2293 2291 /*recover pointer: */ 2294 2292 count=*pcount; 2295 2293 2296 2294 /*will be needed later: */ 2297 2295 for(i=0;i<6;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch. … … 2502 2500 /*FUNCTION Penta::RequestedOutput{{{*/ 2503 2501 void Penta::RequestedOutput(int output_enum,int step,IssmDouble time){ 2504 2502 2505 2503 if(IsInput(output_enum)){ 2506 2504 /*just transfer this input to results, and we are done: */ … … 2677 2675 GetInputListOnVertices(&a_neg[0],SurfaceforcingsANegEnum); 2678 2676 GetInputListOnVertices(&b_neg[0],SurfaceforcingsBNegEnum); 2679 2677 2680 2678 /*Recover surface elevatio at vertices: */ 2681 2679 GetInputListOnVertices(&h[0],ThicknessEnum); … … 2685 2683 rho_ice=matpar->GetRhoIce(); 2686 2684 rho_water=matpar->GetRhoFreshwater(); 2687 2685 2688 2686 // loop over all vertices 2689 2687 for(i=0;i<NUMVERTICES;i++){ … … 2791 2789 minz=xyz_list[0][2]; 2792 2790 maxz=xyz_list[0][2]; 2793 2791 2794 2792 for(i=1;i<NUMVERTICES;i++){ 2795 2793 if (xyz_list[i][0]<minx)minx=xyz_list[i][0]; … … 2973 2971 if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[nodes[i]->Sid()])){ 2974 2972 vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL); 2975 2973 2976 2974 /*If node was not on ice shelf, we flipped*/ 2977 2975 if(nodes_on_iceshelf[nodes[i]->Sid()]==0){ … … 3013 3011 for (int iv=0;iv<NUMVERTICES;iv++){ 3014 3012 gauss->GaussVertex(iv); 3015 3013 3016 3014 thickness_input->GetInputValue(&thickness,gauss); 3017 3015 … … 3019 3017 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 3020 3018 GetPhi(&phi, &epsilon[0], viscosity); 3021 3022 3019 3023 3020 viscousheating[iv]=phi*thickness; … … 3274 3271 3275 3272 if(IsOnWater() || !IsOnSurface()) return 0.; 3276 3273 3277 3274 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 3278 3275 … … 3286 3283 smb_input->GetInputAverage(&smb); 3287 3284 Total_Smb=rho_ice*base*smb;// smb on element in kg s-1 3288 3285 3289 3286 /*Process units: */ 3290 3287 Total_Smb=UnitConversion(Total_Smb,IuToExtEnum,TotalSmbEnum);// smb on element in GigaTon yr-1 3291 3288 3292 3289 /*Return: */ 3293 3290 return Total_Smb; … … 3299 3296 /*FUNCTION Penta::CreateKMatrixEnthalpy {{{*/ 3300 3297 ElementMatrix* Penta::CreateKMatrixEnthalpy(void){ 3301 3298 3302 3299 /*compute all stiffness matrices for this element*/ 3303 3300 ElementMatrix* Ke1=CreateKMatrixEnthalpyVolume(); 3304 3301 ElementMatrix* Ke2=CreateKMatrixEnthalpyShelf(); 3305 3302 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); 3306 3303 3307 3304 /*clean-up and return*/ 3308 3305 delete Ke1; … … 3504 3501 3505 3502 gauss->GaussPoint(ig); 3506 3503 3507 3504 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); 3508 3505 GetNodalFunctionsP1(&basis[0], gauss); 3509 3506 3510 3507 D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(rho_ice*heatcapacity); 3511 3508 if(reCast<bool,IssmDouble>(dt)) D_scalar=dt*D_scalar; … … 3516 3513 &Ke->values[0],1); 3517 3514 } 3518 3515 3519 3516 /*Clean up and return*/ 3520 3517 delete gauss; … … 3536 3533 /*FUNCTION Penta::CreateKMatrixThermal {{{*/ 3537 3534 ElementMatrix* Penta::CreateKMatrixThermal(void){ 3538 3535 3539 3536 /*compute all stiffness matrices for this element*/ 3540 3537 ElementMatrix* Ke1=CreateKMatrixThermalVolume(); 3541 3538 ElementMatrix* Ke2=CreateKMatrixThermalShelf(); 3542 3539 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); 3543 3540 3544 3541 /*clean-up and return*/ 3545 3542 delete Ke1; … … 3703 3700 ElementMatrix* Penta::CreateKMatrixThermalShelf(void){ 3704 3701 3705 3706 3702 /*Constants*/ 3707 3703 const int numdof=NDOF1*NUMVERTICES; … … 3737 3733 3738 3734 gauss->GaussPoint(ig); 3739 3735 3740 3736 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); 3741 3737 GetNodalFunctionsP1(&basis[0], gauss); 3742 3738 3743 3739 D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice); 3744 3740 if(reCast<bool,IssmDouble>(dt)) D_scalar=dt*D_scalar; … … 3749 3745 &Ke->values[0],1); 3750 3746 } 3751 3747 3752 3748 /*Clean up and return*/ 3753 3749 delete gauss; … … 4409 4405 GetInputListOnVertices(&pressure[0],PressureEnum); 4410 4406 Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); 4411 4407 4412 4408 this->inputs->GetInputValue(&converged,ConvergedEnum); 4413 4409 if(converged){ … … 4418 4414 //if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector"); 4419 4415 } 4420 4416 4421 4417 this->inputs->AddInput(new PentaP1Input(EnthalpyEnum,values)); 4422 4418 this->inputs->AddInput(new PentaP1Input(WaterfractionEnum,waterfraction)); … … 4473 4469 input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); 4474 4470 } 4475 4471 4476 4472 else{ 4477 4473 input=inputs->GetInput(enum_type); … … 5533 5529 new_input = new PentaP1Input(control_enum,values); 5534 5530 5535 5536 5531 if(control_enum==MaterialsRheologyBbarEnum){ 5537 5532 input=(Input*)material->inputs->GetInput(control_enum); _assert_(input); … … 5553 5548 /*FUNCTION Penta::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);{{{*/ 5554 5549 void Penta::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){ 5555 5550 5556 5551 int i,j; 5557 5552 … … 5580 5575 IssmDouble surface[6]; 5581 5576 IssmDouble bed[6]; 5582 5577 5583 5578 /*retrieve inputs: */ 5584 5579 GetInputListOnVertices(&thickness_init[0],ThicknessEnum); … … 5665 5660 /*FUNCTION Penta::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type);{{{*/ 5666 5661 void Penta::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){ 5667 5662 5668 5663 int i,j,t; 5669 5664 TransientInput* transientinput=NULL; … … 5679 5674 5680 5675 case VertexEnum: 5681 5676 5682 5677 /*Create transient input: */ 5683 5678 5684 5679 parameters->FindParam(&yts,ConstantsYtsEnum); 5685 5680 … … 5752 5747 /*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattyn{{{*/ 5753 5748 ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattyn(void){ 5754 5749 5755 5750 /*compute all stiffness matrices for this element*/ 5756 5751 ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealPattynViscous(); 5757 5752 ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealPattynFriction(); 5758 5753 ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); 5759 5754 5760 5755 /*clean-up and return*/ 5761 5756 delete Ke1; … … 5864 5859 const int numdof = NDOF2 *NUMVERTICES; 5865 5860 const int numdoftotal = NDOF4 *NUMVERTICES; 5866 5861 5867 5862 /*Intermediaries */ 5868 5863 int i,j,ig,analysis_type; … … 5933 5928 DL_scalar=alpha2*gauss->weight*Jdet2d; 5934 5929 for (i=0;i<2;i++) DL[i][i]=DL_scalar; 5935 5930 5936 5931 /* Do the triple producte tL*D*L: */ 5937 5932 TripleMultiply( &L[0][0],2,numdof,1, … … 6165 6160 DLStokesMacAyeal[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2]; 6166 6161 DLStokesMacAyeal[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2]; 6167 6162 6168 6163 TripleMultiply( &LMacAyealStokes[0][0],8,numdof2dm,1, 6169 6164 &DLMacAyealStokes[0][0],8,8,0, … … 6700 6695 /*Constants*/ 6701 6696 const int numdof = NDOF2*NUMVERTICES; 6702 6697 6703 6698 /*Intermediaries */ 6704 6699 int i,j,ig; … … 6751 6746 alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT); 6752 6747 } 6753 6748 6754 6749 DL_scalar=alpha2*gauss->weight*Jdet; 6755 6750 for (i=0;i<2;i++) DL[i][i]=DL_scalar; 6756 6751 6757 6752 TripleMultiply( &L[0][0],2,numdof,1, 6758 6753 &DL[0][0],2,2,0, … … 6926 6921 /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/ 6927 6922 //TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); 6928 6923 6929 6924 /*Clean up and return*/ 6930 6925 delete gauss; … … 6935 6930 /*FUNCTION Penta::CreateKMatrixDiagnosticVert {{{*/ 6936 6931 ElementMatrix* Penta::CreateKMatrixDiagnosticVert(void){ 6937 6932 6938 6933 /*compute all stiffness matrices for this element*/ 6939 6934 ElementMatrix* Ke1=CreateKMatrixDiagnosticVertVolume(); … … 7816 7811 ElementVector* Penta::CreatePVectorDiagnosticVertBase(void){ 7817 7812 7818 7819 7813 /*Constants*/ 7820 7814 const int numdof=NDOF1*NUMVERTICES; … … 8415 8409 penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id); 8416 8410 } 8417 8411 8418 8412 /*Free ressources:*/ 8419 8413 xDelete<int>(doflist); … … 8687 8681 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattyn {{{*/ 8688 8682 void Penta::InputUpdateFromSolutionDiagnosticPattyn(IssmDouble* solution){ 8689 8683 8690 8684 const int numdof=NDOF2*NUMVERTICES; 8691 8685 … … 8855 8849 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHutter {{{*/ 8856 8850 void Penta::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){ 8857 8851 8858 8852 const int numdof=NDOF2*NUMVERTICES; 8859 8853 … … 8920 8914 8921 8915 const int numdof=NDOF1*NUMVERTICES; 8922 8916 8923 8917 int i; 8924 8918 int approximation; … … 8937 8931 int* doflist = NULL; 8938 8932 8939 8940 8933 /*Get the approximation and do nothing if the element in Stokes or None*/ 8941 8934 inputs->GetInputValue(&approximation,ApproximationEnum); … … 9022 9015 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticStokes {{{*/ 9023 9016 void Penta::InputUpdateFromSolutionDiagnosticStokes(IssmDouble* solution){ 9024 9017 9025 9018 const int numdof=NDOF4*NUMVERTICES; 9026 9019 … … 9062 9055 for(i=0;i<NUMVERTICES;i++) pressure[i]=pressure[i]*stokesreconditioning; 9063 9056 for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5); 9064 9057 9065 9058 /*Now, we have to move the previous inputs to old 9066 9059 * status, otherwise, we'll wipe them off: */ … … 9132 9125 /*}}}*/ 9133 9126 #endif 9134 -
issm/trunk-jpl/src/c/classes/objects/Elements/PentaHook.cpp
r13129 r13622 49 49 /*intermediary: */ 50 50 int matpar_id; 51 51 52 52 /*retrieve parameters: */ 53 53 iomodel->Constant(&matpar_id,MeshNumberofelementsEnum); matpar_id++; -
issm/trunk-jpl/src/c/classes/objects/Elements/PentaRef.cpp
r13036 r13622 699 699 IssmDouble l1l2l3[NUMNODESP1_2d]; 700 700 701 702 701 /*Get l1l2l3 in actual coordinate system: */ 703 702 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; … … 812 811 IssmDouble l1l2l3[NUMNODESP1_2d]; 813 812 814 815 813 /*Get l1l2l3 in actual coordinate system: */ 816 814 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; -
issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp
r13521 r13622 42 42 :TriaRef(nummodels) 43 43 ,TriaHook(nummodels,index+1,iomodel){ 44 44 45 45 int i; 46 46 /*id: */ … … 175 175 _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); 176 176 /*}}}*/ 177 177 178 178 /*Skip if water element*/ 179 179 if(IsOnWater()) return; … … 499 499 500 500 gauss->GaussPoint(ig); 501 501 502 502 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 503 503 D=gauss->weight*Jdet; … … 528 528 _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); 529 529 /*}}}*/ 530 530 531 531 /*Skip if water element*/ 532 532 if(IsOnWater()) return; … … 691 691 /*Constants*/ 692 692 const int numdof=NDOF1*NUMVERTICES; 693 693 694 694 /*Intermediaries */ 695 695 int i,j,ig; … … 714 714 slope_input=inputs->GetInput(BedEnum); _assert_(slope_input); 715 715 } 716 716 717 717 /* Start looping on the number of gaussian points: */ 718 718 gauss=new GaussTria(2); … … 820 820 sigma_xy[iv]=2*viscosity*epsilon[2]; 821 821 } 822 822 823 823 /*Add Stress tensor components into inputs*/ 824 824 this->inputs->AddInput(new TriaP1Input(StressTensorxxEnum,&sigma_xx[0])); … … 835 835 /*FUNCTION Tria::Configure {{{*/ 836 836 void Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ 837 837 838 838 /*go into parameters and get the analysis_counter: */ 839 839 int analysis_counter; … … 894 894 _printLine_("neighboor sids: "); 895 895 _printLine_(" " << horizontalneighborsids[0] << " " << horizontalneighborsids[1] << " " << horizontalneighborsids[2]); 896 896 897 897 return; 898 898 } … … 946 946 this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,(finaltime-(21000*yts))); 947 947 this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,time); 948 948 949 949 /*Compute the temperature and precipitation*/ 950 950 for(int iv=0;iv<NUMVERTICES;iv++){ … … 1027 1027 x2=xyz_list[1][0]; y2=xyz_list[1][1]; 1028 1028 x3=xyz_list[2][0]; y3=xyz_list[2][1]; 1029 1029 1030 1030 _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0); 1031 1031 return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2; … … 1284 1284 /*FUNCTION Tria::Id {{{*/ 1285 1285 int Tria::Id(){ 1286 1286 1287 1287 return id; 1288 1288 … … 1291 1291 /*FUNCTION Tria::Sid {{{*/ 1292 1292 int Tria::Sid(){ 1293 1293 1294 1294 return sid; 1295 1295 … … 1404 1404 * object out of the input, with the additional step and time information: */ 1405 1405 this->results->AddObject((Object*)input->SpawnResult(step,time)); 1406 1406 1407 1407 #ifdef _HAVE_CONTROL_ 1408 1408 if(input->ObjectEnum()==ControlInputEnum){ … … 1452 1452 IssmDouble yts; 1453 1453 int num_cm_responses; 1454 1454 1455 1455 /*Get parameters: */ 1456 1456 iomodel->Constant(&yts,ConstantsYtsEnum); … … 1719 1719 /*Check that name is an element input*/ 1720 1720 if (!IsInput(name)) return; 1721 1721 1722 1722 if ((code==5) || (code==1)){ //boolean 1723 1723 this->inputs->AddInput(new BoolInput(name,reCast<bool>(scalar))); … … 1746 1746 int numberofelements; 1747 1747 IssmDouble yts; 1748 1749 1748 1750 1749 /*Fetch parameters: */ … … 1861 1860 /*FUNCTION Tria::IsOnBed {{{*/ 1862 1861 bool Tria::IsOnBed(){ 1863 1862 1864 1863 bool onbed; 1865 1864 inputs->GetInputValue(&onbed,MeshElementonbedEnum); … … 1984 1983 rho_ice=matpar->GetRhoIce(); 1985 1984 density=rho_ice/rho_water; 1986 1985 1987 1986 /*go through vertices, and update inputs, considering them to be TriaVertex type: */ 1988 1987 for(i=0;i<NUMVERTICES;i++){ … … 2025 2024 } 2026 2025 } 2027 2026 2028 2027 /*Add basal melting rate if element just ungrounded*/ 2029 2028 if(!this->IsFloating() && elementonshelf==true){ … … 2081 2080 /*recover pointer: */ 2082 2081 row=*prow; 2083 2082 2084 2083 for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch. 2085 2084 … … 2283 2282 GetInputListOnVertices(&a_neg[0],SurfaceforcingsANegEnum); 2284 2283 GetInputListOnVertices(&b_neg[0],SurfaceforcingsBNegEnum); 2285 2284 2286 2285 /*Recover surface elevatio at vertices: */ 2287 2286 GetInputListOnVertices(&h[0],ThicknessEnum); … … 2291 2290 rho_ice=matpar->GetRhoIce(); 2292 2291 rho_water=matpar->GetRhoFreshwater(); 2293 2292 2294 2293 // loop over all vertices 2295 2294 for(i=0;i<NUMVERTICES;i++){ … … 2319 2318 /*FUNCTION Tria::SetCurrentConfiguration {{{*/ 2320 2319 void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ 2321 2320 2322 2321 /*go into parameters and get the analysis_counter: */ 2323 2322 int analysis_counter; … … 2540 2539 if (reCast<bool>(vertices_potentially_ungrounding[nodes[i]->Sid()])){ 2541 2540 vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL); 2542 2541 2543 2542 /*If node was not on ice shelf, we flipped*/ 2544 2543 if(nodes_on_iceshelf[nodes[i]->Sid()]==0){ … … 2849 2848 smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1 2850 2849 Total_Smb=rho_ice*base*smb; // smb on element in kg s-1 2851 2850 2852 2851 /*Process units: */ 2853 2852 Total_Smb=UnitConversion(Total_Smb,IuToExtEnum,TotalSmbEnum); // smb on element in GigaTon yr-1 2854 2853 2855 2854 /*Return: */ 2856 2855 return Total_Smb; … … 2867 2866 ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyealFriction(); 2868 2867 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); 2869 2868 2870 2869 /*clean-up and return*/ 2871 2870 delete Ke1; … … 2991 2990 DL_scalar=alpha2*gauss->weight*Jdet; 2992 2991 for (i=0;i<2;i++) DL[i][i]=DL_scalar; 2993 2992 2994 2993 TripleMultiply( &L[0][0],2,numdof,1, 2995 2994 &DL[0][0],2,2,0, … … 3278 3277 /*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHoriz {{{*/ 3279 3278 void Tria::InputUpdateFromSolutionDiagnosticHoriz(IssmDouble* solution){ 3280 3279 3281 3280 const int numdof=NDOF2*NUMVERTICES; 3282 3281 … … 3291 3290 IssmDouble pressure[NUMVERTICES]; 3292 3291 IssmDouble thickness[NUMVERTICES]; 3293 3292 3294 3293 /*Get dof list: */ 3295 3294 GetDofList(&doflist,NoneApproximationEnum,GsetEnum); … … 3341 3340 /*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHutter {{{*/ 3342 3341 void Tria::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){ 3343 3342 3344 3343 const int numdof=NDOF2*NUMVERTICES; 3345 3344 3346 3345 int i; 3347 3346 int* doflist=NULL; … … 3354 3353 IssmDouble pressure[NUMVERTICES]; 3355 3354 IssmDouble thickness[NUMVERTICES]; 3356 3355 3357 3356 /*Get dof list: */ 3358 3357 GetDofList(&doflist,NoneApproximationEnum,GsetEnum); … … 3800 3799 /*Build alpha_complement_list: */ 3801 3800 friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum); 3802 3801 3803 3802 dragcoefficient_input->GetInputValue(&drag, gauss); 3804 3803 adjointx_input->GetInputValue(&lambda, gauss); … … 3812 3811 grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i]; 3813 3812 } 3814 3813 3815 3814 /*Add gradje_g_gaussian vector to gradje_g: */ 3816 3815 for(i=0;i<NUMVERTICES;i++){ … … 3930 3929 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 3931 3930 GetNodalFunctions(basis, gauss); 3932 3931 3933 3932 adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); 3934 3933 thickness_input->GetInputValue(&thickness, gauss); … … 4197 4196 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input); 4198 4197 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 4199 4198 4200 4199 /* Start looping on the number of gaussian points: */ 4201 4200 gauss=new GaussTria(4); … … 5303 5302 /* compute VelocityFactor */ 5304 5303 VelocityFactor= n_man*CR*CR*rho_water*g/mu_water; 5305 5304 5306 5305 gauss=new GaussTria(); 5307 5306 for (int iv=0;iv<NUMVERTICES;iv++){ … … 5476 5475 else for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i]; 5477 5476 } 5478 5477 5479 5478 /*Clean up and return*/ 5480 5479 delete gauss; … … 5548 5547 /*FUNCTION Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);{{{*/ 5549 5548 void Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){ 5550 5549 5551 5550 int i,j; 5552 5551 … … 5574 5573 IssmDouble surface[3]; 5575 5574 IssmDouble bed[3]; 5576 5575 5577 5576 /*retrieve inputs: */ 5578 5577 GetInputListOnVertices(&thickness_init[0],ThicknessEnum); … … 5647 5646 /*FUNCTION Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type);{{{*/ 5648 5647 void Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){ 5649 5648 5650 5649 int i,j,t; 5651 5650 TransientInput* transientinput=NULL; … … 5661 5660 5662 5661 case VertexEnum: 5663 5662 5664 5663 /*Create transient input: */ 5665 5664 5666 5665 parameters->FindParam(&yts,ConstantsYtsEnum); 5667 5666 for(t=0;t<ncols;t++){ //ncols is the number of times … … 5889 5888 /*Constants*/ 5890 5889 const int numdof=NDOF1*NUMVERTICES; 5891 5890 5892 5891 /*Intermediaries */ 5893 5892 int i,j,ig; … … 5905 5904 Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); 5906 5905 Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); 5907 5906 5908 5907 /* Start looping on the number of gaussian points: */ 5909 5908 gauss=new GaussTria(2); -
issm/trunk-jpl/src/c/classes/objects/Elements/TriaHook.cpp
r13129 r13622 50 50 /*retrieve parameters: */ 51 51 iomodel->Constant(&matpar_id,MeshNumberofelementsEnum); matpar_id++; 52 52 53 53 this->numanalyses=in_numanalyses; 54 54 this->hnodes= new Hook*[in_numanalyses]; -
issm/trunk-jpl/src/c/classes/objects/Elements/TriaRef.cpp
r13036 r13622 319 319 y3=*(xyz_list+NUMNODES*2+1); 320 320 321 322 321 *(J+NDOF2*0+0)=0.5*(x2-x1); 323 322 *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2); -
issm/trunk-jpl/src/c/classes/objects/IndependentObject.cpp
r13612 r13622 77 77 } /*}}}*/ 78 78 79 80 79 /*IndependentObject methods: */ 81 80 /*FUNCTION IndependentObject::FetchIndependent{{{*/ … … 84 83 int my_rank; 85 84 FILE* fid=NULL; 86 85 87 86 /*recover my_rank:*/ 88 87 my_rank=IssmComm::GetRank(); 89 88 90 89 #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents: 91 90 92 91 /*Branch according to the type of variable: */ 93 92 if(type==0){ /*scalar: {{{*/ … … 111 110 /*Now, before we even broadcast this to other nodes, declare the scalar as an independent variable!: */ 112 111 scalar<<=pscalar; 113 112 114 113 #ifdef _HAVE_MPI_ 115 114 MPI_Bcast(&scalar,1,MPI_DOUBLE,0,IssmComm::GetComm()); 116 115 #endif 117 118 116 119 117 /*Ok, we are almost done. scalar is now an independent variable. We don't want this variable to be fetched again in the … … 122 120 space for it: */ 123 121 scalar_slot=xNew<IssmDouble>(1); *scalar_slot=scalar; 124 122 125 123 iomodel->data[name]=scalar_slot; 126 124 iomodel->independents[name]=true; 127 125 128 129 126 } /*}}}*/ 130 127 else if(type==1){ /* vector: {{{*/ 131 128 132 129 FILE* fid=NULL; 133 130 int M,N; … … 136 133 int code=0; 137 134 int vector_type=0; 138 135 139 136 /*Set file pointer to beginning of the data: */ 140 137 fid=iomodel->SetFilePointerToData(&code,&vector_type,name); 141 138 if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(name)); 142 139 143 140 /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ 144 141 /*numberofelements: */ … … 165 162 if(my_rank==0){ 166 163 if(fread(buffer,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix "); 167 164 168 165 /*Now, before we even broadcast this to other nodes, declare the whole matrix as a independent variable!: */ 169 166 for (int i=0;i<M*N;++i) matrix[i]<<=buffer[i]; /*we use the <<= ADOLC overloaded operator to declare the independency*/ … … 172 169 MPI_Bcast(matrix,M*N,MPI_DOUBLE,0,IssmComm::GetComm()); 173 170 #endif 174 171 175 172 xDelete<IssmPDouble>(buffer); 176 173 } … … 205 202 206 203 int i; 207 204 208 205 /*Branch according to the type of variable: */ 209 206 if(type==0){ /*scalar:*/ -
issm/trunk-jpl/src/c/classes/objects/Inputs/BoolInput.cpp
r13414 r13622 62 62 /*FUNCTION BoolInput::copy{{{*/ 63 63 Object* BoolInput::copy() { 64 64 65 65 return new BoolInput(this->enum_type,this->value); 66 66 67 67 } 68 68 /*}}}*/ 69 69 70 70 /*BoolInput management*/ 71 71 /*FUNCTION BoolInput::InstanceEnum{{{*/ … … 93 93 /*FUNCTION BoolInput::SpawnResult{{{*/ 94 94 ElementResult* BoolInput::SpawnResult(int step, IssmDouble time){ 95 95 96 96 return new BoolElementResult(this->enum_type,this->value,step,time); 97 97 -
issm/trunk-jpl/src/c/classes/objects/Inputs/ControlInput.cpp
r13414 r13622 93 93 /*FUNCTION ControlInput::copy{{{*/ 94 94 Object* ControlInput::copy() { 95 95 96 96 ControlInput* output=NULL; 97 97 … … 109 109 } 110 110 /*}}}*/ 111 111 112 112 /*ControlInput management*/ 113 113 /*FUNCTION ControlInput::InstanceEnum{{{*/ -
issm/trunk-jpl/src/c/classes/objects/Inputs/DatasetInput.cpp
r13414 r13622 64 64 /*FUNCTION DatasetInput::copy{{{*/ 65 65 Object* DatasetInput::copy() { 66 66 67 67 DatasetInput* output=NULL; 68 68 … … 89 89 } 90 90 /*}}}*/ 91 91 92 92 /*DatasetInput management*/ 93 93 /*FUNCTION DatasetInput::InstanceEnum{{{*/ … … 111 111 if(index<0 || index > inputs->Size()-1) _error_("index requested (" << index << ") exceeds dataset size (" << inputs->Size() << ")"); 112 112 Input* input=(Input*)this->inputs->GetObjectByOffset(index); 113 113 114 114 input->GetInputValue(pvalue,gauss); 115 115 } -
issm/trunk-jpl/src/c/classes/objects/Inputs/DoubleInput.cpp
r13414 r13622 62 62 /*FUNCTION DoubleInput::copy{{{*/ 63 63 Object* DoubleInput::copy() { 64 64 65 65 return new DoubleInput(this->enum_type,this->value); 66 66 67 67 } 68 68 /*}}}*/ 69 69 70 70 /*DoubleInput management*/ 71 71 /*FUNCTION DoubleInput::InstanceEnum{{{*/ -
issm/trunk-jpl/src/c/classes/objects/Inputs/IntInput.cpp
r13414 r13622 57 57 /*FUNCTION IntInput::copy{{{*/ 58 58 Object* IntInput::copy() { 59 59 60 60 return new IntInput(this->enum_type,this->value); 61 61 … … 92 92 /*FUNCTION IntInput::SpawnResult{{{*/ 93 93 ElementResult* IntInput::SpawnResult(int step, IssmDouble time){ 94 94 95 95 _error_("not supported yet!"); 96 96 -
issm/trunk-jpl/src/c/classes/objects/Inputs/PentaP1Input.cpp
r13414 r13622 71 71 } 72 72 /*}}}*/ 73 73 74 74 /*PentaP1Input management*/ 75 75 /*FUNCTION PentaP1Input::copy{{{*/ 76 76 Object* PentaP1Input::copy() { 77 77 78 78 return new PentaP1Input(this->enum_type,this->values); 79 79 … … 357 357 /*FUNCTION PentaP1Input::ConstrainMin{{{*/ 358 358 void PentaP1Input::ConstrainMin(IssmDouble minimum){ 359 359 360 360 int i; 361 361 const int numnodes=6; … … 425 425 /*FUNCTION PentaP1Input::Scale{{{*/ 426 426 void PentaP1Input::Scale(IssmDouble scale_factor){ 427 427 428 428 int i; 429 429 const int numnodes=6; … … 464 464 int i; 465 465 const int numnodes=6; 466 466 467 467 if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min; 468 468 if(!xIsNan<IssmDouble>(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max; -
issm/trunk-jpl/src/c/classes/objects/Inputs/TransientInput.cpp
r13577 r13622 103 103 } 104 104 /*}}}*/ 105 105 106 106 /*TransientInput management*/ 107 107 /*FUNCTION TransientInput::InstanceEnum{{{*/ … … 216 216 /*Retrieve interpolated values for this time step: */ 217 217 Input* input=GetTimeInput(time); 218 218 219 219 /*Call input function*/ 220 220 input->GetInputDerivativeValue(p,xyz_list,gauss); … … 231 231 /*FUNCTION TransientInput::GetInputAverage{{{*/ 232 232 void TransientInput::GetInputAverage(IssmDouble* pvalue){ 233 233 234 234 IssmDouble time; 235 235 … … 242 242 /*Call input function*/ 243 243 input->GetInputAverage(pvalue); 244 244 245 245 delete input; 246 246 … … 296 296 /*Retrieve interpolated values for this time step: */ 297 297 Input* input=GetTimeInput(time); 298 298 299 299 /*Call input function*/ 300 300 input->SquareMin(psquaremin,process_units,parameters); 301 301 302 302 delete input; 303 303 … … 318 318 /*Call input function*/ 319 319 infnorm=input->InfinityNorm(); 320 320 321 321 /*Clean-up and return*/ 322 322 delete input; … … 335 335 /*Retrieve interpolated values for this time step: */ 336 336 Input* input=GetTimeInput(time); 337 337 338 338 /*Call input function*/ 339 339 max=input->Max(); 340 340 341 341 delete input; 342 342 … … 400 400 /*Call input function*/ 401 401 minabs=input->MinAbs(); 402 402 403 403 /*Clean-up and return*/ 404 404 delete input; … … 416 416 /*Retrieve interpolated values for this time step: */ 417 417 Input* input=GetTimeInput(time); 418 418 419 419 /*Call input function*/ 420 420 input->GetVectorFromInputs(vector,doflist); 421 421 422 422 delete input; 423 423 … … 430 430 int found; 431 431 int offset; 432 432 433 433 Input *input = NULL; 434 434 Input *input1 = NULL; 435 435 Input *input2 = NULL; 436 436 437 437 /*go through the timesteps, and figure out which interval we 438 438 *fall within. Then interpolate the values on this interval: */ -
issm/trunk-jpl/src/c/classes/objects/Inputs/TriaP1Input.cpp
r13414 r13622 73 73 /*FUNCTION TriaP1Input::copy{{{*/ 74 74 Object* TriaP1Input::copy() { 75 75 76 76 return new TriaP1Input(this->enum_type,this->values); 77 77 78 78 } 79 79 /*}}}*/ 80 80 81 81 /*TriaP1Input management*/ 82 82 /*FUNCTION TriaP1Input::InstanceEnum{{{*/ … … 208 208 /*FUNCTION TriaP1Input::ContrainMin{{{*/ 209 209 void TriaP1Input::ConstrainMin(IssmDouble minimum){ 210 210 211 211 int i; 212 212 const int numnodes=3; … … 276 276 /*FUNCTION TriaP1Input::Scale{{{*/ 277 277 void TriaP1Input::Scale(IssmDouble scale_factor){ 278 278 279 279 int i; 280 280 const int numnodes=3; … … 327 327 int i; 328 328 const int numnodes=3; 329 329 330 330 if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min; 331 331 if(!xIsNan<IssmDouble>(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max; -
issm/trunk-jpl/src/c/classes/objects/KML/KML_Document.cpp
r13056 r13622 125 125 } 126 126 /*}}}*/ 127 -
issm/trunk-jpl/src/c/classes/objects/KML/KML_LatLonBox.cpp
r13056 r13622 43 43 /*FUNCTION KML_LatLonBox::Echo {{{*/ 44 44 void KML_LatLonBox::Echo(){ 45 46 45 47 46 _printLine_("KML_LatLonBox:"); -
issm/trunk-jpl/src/c/classes/objects/KML/KML_Polygon.cpp
r13036 r13622 250 250 } 251 251 252 253 252 else if (!strncmp(kstri,"<",1)) 254 253 KML_Geometry::Read(fid,kstri); -
issm/trunk-jpl/src/c/classes/objects/KML/KML_StyleSelector.cpp
r13056 r13622 93 93 } 94 94 /*}}}*/ 95 -
issm/trunk-jpl/src/c/classes/objects/KML/KML_SubStyle.cpp
r13036 r13622 93 93 } 94 94 /*}}}*/ 95 -
issm/trunk-jpl/src/c/classes/objects/KML/KML_Unknown.cpp
r13056 r13622 79 79 valuei=xNew<char>(strlen(value)+1); 80 80 memcpy(valuei,value,(strlen(value)+1)*sizeof(char)); 81 81 82 82 vtoken=strtok(valuei,nl); 83 83 if(flag) _pprintString_(indent << " value: \"" << vtoken); 84 84 85 85 while (vtoken=strtok(NULL,nl)) 86 86 if(flag) _pprintString_("\n" << indent << " " << vtoken); … … 110 110 valuei=xNew<char>(strlen(value)+1); 111 111 memcpy(valuei,value,(strlen(value)+1)*sizeof(char)); 112 112 113 113 vtoken=strtok(valuei,nl); 114 114 fprintf(filout,"%s %s\n",indent,vtoken); 115 115 116 116 while (vtoken=strtok(NULL,nl)) 117 117 fprintf(filout,"%s %s\n",indent,vtoken); -
issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.cpp
r13414 r13622 52 52 int icefront_node_ids[NUMVERTICESQUA]; //initialize with largest size 53 53 int icefront_fill; 54 54 55 55 /*find parameters: */ 56 56 iomodel->Constant(&dim,MeshDimensionEnum); … … 90 90 /*Fill*/ 91 91 icefront_fill=reCast<int>(iomodel->Data(DiagnosticIcefrontEnum)[segment_width*i+segment_width-1]); 92 92 93 93 /*Ok, we have everything to build the object: */ 94 94 this->id=icefront_id; … … 104 104 this->inputs->AddInput(new IntInput(FillEnum,icefront_fill)); 105 105 this->inputs->AddInput(new IntInput(TypeEnum,in_icefront_type)); 106 106 107 107 //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this. 108 108 this->parameters=NULL; … … 111 111 this->matpar= NULL; 112 112 } 113 114 113 115 114 /*}}}*/ … … 166 165 /*FUNCTION Icefront::copy {{{*/ 167 166 Object* Icefront::copy() { 168 167 169 168 Icefront* icefront=NULL; 170 169 … … 666 665 int* doflist=NULL; 667 666 668 669 667 /*recover type: */ 670 668 inputs->GetInputValue(&type,TypeEnum); … … 672 670 /*Some checks for debugging*/ 673 671 _assert_(nodes); 674 672 675 673 /*How many nodes? :*/ 676 674 if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum) … … 678 676 else 679 677 numberofnodes=4; 680 678 681 679 /*Figure out size of doflist: */ 682 680 for(i=0;i<numberofnodes;i++){ -
issm/trunk-jpl/src/c/classes/objects/Loads/Numericalflux.cpp
r13414 r13622 175 175 _printLine_(" inputs"); 176 176 inputs->DeepEcho(); 177 177 178 178 } 179 179 /*}}}*/ … … 192 192 /*FUNCTION Numericalflux::copy {{{*/ 193 193 Object* Numericalflux::copy() { 194 194 195 195 Numericalflux* numericalflux=NULL; 196 196 … … 404 404 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j]; 405 405 } 406 406 407 407 /*Clean up and return*/ 408 408 delete gauss; -
issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.cpp
r13414 r13622 19 19 #include "../../../Container/Container.h" 20 20 /*}}}*/ 21 21 22 22 /*Element macros*/ 23 23 #define NUMVERTICES 1 … … 34 34 this->hmatpar=NULL; 35 35 this->matpar=NULL; 36 36 37 37 /*not active, not zigzagging: */ 38 38 active=0; … … 64 64 this->id=id; 65 65 this->analysis_type=in_analysis_type; 66 66 67 67 /*hooks: */ 68 68 pengrid_node_id=iomodel->nodecounter+index+1; … … 99 99 } 100 100 /*}}}*/ 101 101 102 102 /*Object virtual functions definitions:*/ 103 103 /*FUNCTION Pengrid::Echo {{{*/ … … 134 134 /*FUNCTION Icefront::copy {{{*/ 135 135 Object* Pengrid::copy() { 136 136 137 137 Pengrid* pengrid=NULL; 138 138 … … 395 395 /*recover pointers: */ 396 396 Penta* penta=(Penta*)element; 397 397 398 398 /*check that pengrid is not a clone (penalty to be added only once)*/ 399 399 if (node->IsClone()){ … … 409 409 //Recover our data: 410 410 parameters->FindParam(&penalty_lock,ThermalPenaltyLockEnum); 411 411 412 412 //Compute pressure melting point 413 413 t_pmp=matpar->TMeltingPoint(pressure); … … 421 421 new_active=0; 422 422 } 423 424 423 425 424 //Figure out stability of this penalty … … 450 449 /*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticStokes {{{*/ 451 450 ElementMatrix* Pengrid::PenaltyCreateKMatrixDiagnosticStokes(IssmDouble kmax){ 452 451 453 452 const int numdof = NUMVERTICES *NDOF4; 454 453 IssmDouble slope[2]; … … 499 498 penta->GetInputValue(&temperature,node,TemperatureEnum); 500 499 parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); 501 500 502 501 /*Compute pressure melting point*/ 503 502 t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; … … 533 532 /*FUNCTION Pengrid::PenaltyCreatePVectorMelting {{{*/ 534 533 ElementVector* Pengrid::PenaltyCreatePVectorMelting(IssmDouble kmax){ 535 534 536 535 const int numdof=NUMVERTICES*NDOF1; 537 536 IssmDouble pressure; -
issm/trunk-jpl/src/c/classes/objects/Loads/Penpair.cpp
r13414 r13622 34 34 /*FUNCTION Penpair::creation {{{*/ 35 35 Penpair::Penpair(int penpair_id, int* penpair_node_ids,int in_analysis_type){ 36 36 37 37 this->id=penpair_id; 38 38 this->analysis_type=in_analysis_type; … … 40 40 this->parameters=NULL; 41 41 this->nodes=NULL; 42 42 43 43 return; 44 44 } … … 61 61 _printLine_(" analysis_type: " << EnumToStringx(analysis_type)); 62 62 hnodes->Echo(); 63 63 64 64 return; 65 65 } … … 87 87 /*FUNCTION Penpair::copy {{{*/ 88 88 Object* Penpair::copy() { 89 89 90 90 Penpair* penpair=NULL; 91 91 … … 107 107 } 108 108 /*}}}*/ 109 109 110 110 /*Load virtual functions definitions:*/ 111 111 /*FUNCTION Penpair::Configure {{{*/ … … 264 264 /*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticMacAyealPattyn {{{*/ 265 265 ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticMacAyealPattyn(IssmDouble kmax){ 266 266 267 267 const int numdof=NUMVERTICES*NDOF2; 268 268 IssmDouble penalty_offset; … … 291 291 /*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticStokes {{{*/ 292 292 ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticStokes(IssmDouble kmax){ 293 293 294 294 const int numdof=NUMVERTICES*NDOF4; 295 295 IssmDouble penalty_offset; … … 311 311 Ke->values[5*numdof+1]=-kmax*pow((IssmDouble)10.0,penalty_offset); 312 312 Ke->values[5*numdof+5]=+kmax*pow((IssmDouble)10.0,penalty_offset); 313 313 314 314 Ke->values[2*numdof+2]=+kmax*pow((IssmDouble)10.0,penalty_offset); 315 315 Ke->values[2*numdof+6]=-kmax*pow((IssmDouble)10.0,penalty_offset); -
issm/trunk-jpl/src/c/classes/objects/Loads/Riftfront.cpp
r13414 r13622 97 97 //intialize inputs, and add as many inputs per element as requested: 98 98 this->inputs=new Inputs(); 99 99 100 100 riftfront_type=SegmentRiftfrontEnum; 101 101 riftfront_fill = reCast<int,IssmDouble>(*(iomodel->Data(RiftsRiftstructEnum)+RIFTINFOSIZE*i+7)); … … 109 109 this->inputs->AddInput(new DoubleInput(FractionIncrementEnum,riftfront_fractionincrement)); 110 110 this->inputs->AddInput(new BoolInput(SegmentOnIceShelfEnum,riftfront_shelf)); 111 111 112 112 //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this. 113 113 this->parameters=NULL; … … 115 115 this->elements= NULL; 116 116 this->matpar= NULL; 117 117 118 118 } 119 119 /*}}}*/ … … 137 137 IssmDouble friction,fractionincrement; 138 138 139 140 139 /*recover some inputs first: */ 141 140 input=(Input*)this->inputs->GetInput(FillEnum); input->GetInputValue(&fill); … … 165 164 _printLine_(" state: " << state); 166 165 _printLine_(" frozen: " << (frozen ? "true":"false")); 167 166 168 167 } 169 168 /*}}}*/ … … 195 194 /*FUNCTION Riftfront::copy {{{*/ 196 195 Object* Riftfront::copy() { 197 196 198 197 Riftfront* riftfront=NULL; 199 198 … … 234 233 riftfront->length=this->length; 235 234 riftfront->fraction=this->fraction; 236 235 237 236 return riftfront; 238 237 239 238 } 240 239 /*}}}*/ 241 240 242 241 /*Update virtual functions definitions:*/ 243 242 /*FUNCTION Riftfront::InputUpdateFromConstant(bool constant,int name) {{{*/ … … 275 274 } 276 275 /*}}}*/ 277 278 276 279 277 /*Load virtual functions definitions:*/ … … 596 594 if(this->state==OpenEnum)this->active=0; 597 595 if(this->state==ClosedEnum)this->active=1; 598 596 599 597 /*this segment is like frozen, no instability here: */ 600 598 *punstable=0; 601 599 return 1; 602 600 } 603 604 601 605 602 /*recover parameters: */ … … 723 720 /*If we are zigzag locked, same thing: */ 724 721 if(this->counter>this->penalty_lock)penetration=-1; 725 722 726 723 /*assign output pointer: */ 727 724 *ppenetration=penetration; … … 762 759 /*Now, we return penetration only if we are active!: */ 763 760 if(this->active==0)penetration=0; 764 761 765 762 /*assign output pointer: */ 766 763 *ppenetration=penetration; … … 770 767 /*FUNCTION Riftfront::PotentialUnstableConstraint {{{*/ 771 768 int Riftfront::PotentialUnstableConstraint(int* punstable){ 772 773 769 774 770 const int numnodes = 2; -
issm/trunk-jpl/src/c/classes/objects/Materials/Matdamageice.cpp
r13414 r13622 15 15 #include "../../../shared/shared.h" 16 16 #include "../../../include/include.h" 17 17 18 18 /*Matdamageice constructors and destructor*/ 19 19 /*FUNCTION Matdamageice::Matdamageice(){{{*/ … … 299 299 * return g, initial viscosity. 300 300 */ 301 301 302 302 /*output: */ 303 303 IssmDouble viscosity3d; … … 341 341 else{ 342 342 e=(n-1)/2/n; 343 343 344 344 viscosity3d=B/(2*pow(A,e)); 345 345 } … … 370 370 * return g, initial viscosity. 371 371 */ 372 372 373 373 /*output: */ 374 374 IssmDouble viscosity3d; … … 387 387 Z=GetZ(); 388 388 B=Z*GetB(); 389 389 390 390 if (n==1){ 391 391 /*Viscous behaviour! viscosity3d=B: */ … … 440 440 * return mu20, initial viscosity. 441 441 */ 442 442 443 443 /*output: */ 444 444 IssmDouble viscosity_complement; … … 468 468 else{ 469 469 e=(n-1)/(2*n); 470 470 471 471 viscosity_complement=1/(2*pow(A,e)); 472 472 } … … 480 480 _assert_(n>0); 481 481 _assert_(viscosity_complement>0); 482 482 483 483 /*Return: */ 484 484 *pviscosity_complement=viscosity_complement; … … 496 496 * return mu20, initial viscosity. 497 497 */ 498 498 499 499 /*output: */ 500 500 IssmDouble viscosity_complement; … … 524 524 else{ 525 525 e=(n-1)/(2*n); 526 526 527 527 viscosity_complement=B/(2*pow(A,e)); 528 528 } … … 536 536 _assert_(n>0); 537 537 _assert_(viscosity_complement>0); 538 538 539 539 /*Return: */ 540 540 *pviscosity_complement=viscosity_complement; … … 700 700 default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); 701 701 } 702 703 704 702 705 703 } -
issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp
r13414 r13622 15 15 #include "../../../shared/shared.h" 16 16 #include "../../../include/include.h" 17 17 18 18 /*Matice constructors and destructor*/ 19 19 /*FUNCTION Matice::Matice(){{{*/ … … 277 277 * return g, initial viscosity. 278 278 */ 279 279 280 280 /*output: */ 281 281 IssmDouble viscosity3d; … … 318 318 else{ 319 319 e=(n-1)/2/n; 320 320 321 321 viscosity3d=B/(2*pow(A,e)); 322 322 } … … 347 347 * return g, initial viscosity. 348 348 */ 349 349 350 350 /*output: */ 351 351 IssmDouble viscosity3d; … … 363 363 B=GetB(); 364 364 n=GetN(); 365 365 366 366 if (n==1){ 367 367 /*Viscous behaviour! viscosity3d=B: */ … … 416 416 * return mu20, initial viscosity. 417 417 */ 418 418 419 419 /*output: */ 420 420 IssmDouble viscosity_complement; … … 444 444 else{ 445 445 e=(n-1)/(2*n); 446 446 447 447 viscosity_complement=1/(2*pow(A,e)); 448 448 } … … 456 456 _assert_(n>0); 457 457 _assert_(viscosity_complement>0); 458 458 459 459 /*Return: */ 460 460 *pviscosity_complement=viscosity_complement; … … 621 621 } 622 622 623 624 625 623 } 626 624 /*}}}*/ -
issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp
r13414 r13622 15 15 #include "../../../include/include.h" 16 16 #include "../../../EnumDefinitions/EnumDefinitions.h" 17 17 18 18 /*Matpar constructors and destructor*/ 19 19 /*FUNCTION Matpar::Matpar() {{{*/ … … 39 39 iomodel->Constant(&this->thermal_exchange_velocity,MaterialsThermalExchangeVelocityEnum); 40 40 iomodel->Constant(&this->g,ConstantsGEnum); 41 41 42 42 iomodel->Constant(&this->hydro_CR,HydrologyCREnum); 43 43 iomodel->Constant(&this->kn,HydrologyKnEnum); … … 244 244 /*FUNCTION Matpar::GetRhoIce {{{*/ 245 245 IssmDouble Matpar::GetRhoIce(){ 246 246 247 247 return rho_ice; 248 248 } … … 323 323 /*Ouput*/ 324 324 IssmDouble temperature,waterfraction; 325 325 326 326 if(enthalpy<PureIceEnthalpy(pressure)){ 327 327 temperature=referencetemperature+enthalpy/heatcapacity; … … 343 343 /*Ouput*/ 344 344 IssmDouble enthalpy; 345 345 346 346 if(temperature<TMeltingPoint(pressure)){ 347 347 enthalpy=heatcapacity*(temperature-referencetemperature); -
issm/trunk-jpl/src/c/classes/objects/Node.cpp
r13612 r13622 167 167 _printLine_(" inputs: " << inputs); 168 168 169 170 169 } 171 170 /*}}}*/ … … 182 181 _printLine_(" inputs"); 183 182 184 185 183 } 186 184 /*}}}*/ … … 232 230 int count=0; 233 231 int count2=0; 234 232 235 233 if(approximation_enum==NoneApproximationEnum){ 236 234 if(setenum==GsetEnum)for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=indexing.gdoflist[i]; … … 294 292 int count=0; 295 293 int count2=0; 296 294 297 295 if(approximation_enum==NoneApproximationEnum){ 298 296 if(setenum==GsetEnum)for(i=0;i<this->indexing.gsize;i++) outdoflist[i]=i; … … 468 466 /*g set: */ 469 467 pv_g->SetValue(indexing.gdoflist[i],gvalue,INS_VAL); 470 468 471 469 /*f set: */ 472 470 value=(IssmDouble)this->indexing.f_set[i]; … … 478 476 479 477 } 480 481 478 482 479 } … … 500 497 } 501 498 } 502 499 503 500 /*Add values into constraint vector: */ 504 501 ys->SetValues(this->indexing.ssize,this->indexing.sdoflist,values,INS_VAL); … … 508 505 xDelete<IssmDouble>(values); 509 506 510 511 507 } 512 508 /*}}}*/ … … 533 529 /*FUNCTION Node::FreezeDof{{{*/ 534 530 void Node::FreezeDof(int dof){ 535 531 536 532 DofInSSet(dof-1); //with 0 displacement for this dof. 537 533 … … 562 558 /*Get number of degrees of freedom in a node, for a certain set (g,f or s-set) 563 559 *and for a certain approximation type: */ 564 560 565 561 int i; 566 562 int numdofs=0; … … 639 635 /*FUNCTION Node::IsClone {{{*/ 640 636 int Node::IsClone(){ 641 637 642 638 return indexing.clone; 643 639 … … 668 664 /*FUNCTION Node::IsFloating {{{*/ 669 665 int Node::IsFloating(){ 670 666 671 667 bool onshelf; 672 668 … … 872 868 if(setenum==FsetEnum) this->indexing.InitSet(setenum); 873 869 if(setenum==SsetEnum) this->indexing.InitSet(setenum); 874 870 875 871 /*For clone nodfs, don't distribute dofs, we will get them from another cpu in UpdateCloneDofs!*/ 876 872 if(indexing.clone){ … … 905 901 /*FUNCTION Node::OffsetDofs{{{*/ 906 902 void Node::OffsetDofs(int dofcount,int setenum){ 907 903 908 904 int i; 909 905 910 906 if(indexing.clone){ 911 907 /*This node is a clone, don't off_set the dofs!: */ … … 980 976 981 977 int my_rank; 982 978 983 979 /*recover my_rank:*/ 984 980 my_rank=IssmComm::GetRank(); -
issm/trunk-jpl/src/c/classes/objects/Params/BoolParam.cpp
r13414 r13622 65 65 /*FUNCTION BoolParam::copy{{{*/ 66 66 Object* BoolParam::copy() { 67 67 68 68 return new BoolParam(this->enum_type,this->value); 69 69 -
issm/trunk-jpl/src/c/classes/objects/Params/DataSetParam.cpp
r13608 r13622 66 66 /*FUNCTION DataSetParam::copy{{{*/ 67 67 Object* DataSetParam::copy() { 68 68 69 69 return new DataSetParam(this->enum_type,this->value); 70 70 -
issm/trunk-jpl/src/c/classes/objects/Params/DoubleMatArrayParam.cpp
r13414 r13622 77 77 xDelete<IssmDouble>(matrix); 78 78 } 79 79 80 80 xDelete<IssmDouble*>(array); 81 81 return; … … 100 100 int m,n; 101 101 IssmDouble* matrix=NULL; 102 102 103 103 _printLine_("DoubleMatArrayParam:"); 104 104 _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")"); … … 130 130 /*FUNCTION DoubleMatArrayParam::copy{{{*/ 131 131 Object* DoubleMatArrayParam::copy() { 132 132 133 133 return new DoubleMatArrayParam(this->enum_type,this->array, this->M, this->mdim_array,this->ndim_array); 134 134 … … 150 150 int* out_ndim_array=NULL; 151 151 152 153 152 out_M=this->M; 154 153 if(out_M){ … … 181 180 } 182 181 183 184 182 /*Assign output pointers:*/ 185 183 if(pout_M) *pout_M=out_M; … … 216 214 this->mdim_array=xNew<int>(M); 217 215 this->ndim_array=xNew<int>(M); 218 216 219 217 xMemCpy<int>(this->mdim_array,in_mdim_array,M); 220 218 xMemCpy<int>(this->ndim_array,in_ndim_array,M); -
issm/trunk-jpl/src/c/classes/objects/Params/DoubleMatParam.cpp
r13414 r13622 58 58 59 59 int i,j; 60 60 61 61 _printLine_("DoubleMatParam:"); 62 62 _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")"); … … 81 81 /*FUNCTION DoubleMatParam::copy{{{*/ 82 82 Object* DoubleMatParam::copy() { 83 83 84 84 return new DoubleMatParam(this->enum_type,this->value,this->M,this->N); 85 85 -
issm/trunk-jpl/src/c/classes/objects/Params/DoubleParam.cpp
r13414 r13622 62 62 /*FUNCTION DoubleParam::copy{{{*/ 63 63 Object* DoubleParam::copy() { 64 64 65 65 return new DoubleParam(this->enum_type,this->value); 66 66 -
issm/trunk-jpl/src/c/classes/objects/Params/DoubleVecParam.cpp
r13414 r13622 57 57 58 58 int i; 59 59 60 60 _printLine_("DoubleVecParam:"); 61 61 _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")"); … … 78 78 /*FUNCTION DoubleVecParam::copy{{{*/ 79 79 Object* DoubleVecParam::copy() { 80 80 81 81 return new DoubleVecParam(this->enum_type,this->values,this->M); 82 82 -
issm/trunk-jpl/src/c/classes/objects/Params/FileParam.cpp
r13414 r13622 65 65 /*FUNCTION FileParam::copy{{{*/ 66 66 Object* FileParam::copy() { 67 67 68 68 return new FileParam(this->enum_type,this->value); 69 69 -
issm/trunk-jpl/src/c/classes/objects/Params/IntMatParam.cpp
r13414 r13622 58 58 59 59 int i,j; 60 60 61 61 _printLine_("IntMatParam:"); 62 62 _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")"); … … 81 81 /*FUNCTION IntMatParam::copy{{{*/ 82 82 Object* IntMatParam::copy() { 83 83 84 84 return new IntMatParam(this->enum_type,this->value,this->M,this->N); 85 85 -
issm/trunk-jpl/src/c/classes/objects/Params/IntParam.cpp
r13414 r13622 65 65 /*FUNCTION IntParam::copy{{{*/ 66 66 Object* IntParam::copy() { 67 67 68 68 return new IntParam(this->enum_type,this->value); 69 69 -
issm/trunk-jpl/src/c/classes/objects/Params/IntVecParam.cpp
r13414 r13622 73 73 74 74 int i; 75 75 76 76 _printLine_("IntVecParam:"); 77 77 _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")"); … … 94 94 /*FUNCTION IntVecParam::copy{{{*/ 95 95 Object* IntVecParam::copy() { 96 96 97 97 return new IntVecParam(this->enum_type,this->values,this->M); 98 98 -
issm/trunk-jpl/src/c/classes/objects/Params/MatrixParam.cpp
r13414 r13622 73 73 /*FUNCTION MatrixParam::copy{{{*/ 74 74 Object* MatrixParam::copy() { 75 75 76 76 return new MatrixParam(this->enum_type,this->value); 77 77 … … 97 97 /*FUNCTION MatrixParam::SetValue{{{*/ 98 98 void MatrixParam::SetValue(Matrix<IssmDouble>* matrix){ 99 99 100 100 /*avoid leak: */ 101 101 xdelete(&value); 102 102 103 103 /*copy: */ 104 104 value=matrix->Duplicate(); -
issm/trunk-jpl/src/c/classes/objects/Params/StringArrayParam.cpp
r13414 r13622 46 46 } 47 47 else value=NULL; 48 48 49 49 } 50 50 /*}}}*/ 51 51 /*FUNCTION StringArrayParam::~StringArrayParam(){{{*/ 52 52 StringArrayParam::~StringArrayParam(){ 53 53 54 54 int i; 55 55 56 56 char* string=NULL; 57 57 for(i=0;i<this->numstrings;i++){ … … 95 95 /*FUNCTION StringArrayParam::copy{{{*/ 96 96 Object* StringArrayParam::copy() { 97 97 98 98 return new StringArrayParam(this->enum_type,this->value,this->numstrings); 99 99 … … 104 104 /*FUNCTION StringArrayParam::GetParameterValue{{{*/ 105 105 void StringArrayParam::GetParameterValue(char*** pstringarray,int* pM){ 106 106 107 107 int i; 108 108 char** outstrings=NULL; … … 140 140 /*FUNCTION StringArrayParam::SetValue{{{*/ 141 141 void StringArrayParam::SetValue(char** stringarray,int M){ 142 142 143 143 int i; 144 144 char *string = NULL; -
issm/trunk-jpl/src/c/classes/objects/Params/StringParam.cpp
r13414 r13622 33 33 xMemCpy<char>(value,in_value,(strlen(in_value)+1)); 34 34 35 36 35 } 37 36 /*}}}*/ … … 67 66 /*FUNCTION StringParam::copy{{{*/ 68 67 Object* StringParam::copy() { 69 68 70 69 return new StringParam(this->enum_type,this->value); 71 70 … … 76 75 /*FUNCTION StringParam::GetParameterValue{{{*/ 77 76 void StringParam::GetParameterValue(char** pstring){ 78 77 79 78 char* outstring=NULL; 80 79 int stringsize; … … 96 95 /*FUNCTION StringParam::SetValue{{{*/ 97 96 void StringParam::SetValue(char* string){ 98 97 99 98 int stringsize; 100 99 101 100 /*avoid leak: */ 102 101 xDelete<char>(this->value); -
issm/trunk-jpl/src/c/classes/objects/Params/TransientParam.cpp
r13414 r13622 62 62 63 63 int i,j; 64 64 65 65 _printLine_("TransientParam:"); 66 66 _printLine_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")"); … … 83 83 /*FUNCTION TransientParam::copy{{{*/ 84 84 Object* TransientParam::copy() { 85 85 86 86 return new TransientParam(this->enum_type,this->values,this->timesteps,this->N); 87 87 -
issm/trunk-jpl/src/c/classes/objects/Params/VectorParam.cpp
r13414 r13622 75 75 /*FUNCTION VectorParam::copy{{{*/ 76 76 Object* VectorParam::copy() { 77 77 78 78 return new VectorParam(this->enum_type,this->value); 79 79 … … 103 103 /*avoid leak: */ 104 104 xdelete(&value); 105 105 106 106 /*copy: */ 107 107 value=vector->Duplicate(); -
issm/trunk-jpl/src/c/classes/objects/Vertex.cpp
r13612 r13622 132 132 /*retrieve current pid*/ 133 133 int pidcount=*ppidcount; 134 134 135 135 /*This vertex is a clone! Don't distribute pids, it will get them from another cpu!*/ 136 136 if(this->clone) return; … … 146 146 /*FUNCTION Vertex::OffsetPids{{{*/ 147 147 void Vertex::OffsetPids(int pidcount){ 148 148 149 149 /*This vertex is a clone, don't offset the pids*/ 150 150 if(this->clone) return; … … 179 179 180 180 int my_rank; 181 181 182 182 /*recover my_rank:*/ 183 183 my_rank=IssmComm::GetRank(); … … 194 194 } 195 195 /*}}}*/ 196 -
issm/trunk-jpl/src/c/io/Disk/WriteLockFile.cpp
r13612 r13622 13 13 /*recover my_rank:*/ 14 14 my_rank=IssmComm::GetRank(); 15 15 16 16 /* output: */ 17 17 FILE* fid=NULL; -
issm/trunk-jpl/src/c/io/Disk/pfopen.cpp
r13606 r13622 16 16 17 17 FILE* fid=NULL; 18 18 19 19 /*Open handle to data on disk: */ 20 20 fid=fopen(filename,format); … … 23 23 return fid; 24 24 } 25 -
issm/trunk-jpl/src/c/io/PrintfFunction.cpp
r13612 r13622 19 19 //variable list of arguments 20 20 va_list args; 21 21 22 22 /*recover my_rank:*/ 23 23 my_rank=IssmComm::GetRank(); … … 54 54 int PrintfFunction(const string & message){ 55 55 int my_rank; 56 56 57 57 /*recover my_rank:*/ 58 58 my_rank=IssmComm::GetRank(); -
issm/trunk-jpl/src/c/matlab/io/CheckNumMatlabArguments.cpp
r13056 r13622 9 9 #endif 10 10 11 12 11 #include "../../shared/Exceptions/exceptions.h" 13 12 #include "../../include/include.h" … … 15 14 16 15 int CheckNumMatlabArguments(int nlhs,int NLHS, int nrhs,int NRHS, const char* __FUNCT__, void (*function)( void )){ 17 16 18 17 /*checks on arguments on the matlab side: */ 19 18 if (nrhs==0 && nlhs==0) { -
issm/trunk-jpl/src/c/matlab/io/FetchMatlabData.cpp
r13378 r13622 50 50 _error_("Input parameter of class " << mxGetClassName(dataref) << " not supported yet"); 51 51 } 52 52 53 53 /*Assign output pointers:*/ 54 54 *pmatrix=outmatrix; … … 97 97 _error_("Input parameter of class " << mxGetClassName(dataref) << " not supported yet"); 98 98 } 99 99 100 100 /*Assign output pointers:*/ 101 101 *pmatrix=outmatrix; … … 252 252 _error_("Input parameter of class " << mxGetClassName(dataref) << " not supported yet"); 253 253 } 254 254 255 255 /*Assign output pointers:*/ 256 256 *pmatrix=outmatrix; … … 390 390 char* outstring=NULL; 391 391 392 393 392 /*Ok, the string should be coming directly from the matlab workspace: */ 394 393 if (!mxIsClass(dataref,"char")){ … … 398 397 /*Recover the string:*/ 399 398 int stringlen; 400 399 401 400 stringlen = mxGetM(dataref)*mxGetN(dataref)+1; 402 401 outstring =xNew<char>(stringlen); … … 430 429 _error_("Input parameter of class " << mxGetClassName(dataref) << " not supported yet"); 431 430 } 432 431 433 432 /*Assign output pointers:*/ 434 433 *pmatrix=outmatrix; -
issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToDoubleMatrix.cpp
r13056 r13622 31 31 rows=mxGetM(mxmatrix); 32 32 cols=mxGetN(mxmatrix); 33 33 34 34 if(rows*cols){ 35 35 matrix=xNewZeroInit<double>(rows*cols); … … 55 55 rows=mxGetM(mxmatrix); 56 56 cols=mxGetN(mxmatrix); 57 57 58 58 /*Create serial matrix: */ 59 59 if(rows*cols){ -
issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToMatrix.cpp
r13216 r13622 19 19 #include "../../include/include.h" 20 20 #include "../../toolkits/toolkits.h" 21 21 22 22 /*}}}*/ 23 23 … … 35 35 matrix->smatrix=MatlabMatrixToSeqMat(mxmatrix); 36 36 #endif 37 37 38 38 return matrix; 39 39 } -
issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToPetscMat.cpp
r12850 r13622 2 2 * \brief: convert a sparse or dense matlab matrix to a serial Petsc matrix: 3 3 */ 4 5 4 6 5 #ifdef HAVE_CONFIG_H … … 51 50 int *idxm = NULL; 52 51 int *idxn = NULL; 53 52 54 53 /*Ok, first check if we are dealing with a sparse or full matrix: */ 55 54 if (mxIsSparse(mxmatrix)){ -
issm/trunk-jpl/src/c/matlab/io/MatlabNArrayToNArray.cpp
r13378 r13622 2 2 * \brief: convert a sparse or dense matlab n-dimensional array to cpp n-dimensional array 3 3 */ 4 5 4 6 5 #ifdef HAVE_CONFIG_H … … 68 67 /*Dealing with dense matrix: recover pointer and size: */ 69 68 mxmatrix_ptr=(double*)mxGetPr(mxmatrix); 70 69 71 70 /*Create serial matrix: */ 72 71 matrix=xNewZeroInit<double>(numel); … … 141 140 /*Dealing with dense matrix: recover pointer and size: */ 142 141 mxmatrix_ptr=(bool*)mxGetData(mxmatrix); 143 142 144 143 /*Create serial matrix: */ 145 144 matrix=xNew<bool>(numel); … … 212 211 /*Dealing with dense matrix: recover pointer and size: */ 213 212 mxmatrix_ptr=mxGetChars(mxmatrix); 214 213 215 214 /*Create serial matrix: */ 216 215 matrix=xNew<char>(numel+1); -
issm/trunk-jpl/src/c/matlab/io/MatlabVectorToDoubleVector.cpp
r13056 r13622 2 2 * \brief: convert a sparse or dense matlab vector to a serial vector: 3 3 */ 4 5 4 6 5 #ifdef HAVE_CONFIG_H … … 9 8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 10 9 #endif 11 12 10 13 11 #include <string.h> … … 44 42 cols=mxGetN(mxvector); 45 43 nnz=mxGetNzmax(mxvector); 46 44 47 45 /*Check that input is actualy a vector*/ 48 46 if (cols!=1) _error_("input vector of size " << rows << "x" << cols << " should have only one column"); -
issm/trunk-jpl/src/c/matlab/io/MatlabVectorToPetscVec.cpp
r12850 r13622 2 2 * \brief: convert a sparse or dense matlab vector to a serial Petsc vector: 3 3 */ 4 5 4 6 5 #ifdef HAVE_CONFIG_H … … 25 24 int dummy; 26 25 PetscVec* vector=new PetscVec(); 27 26 28 27 MatlabVectorToPetscVec(&vector->vector,&dummy, mxvector); 29 28 … … 51 50 /*petsc indices: */ 52 51 int* idxm=NULL; 53 52 54 53 /*Ok, first check if we are dealing with a sparse or full vector: */ 55 54 if (mxIsSparse(mxvector)){ -
issm/trunk-jpl/src/c/matlab/io/MatlabVectorToVector.cpp
r13216 r13622 19 19 #include "../../include/include.h" 20 20 #include "../../toolkits/toolkits.h" 21 21 22 22 /*}}}*/ 23 23 … … 35 35 vector->svector=MatlabVectorToSeqVec(mxvector); 36 36 #endif 37 37 38 38 return vector; 39 39 } -
issm/trunk-jpl/src/c/matlab/io/PrintfFunction.cpp
r13612 r13622 19 19 int string_size; 20 20 int my_rank; 21 21 22 22 /*recover my_rank:*/ 23 23 my_rank=IssmComm::GetRank(); -
issm/trunk-jpl/src/c/matlab/io/WriteMatlabData.cpp
r13447 r13622 20 20 mxArray *dataref = NULL; 21 21 double *tmatrix = NULL; 22 22 23 23 if(matrix){ 24 24 /*create the matlab matrixwith Matlab's memory manager */ … … 67 67 /*FUNCTION WriteData(mxArray** pdataref,double* vector, int M){{{*/ 68 68 void WriteData(mxArray** pdataref,double* vector, int M){ 69 69 70 70 mxArray* dataref = NULL; 71 71 double* vector_matlab = NULL; … … 225 225 /*FUNCTION WriteData(mxArray** pdataref,SeqMat<double>* matrix){{{*/ 226 226 void WriteData(mxArray** pdataref,SeqMat<double>* matrix){ 227 227 228 228 int i,j; 229 229 int rows,cols; … … 231 231 double *matrix_ptr = NULL; 232 232 double *tmatrix_ptr = NULL; 233 233 234 234 if(matrix){ 235 235 236 236 matrix_ptr=matrix->ToSerial(); 237 237 matrix->GetSize(&rows,&cols); … … 244 244 } 245 245 } 246 246 247 247 /*create matlab matrix: */ 248 248 dataref=mxCreateDoubleMatrix(0,0,mxREAL); … … 263 263 /*FUNCTION WriteData(mxArray** pdataref,SeqVec<double>* vector){{{*/ 264 264 void WriteData(mxArray** pdataref,SeqVec<double>* vector){ 265 265 266 266 mxArray* dataref=NULL; 267 267 double* vector_ptr=NULL; 268 268 double* vector_matlab=NULL; 269 269 int rows; 270 270 271 271 if(vector){ 272 272 /*call toolkit routine: */ 273 273 vector_ptr=vector->ToMPISerial(); 274 274 vector->GetSize(&rows); 275 275 276 276 /*now create the matlab vector with Matlab's memory manager */ 277 277 vector_matlab=(double*)mxMalloc(rows*sizeof(double)); … … 299 299 mxArray* field = NULL; 300 300 301 302 301 /*Convert field*/ 303 302 WriteData(&field,fieldpointer,M,N); -
issm/trunk-jpl/src/c/matlab/io/mxGetAssignedField.cpp
r12013 r13622 20 20 //output 21 21 mxArray* mxfield=NULL; 22 22 23 23 //input 24 24 mxArray *inputs[2]; -
issm/trunk-jpl/src/c/modules/AddExternalResultx/AddExternalResultx.cpp
r13380 r13622 8 8 #include "../../io/io.h" 9 9 #include "../../classes/objects/objects.h" 10 10 11 11 void AddExternalResultx( DataSet* results, int enumtype, IssmDouble value){ 12 12 /* Add new result in into results*/ -
issm/trunk-jpl/src/c/modules/AverageFilterx/AverageFilterx.cpp
r12450 r13622 11 11 12 12 int AverageFilterx(double** pimageout,double* image, int lines,int samps,int smooth){ 13 13 14 14 double temp; 15 15 int numvalues; -
issm/trunk-jpl/src/c/modules/AverageOntoPartitionx/AverageOntoPartitionx.cpp
r13216 r13622 18 18 #include "../modules.h" 19 19 20 21 20 void AverageOntoPartitionx(double** paverage, Elements* elements, Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,double* vertex_response){ 22 21 -
issm/trunk-jpl/src/c/modules/Chacox/Chacox.cpp
r12511 r13622 21 21 { 22 22 #ifdef _HAVE_CHACO_ //only works if Chaco library has been compiled in. 23 24 23 25 24 extern int Using_Main; /* is main routine being called? */ 26 25 extern char *PARAMS_FILENAME; /* name of file with parameter updates */ … … 62 61 int i,tvwgt; 63 62 double tgoal; 64 65 66 63 67 64 if (DEBUG_TRACE > 0) { … … 187 184 _printLine_("<Leaving main>"); 188 185 } 189 186 190 187 return(0); 191 192 188 193 189 #else //ifdef _HAVE_CHACO_ -
issm/trunk-jpl/src/c/modules/Chacox/chaco_seconds.cpp
r11236 r13622 8 8 #endif 9 9 10 double chaco_seconds(void){ 10 11 11 12 double chaco_seconds(void){13 14 12 double curtime; 15 13 … … 22 20 *is defined in the <sys/time.h> and <sys/resource.h> header files. Leaving it 23 21 *for reference in case we have a problem here in the future*/ 24 22 25 23 getrusage(RUSAGE_SELF, &rusage); 26 24 curtime = ((rusage.ru_utime.tv_sec + rusage.ru_stime.tv_sec) + -
issm/trunk-jpl/src/c/modules/Chacox/input_parse.cpp
r12511 r13622 4 4 5 5 #include "./Chacox.h" 6 7 6 8 7 #undef __FUNCT__ -
issm/trunk-jpl/src/c/modules/Chacox/user_params.cpp
r9320 r13622 14 14 #endif 15 15 16 17 16 #ifdef _HAVE_CHACO_ //only works if dakota library has been compiled in. 18 17 … … 21 20 #define TRUE 1 22 21 #define FALSE 0 23 24 22 25 23 /* Input and ouput control parameters */ … … 34 32 int PROMPT = FALSE; /* Prompt for input? (TRUE/FALSE) */ 35 33 int PRINT_HEADERS = FALSE; /* Print pretty output headers (TRUE/FALSE) */ 36 37 34 38 35 /* Eigenvector calculation parameters */ … … 57 54 int TIME_KERNELS = FALSE; /* Time numerical kernels? (TRUE/FALSE) */ 58 55 59 60 56 /* Other parameters for spectral methods */ 61 57 … … 69 65 int OPT3D_NTRIES = 5; /* # local opts to look for global min in opt3d */ 70 66 71 72 67 /* Kernighan--Lin/Fiduccia--Mattheyses parameters */ 73 68 … … 78 73 int KL_UNDO_LIST = TRUE; /* Only resort changed vtxs? (TRUE/FALSE) */ 79 74 double KL_IMBALANCE = 0.0; /* Fractional imbalance allowed by KL */ 80 81 75 82 76 /* Coarsening parameters */ … … 91 85 int KL_ONLY_BNDY = TRUE; /* Start moving vtxs on boundary? (TRUE/FALSE) */ 92 86 93 94 87 /* Parameters for post-processing options */ 95 88 … … 97 90 int INTERNAL_VERTICES = FALSE; /* ... to up internal vtxs? (TRUE/FALSE) */ 98 91 int REFINE_MAP = FALSE; /* ... to reduce hops? (TRUE/FALSE) */ 99 100 92 101 93 /* Architecture and simulator parameters */ … … 129 121 char *PARAMS_FILENAME = "User_Params"; /* File of parameter changes */ 130 122 131 132 123 /* Parameters that control debugging output */ 133 124 … … 151 142 int DEBUG_MACH_PARAMS = 0;/* Print computed machine params? (0..1) */ 152 143 153 154 144 #endif //ifdef _HAVE_CHACO_ -
issm/trunk-jpl/src/c/modules/ComputeBasalStressx/ComputeBasalStressx.cpp
r13216 r13622 38 38 /*Assign output pointers: */ 39 39 *psigma=sigma; 40 40 41 41 } -
issm/trunk-jpl/src/c/modules/ComputeStrainRatex/ComputeStrainRatex.cpp
r13216 r13622 38 38 /*Assign output pointers: */ 39 39 *peps=eps; 40 40 41 41 } -
issm/trunk-jpl/src/c/modules/ConfigureObjectsx/ConfigureObjectsx.cpp
r12515 r13622 24 24 /*Get analysis type: */ 25 25 parameters->FindParam(&configuration_type,ConfigurationTypeEnum); 26 26 27 27 if(VerboseMProcessor()) _pprintLine_(" Configuring elements..."); 28 28 for (i=0;i<elements->Size();i++){ … … 44 44 } 45 45 } 46 46 47 47 if(VerboseMProcessor()) _pprintLine_(" Configuring materials..."); 48 48 for (i=0;i<materials->Size();i++){ -
issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStatex.cpp
r13412 r13622 41 41 converged=1; 42 42 } 43 43 44 44 /*Assign output pointers: */ 45 45 *pconverged=converged; -
issm/trunk-jpl/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp
r13590 r13622 15 15 int RiftIsPresent(Loads* loads,int configuration_type){ 16 16 17 18 17 int i; 19 18 20 19 int found=0; 21 20 int mpi_found=0; … … 50 49 RiftConstrain(&num_unstable_constraints,loads,configuration_type); 51 50 if(num_unstable_constraints==0)converged=1; 52 51 53 52 if(RiftIsFrozen(loads,configuration_type)){ 54 53 converged=1; … … 69 68 70 69 int i; 71 70 72 71 /* generic object pointer: */ 73 72 Riftfront* riftfront=NULL; … … 77 76 int sum_num_unstable_constraints; 78 77 int num_unstable_constraints=0; 79 78 80 79 /*Enforce constraints: */ 81 80 for (i=0;i<loads->Size();i++){ … … 100 99 num_unstable_constraints=sum_num_unstable_constraints; 101 100 #endif 102 101 103 102 /*Assign output pointers: */ 104 103 *pnum_unstable_constraints=num_unstable_constraints; … … 110 109 111 110 int i; 112 111 113 112 /* generic object pointer: */ 114 113 Load* load=NULL; … … 121 120 122 121 if (RiftfrontEnum==loads->GetEnum(i)){ 123 122 124 123 load=(Load*)loads->GetObjectByOffset(i); 125 124 if(load->InAnalysis(configuration_type)){ … … 133 132 } 134 133 } 135 134 136 135 /*Is there just one found? that would mean we have frozen! : */ 137 136 #ifdef _HAVE_MPI_ … … 148 147 149 148 int i; 150 149 151 150 /* generic object pointer: */ 152 151 Load* load=NULL; … … 160 159 load=(Load*)loads->GetObjectByOffset(i); 161 160 if(load->InAnalysis(configuration_type)){ 162 161 163 162 riftfront=(Riftfront*)load; 164 163 riftfront->FreezeConstraints(); … … 176 175 177 176 int i; 178 177 179 178 Riftfront* riftfront=NULL; 180 179 int found=0; … … 207 206 int RiftIsPreStable(Loads* loads){ 208 207 209 210 208 int i; 211 209 212 210 Riftfront* riftfront=NULL; 213 211 int found=0; … … 246 244 void RiftSetPreStable(Loads* loads){ 247 245 248 249 246 int i; 250 247 251 248 Riftfront* riftfront=NULL; 252 249 int found=0; … … 268 265 269 266 int i; 270 267 271 268 /* generic object pointer: */ 272 269 Riftfront* riftfront=NULL; … … 275 272 int sum_num_unstable_constraints; 276 273 int num_unstable_constraints=0; 277 274 278 275 /*Enforce constraints: */ 279 276 for (i=0;i<loads->Size();i++){ … … 294 291 num_unstable_constraints=sum_num_unstable_constraints; 295 292 #endif 296 293 297 294 /*Assign output pointers: */ 298 295 *pnum_unstable_constraints=num_unstable_constraints; … … 304 301 305 302 int i; 306 303 307 304 /* generic object pointer: */ 308 305 Riftfront* riftfront=NULL; … … 346 343 347 344 int i; 348 345 349 346 /* generic object pointer: */ 350 347 Riftfront* riftfront=NULL; -
issm/trunk-jpl/src/c/modules/ConstraintsStatex/ThermalIsPresent.cpp
r13590 r13622 27 27 } 28 28 } 29 29 30 30 #ifdef _HAVE_MPI_ 31 31 MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,IssmComm::GetComm() ); -
issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.cpp
r13229 r13622 7 7 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 8 8 #endif 9 10 9 11 10 #include "./ContourToMeshx.h" … … 27 26 num=_NUMTHREADS_; 28 27 #endif 29 30 28 31 29 /*output: */ -
issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp
r13229 r13622 12 12 13 13 void* ContourToMeshxt(void* vpthread_handle){ 14 14 15 15 int noerr=1; 16 16 -
issm/trunk-jpl/src/c/modules/ContourToNodesx/ContourToNodesx.cpp
r13229 r13622 34 34 /*Assign output pointers: */ 35 35 *pflags=flags; 36 36 37 37 return 1; 38 38 } -
issm/trunk-jpl/src/c/modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp
r13073 r13622 16 16 IssmDouble *scalar_list = NULL; 17 17 IssmDouble scalar; 18 19 18 20 19 /*Retrieve some parameters*/ -
issm/trunk-jpl/src/c/modules/CostFunctionx/CostFunctionx.cpp
r13073 r13622 21 21 /*output: */ 22 22 IssmDouble J,Jplus; 23 23 24 24 /*Recover parameters*/ 25 25 parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); -
issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp
r13216 r13622 11 11 12 12 void CreateJacobianMatrixx(Matrix<IssmDouble>** pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,IssmDouble kmax){ 13 13 14 14 int i,connectivity; 15 15 int numberofdofspernode; … … 30 30 /*Initialize Jacobian Matrix*/ 31 31 Jff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode); 32 32 33 33 /*Create and assemble matrix*/ 34 34 for(i=0;i<elements->Size();i++){ -
issm/trunk-jpl/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp
r13216 r13622 13 13 14 14 int i; 15 15 16 16 /*intermediary: */ 17 17 int numberofdofs; -
issm/trunk-jpl/src/c/modules/DakotaResponsesx/DakotaResponsesx.cpp
r13612 r13622 16 16 #include "../../toolkits/toolkits.h" 17 17 #include "../modules.h" 18 19 18 20 19 void DakotaResponsesx(double* d_responses,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){ … … 43 42 //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. 44 43 //because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */ 45 44 46 45 for(i=0;i<numresponsedescriptors;i++){ 47 46 48 47 flag=DescriptorIndex(root,&index,responses_descriptors[i]); 49 48 … … 77 76 //Responsex(responses_pointer,elements,nodes, vertices,loads,materials, parameters,root,process_units); 78 77 Responsex(&femmodel_response,elements,nodes, vertices,loads,materials, parameters,root,process_units,0);//0 is the index for weights 79 78 80 79 if(my_rank==0){ 81 80 /*plug response: */ … … 93 92 } 94 93 else if (flag==RegularEnum){ 95 94 96 95 /*perfectly normal response function: */ 97 96 Responsex(&femmodel_response,elements,nodes, vertices,loads,materials, parameters,root,process_units,0);//0 is the weight index … … 108 107 } 109 108 110 111 109 /*Synthesize echo: {{{*/ 112 110 if(my_rank==0){ -
issm/trunk-jpl/src/c/modules/Delta18oParameterizationx/Delta18oParameterizationx.cpp
r12748 r13622 15 15 int i; 16 16 Element* element = NULL; 17 17 18 18 /*Compute temperature and precipitation fields: */ 19 19 for(i=0;i<elements->Size();i++){ -
issm/trunk-jpl/src/c/modules/ElementConnectivityx/ElementConnectivityx.cpp
r12450 r13622 41 41 * Once we get the neighbouring elements, figure out if they share a segment with the current element. If so, 42 42 * plug them in the connectivity, unless they are already there.: */ 43 43 44 44 for(n=0;n<nel;n++){ 45 45 … … 47 47 48 48 for(i=0;i<3;i++){ 49 49 50 50 node=(int)*(elements+n*3+i); //already matlab indexed, elements comes directly from the workspace. 51 51 index=node-1; … … 58 58 connectedelement=*(nodeconnectivity+width*index+j); 59 59 connectedelementindex=(int)(connectedelement-1); //go from matlab indexing to c indexing. 60 60 61 61 if(hascommondedge(elements+n*3+0,elements+connectedelementindex*3+0)){ 62 62 /*Ok, this connected element has a commond edge with element, plug it into elementconnectivity, unless … … 80 80 *pelementconnectivity=elementconnectivity; 81 81 } 82 83 82 84 83 int hascommondedge(double* element1,double* element2){ -
issm/trunk-jpl/src/c/modules/ElementResponsex/ElementResponsex.cpp
r13612 r13622 11 11 12 12 void ElementResponsex( IssmDouble* presponse, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,int response_enum,bool process_units){ 13 14 13 15 14 int my_rank; -
issm/trunk-jpl/src/c/modules/Exp2Kmlx/Exp2Kmlx.cpp
r13056 r13622 295 295 return(iret); 296 296 } 297 -
issm/trunk-jpl/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp
r13216 r13622 27 27 gsize=nodes->NumberOfDofs(configuration_type,GsetEnum); 28 28 if (gsize==0) _error_("Allocating a Vec of size 0 as gsize=0 for configuration: " << EnumToStringx(configuration_type)); 29 29 30 30 /*Initialize solution: */ 31 31 solution=new Vector<IssmDouble>(gsize); 32 32 33 33 /*Go through elements and plug solution: */ 34 34 for (i=0;i<elements->Size();i++){ -
issm/trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp
r13216 r13622 37 37 38 38 void GetVectorFromControlInputsx( IssmDouble** pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, const char* data){ 39 39 40 40 /*output: */ 41 41 IssmDouble* vector=NULL; 42 42 43 43 /*intermediary: */ 44 44 Vector<IssmDouble>* vec_vector=NULL; -
issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp
r13216 r13622 42 42 43 43 void GetVectorFromInputsx( IssmDouble** pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name, int type){ 44 44 45 45 /*output: */ 46 46 IssmDouble* vector=NULL; 47 47 48 48 /*intermediary: */ 49 49 Vector<IssmDouble>* vec_vector=NULL; -
issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp
r13216 r13622 19 19 Vector<IssmDouble> *gradient = NULL; 20 20 Vector<IssmDouble> **gradient_list = NULL; 21 21 22 22 /*retrieve some parameters: */ 23 23 parameters->FindParam(&num_controls,InversionNumControlParametersEnum); _assert_(num_controls); -
issm/trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp
r13590 r13622 19 19 Vector<IssmDouble>* vec_old_floatingice = NULL; 20 20 Element* element = NULL; 21 21 22 22 if(VerboseModule()) _pprintLine_(" Migrating grounding line"); 23 23 24 24 /*retrieve parameters: */ 25 25 parameters->FindParam(&migration_style,GroundinglineMigrationEnum); … … 121 121 Element* element = NULL; 122 122 123 124 123 /*recover parameters: */ 125 124 parameters->FindParam(&analysis_type,AnalysisTypeEnum); … … 132 131 nflipped=1; //bootstrap 133 132 while(nflipped){ 134 133 135 134 /*Vector of size number of elements*/ 136 135 vec_elements_neighboring_floatingice=new Vector<IssmDouble>(elements->NumberOfElements(),true); … … 155 154 } 156 155 vec_nodes_on_floatingice->Assemble(); 157 156 158 157 #ifdef _HAVE_MPI_ 159 158 MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,IssmComm::GetComm()); -
issm/trunk-jpl/src/c/modules/HoleFillerx/HoleFillerx.cpp
r12520 r13622 10 10 11 11 int HoleFillerx(double** pimageout,double* image, int lines,int samps,int smooth){ 12 12 13 13 FILE *fp1; 14 14 unsigned long filesize; … … 29 29 float ssw, wsw, wnw, nnw; 30 30 float sum; 31 31 32 32 time_t t1, t2; 33 33 … … 45 45 image2 = xNew<double>(lines*samps); 46 46 memcpy(image2,image,lines*samps*sizeof(double)); 47 47 48 48 for ( i = 0; i < lines; i++ ){ 49 49 for ( j = 0; j < samps; j++ ){ … … 84 84 85 85 again2: 86 86 87 87 #ifdef _DEBUG2_ 88 88 counter=0; … … 99 99 fflush( stdout ); 100 100 #endif 101 101 102 102 afterfirst2: 103 103 … … 128 128 129 129 /* For void edge pixels: */ 130 130 131 131 nsteps = 0.0; ssteps = 0.0; esteps = 0.0; wsteps = 0.0; 132 132 nwsteps = 0.0; nesteps = 0.0; swsteps = 0.0; sesteps = 0.0; … … 139 139 ssw = 0.0; wsw = 0.0; wnw = 0.0; nnw = 0.0; 140 140 141 142 141 /** NSEW **/ 143 142 for ( ii = i - 1; ii >= 0; ii-- ){ /* North */ … … 150 149 if ( *(image2+ii*samps+j) != 0 ){ south = *(image2+ii*samps+j); ssteps = ii - i; break;} 151 150 } 152 151 153 152 for ( jj = j - 1; jj >= 0; jj-- ){ /* West */ 154 153 if ( jj <= 0 ) { west = 0; wsteps = 0; break;} 155 154 if ( *(image2+i*samps+jj) != 0 ){ west = *(image2+i*samps+jj); wsteps = j - jj; break;} 156 155 } 157 156 158 157 for ( jj = j + 1; jj < samps; jj++ ){ /* East */ 159 158 if ( jj >= samps-1 ){ east = 0; esteps = 0; break;} … … 161 160 } 162 161 163 164 162 /** Diagonals **/ 165 163 /* Southeast */ … … 169 167 if ( *(image2+ii*samps+jj) != 0 ){ se = *(image2+ii*samps+jj); sesteps = 1.4142 * k; break;} 170 168 } 171 169 172 170 /* Northeast */ 173 171 for ( k = 1; k < infinit; k++ ){ … … 176 174 if ( *(image2+ii*samps+jj) != 0 ){ ne = *(image2+ii*samps+jj); nesteps = 1.4142 * k; break;} 177 175 } 178 176 179 177 /* Northwest */ 180 178 for ( k = 1; k < infinit; k++ ){ … … 183 181 if ( *(image2+ii*samps+jj) != 0 ){ nw = *(image2+ii*samps+jj); nwsteps = 1.4142 * k; break;} 184 182 } 185 183 186 184 /* Southwest */ 187 185 for ( k = 1; k < infinit; k++ ){ … … 198 196 if ( *(image2+ii*samps+jj) != 0 ){ nne = *(image2+ii*samps+jj); nnesteps = 2.2361 * k; break;} 199 197 } 200 198 201 199 /* ENE */ 202 200 for ( k = 1; k < infinit; k++ ){ … … 205 203 if ( *(image2+ii*samps+jj) != 0 ){ ene = *(image2+ii*samps+jj); enesteps = 2.2361 * k; break;} 206 204 } 207 205 208 206 /* ESE */ 209 207 for ( k = 1; k < infinit; k++ ){ … … 212 210 if ( *(image2+ii*samps+jj) != 0 ){ ese = *(image2+ii*samps+jj); esesteps = 2.2361 * k; break;} 213 211 } 214 212 215 213 /* SSE */ 216 214 for ( k = 1; k < infinit; k++ ){ … … 226 224 if ( *(image2+ii*samps+jj) != 0 ){ ssw = *(image2+ii*samps+jj); sswsteps = 2.2361 * k; break;} 227 225 } 228 226 229 227 /* WSW */ 230 228 for ( k = 1; k < infinit; k++ ){ … … 233 231 if ( *(image2+ii*samps+jj) != 0 ){ wsw = *(image2+ii*samps+jj); wswsteps = 2.2361 * k; break;} 234 232 } 235 233 236 234 /* WNW */ 237 235 for ( k = 1; k < infinit; k++ ){ … … 240 238 if ( *(image2+ii*samps+jj) != 0 ){ wnw = *(image2+ii*samps+jj); wnwsteps = 2.2361 * k; break;} 241 239 } 242 240 243 241 /* NNW */ 244 242 for ( k = 1; k < infinit; k++ ){ … … 247 245 if ( *(image2+ii*samps+jj) != 0 ){ nnw = *(image2+ii*samps+jj); nnwsteps = 2.2361 * k; break;} 248 246 } 249 247 250 248 elev = 0; range = 0; 251 249 /*NSEW*/ … … 259 257 if ( swsteps > 0.5 ){ elev += sw / swsteps; range += 1.00 / swsteps;} 260 258 if ( sesteps > 0.5 ){ elev += se / sesteps; range += 1.00 / sesteps;} 261 259 262 260 /*Other 8*/ 263 261 if ( nnesteps > 0.5 ){ elev += nne / nnesteps; range += 1.00 / nnesteps;} … … 269 267 if ( wnwsteps > 0.5 ){ elev += wnw / wnwsteps; range += 1.00 / wnwsteps;} 270 268 if ( nnwsteps > 0.5 ){ elev += nnw / nnwsteps; range += 1.00 / nnwsteps;} 271 269 272 270 //temp = ( elev / range ) + 0.5 ; 273 271 temp = ( elev / range ); … … 275 273 //if ( temp > 10000 ) temp = 10000; 276 274 //if ( temp < 0 ) temp = 0; 277 275 278 276 #ifdef _DEBUG2_ 279 277 //_printLine_(temp << " " << elev << " " << range << " "); … … 284 282 } 285 283 286 287 288 284 for ( i = 0; i < lines; i++ ){ 289 285 for ( j = 0; j < samps; j++ ){ … … 292 288 } 293 289 294 295 290 for ( i = 0; i < lines; i++ ){ 296 291 for ( j = 0; j < samps; j++ ){ … … 298 293 } 299 294 } 300 295 301 296 if ( smooth == 0 ) goto there2; 302 303 297 304 298 /************************ SMOOTH THE RESULT ***********************/ 305 299 306 300 image4 = xNew<double>(lines*samps); 307 301 memcpy(image4,image3,lines*samps*sizeof(double)); 308 309 302 310 303 for ( i = 0; i < lines; i++ ) { 311 304 for ( j = 0; j < samps; j++ ) { 312 305 if ( *(image4+i*samps+j) != 0 ) { *(image3+i*samps+j) = *(image2+i*samps+j) ; continue; } 313 314 306 315 307 for ( k = 1; k < infinit; k++ ) { /* Find the smallest box size with data */ … … 324 316 } 325 317 } 326 318 327 319 k_nowset: 328 320 k = k / 4; if ( k < 1 ) k = 1; /* Errrrr. Make it fourth size */ … … 343 335 } 344 336 345 346 337 there2: 347 348 338 349 339 /*Allocate output image: */ … … 352 342 353 343 time(&t2); 354 344 355 345 #ifdef _DEBUG2_ 356 346 _printString_( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b" ); … … 358 348 printf ( "\n"); 359 349 #endif 360 350 361 351 end: 362 352 -
issm/trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp
r4967 r13622 10 10 11 11 void InputDuplicatex(Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters,int original_enum, int new_enum){ 12 12 13 13 /*Go through elemnets, and ask to reinitialie the input: */ 14 14 int i; -
issm/trunk-jpl/src/c/modules/InputUpdateFromConstantx/InputUpdateFromConstantx.cpp
r12470 r13622 23 23 load->InputUpdateFromConstant(constant,name); 24 24 } 25 25 26 26 for(i=0;i<materials->Size();i++){ 27 27 Material* material=(Material*)materials->GetObjectByOffset(i); … … 30 30 } 31 31 void InputUpdateFromConstantx( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int constant, int name){ 32 32 33 33 int i; 34 34 … … 43 43 load->InputUpdateFromConstant(constant,name); 44 44 } 45 45 46 46 for(i=0;i<materials->Size();i++){ 47 47 Material* material=(Material*)materials->GetObjectByOffset(i); … … 63 63 load->InputUpdateFromConstant(constant,name); 64 64 } 65 65 66 66 for(i=0;i<materials->Size();i++){ 67 67 Material* material=(Material*)materials->GetObjectByOffset(i); -
issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp
r13590 r13622 14 14 int i,j,k,l; 15 15 int dummy; 16 16 17 17 int numberofvertices; 18 18 int nrows; … … 35 35 36 36 for(i=0;i<numvariables;i++){ 37 37 38 38 descriptor=variables_descriptors[i]; 39 39 40 40 /*From descriptor, figure out if the variable is scaled, indexed, nodal, or just a simple variable: */ 41 41 if (strncmp(descriptor,"scaled_",7)==0){ 42 42 43 43 /*Variable is scaled. Determine root name of variable (ex: scaled_DragCoefficient_1 -> DragCoefficient). Allocate distributed_values and fill the 44 44 * distributed_values with the next npart variables: */ 45 45 46 46 //strcpy(root,strstr(descriptor,"_")+1); *strstr(root,"_")='\0'; 47 47 memcpy(root,strstr(descriptor,"_")+1,(strlen(strstr(descriptor,"_")+1)+1)*sizeof(char)); 48 48 *strstr(root,"_")='\0'; 49 50 49 51 50 distributed_values=xNew<double>(npart); … … 82 81 PetscSynchronizedFlush(IssmComm::GetComm()); 83 82 #endif 84 85 83 86 84 /*Update inputs using the parameter matrix: */ -
issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp
r13216 r13622 23 23 } 24 24 25 26 25 void InputUpdateFromSolutionx( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* solution){ 27 26 -
issm/trunk-jpl/src/c/modules/InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.cpp
r13216 r13622 21 21 } 22 22 23 24 23 void InputUpdateFromVectorDakotax( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,double* vector, int name, int type){ 25 24 … … 46 45 47 46 void InputUpdateFromVectorDakotax( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int* vector, int name, int type){ 48 47 49 48 int i; 50 49 -
issm/trunk-jpl/src/c/modules/InputUpdateFromVectorx/InputUpdateFromVectorx.cpp
r13216 r13622 21 21 } 22 22 23 24 23 void InputUpdateFromVectorx( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* vector, int name, int type){ 25 24 … … 46 45 47 46 void InputUpdateFromVectorx( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int* vector, int name, int type){ 48 47 49 48 int i; 50 49 -
issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp
r13229 r13622 21 21 /*output: */ 22 22 SeqVec<IssmPDouble>* data_mesh=NULL; 23 23 24 24 /*Intermediary*/ 25 25 double* x=NULL; … … 312 312 _assert_(x2>x1 && y2>y1); 313 313 _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1); 314 314 315 315 double xm=(x2-x1)/2; 316 316 double ym=(y2-y1)/2; -
issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp
r13229 r13622 14 14 double* x_prime, double* y_prime, int nods_prime, 15 15 double* default_values,int num_default_values,Contour<IssmPDouble>** contours,int numcontours){ 16 16 17 17 /*Output*/ 18 18 SeqVec<IssmPDouble>* data_prime=NULL; -
issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp
r13229 r13622 18 18 int my_thread = handle->id; 19 19 int num_threads = handle->num; 20 20 21 21 /*recover parameters :*/ 22 22 int interpolation_type = gate->interpolation_type; -
issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp
r13056 r13622 69 69 } 70 70 } 71 71 72 72 /*Get extreme coordinates of the grid*/ 73 73 if (xflip){ -
issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp
r13056 r13622 15 15 int InterpFromMeshToMesh2dx(double** pdata_interp,int* index_data,double* x_data,double* y_data,int nods_data,int nels_data, 16 16 double* data,int M_data,int N_data,double* x_interp,double* y_interp,int N_interp,Options* options){ 17 17 18 18 /*Output*/ 19 19 double* data_interp=NULL; -
issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp
r13056 r13622 11 11 12 12 void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type){ 13 14 13 15 14 /*intermediary: */ … … 48 47 /*static: just create Constraints objects*/ 49 48 count=0; 50 49 51 50 /*Create Constraints from x,y,z: */ 52 51 for (i=0;i<numberofvertices;i++){ 53 52 54 53 /*keep only this partition's nodes:*/ 55 54 if((iomodel->my_vertices[i])){ … … 78 77 /*Create constraints from x,y,z: */ 79 78 for (i=0;i<numberofvertices;i++){ 80 79 81 80 /*keep only this partition's nodes:*/ 82 81 if((iomodel->my_vertices[i])){ -
issm/trunk-jpl/src/c/modules/KMLFileReadx/KMLFileReadx.cpp
r12520 r13622 65 65 return(kfil); 66 66 } 67 -
issm/trunk-jpl/src/c/modules/KMLOverlayx/KMLOverlayx.cpp
r12519 r13622 79 79 kdoc=NULL; 80 80 81 82 81 /* write kml file */ 83 82 … … 95 94 return; 96 95 } 97 -
issm/trunk-jpl/src/c/modules/Kml2Expx/Kml2Expx.cpp
r13056 r13622 69 69 return(iret); 70 70 } 71 -
issm/trunk-jpl/src/c/modules/MassFluxx/MassFluxx.cpp
r13590 r13622 16 16 int element_id; 17 17 bool ispresent=false; 18 18 19 19 /*output: */ 20 20 IssmDouble mass_flux=0; … … 44 44 segments=array[counter-1]; //matlab to "C" indexing 45 45 num_segments=mdims_array[counter-1]; 46 46 47 47 /*Go through segments, and then elements, and figure out which elements belong to a segment. 48 48 * When we find one, use the element to compute the mass flux on the segment: */ … … 72 72 xDelete<int>(ndims_array); 73 73 xDelete<IssmDouble*>(array); 74 74 75 75 /*Assign output pointers: */ 76 76 *pmass_flux=mass_flux; -
issm/trunk-jpl/src/c/modules/MaxAbsVxx/MaxAbsVxx.cpp
r13590 r13622 10 10 11 11 void MaxAbsVxx( IssmDouble* pmaxabsvx, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ 12 12 13 13 int i; 14 14 IssmDouble maxabsvx; … … 20 20 Element* element=(Element*)elements->GetObjectByOffset(i); 21 21 element->MaxAbsVx(&element_maxabsvx,process_units); //go pick up the maximum velocity in the inputs 22 22 23 23 if(i==0)maxabsvx=element_maxabsvx; //initialize maxabsvx 24 24 else{ -
issm/trunk-jpl/src/c/modules/MaxAbsVyx/MaxAbsVyx.cpp
r13590 r13622 10 10 11 11 void MaxAbsVyx( IssmDouble* pmaxabsvy, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ 12 12 13 13 int i; 14 14 IssmDouble maxabsvy; … … 16 16 IssmDouble element_maxabsvy; 17 17 18 19 18 /*Go through elements, and request velocity: */ 20 19 for(i=0;i<elements->Size();i++){ 21 20 Element* element=(Element*)elements->GetObjectByOffset(i); 22 21 element->MaxAbsVy(&element_maxabsvy,process_units); //go pick up the maximum velocity in the inputs 23 22 24 23 if(i==0)maxabsvy=element_maxabsvy; //initialize maxabsvy 25 24 else{ -
issm/trunk-jpl/src/c/modules/MaxAbsVzx/MaxAbsVzx.cpp
r13590 r13622 10 10 11 11 void MaxAbsVzx( IssmDouble* pmaxabsvz, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ 12 12 13 13 int i; 14 14 IssmDouble maxabsvz; … … 20 20 Element* element=(Element*)elements->GetObjectByOffset(i); 21 21 element->MaxAbsVz(&element_maxabsvz,process_units); //go pick up the minimum velocity in the inputs 22 22 23 23 if(i==0)maxabsvz=element_maxabsvz; //initialize maxabsvz 24 24 else{ -
issm/trunk-jpl/src/c/modules/MaxVelx/MaxVelx.cpp
r13590 r13622 11 11 12 12 void MaxVelx( IssmDouble* pmaxvel, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ 13 13 14 14 int i; 15 15 IssmDouble maxvel; … … 21 21 Element* element=(Element*)elements->GetObjectByOffset(i); 22 22 element->MaxVel(&element_maxvel,process_units); //go pick up the maximum velocity in the inputs 23 23 24 24 if(i==0)maxvel=element_maxvel; //initialize maxvel 25 25 else{ -
issm/trunk-jpl/src/c/modules/MaxVxx/MaxVxx.cpp
r13590 r13622 10 10 11 11 void MaxVxx( IssmDouble* pmaxvx, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ 12 12 13 13 int i; 14 14 IssmDouble maxvx; … … 20 20 Element* element=(Element*)elements->GetObjectByOffset(i); 21 21 element->MaxVx(&element_maxvx,process_units); //go pick up the minimum velocity in the inputs 22 22 23 23 if(i==0)maxvx=element_maxvx; //initialize maxvx 24 24 else{ -
issm/trunk-jpl/src/c/modules/MaxVyx/MaxVyx.cpp
r13590 r13622 20 20 Element* element=(Element*)elements->GetObjectByOffset(i); 21 21 element->MaxVy(&element_maxvy,process_units); //go pick up the minimum velocity in the inputs 22 22 23 23 if(i==0)maxvy=element_maxvy; //initialize maxvy 24 24 else{ -
issm/trunk-jpl/src/c/modules/MaxVzx/MaxVzx.cpp
r13590 r13622 11 11 12 12 void MaxVzx( IssmDouble* pmaxvz, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ 13 13 14 14 int i; 15 15 IssmDouble maxvz; … … 21 21 Element* element=(Element*)elements->GetObjectByOffset(i); 22 22 element->MaxVz(&element_maxvz,process_units); //go pick up the minimum velocity in the inputs 23 23 24 24 if(i==0)maxvz=element_maxvz; //initialize maxvz 25 25 else{ -
issm/trunk-jpl/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp
r13216 r13622 40 40 VecMergex(ug,uf,nodes,parameters,FsetEnum); 41 41 } 42 42 43 43 /*Merge s set back into g set: */ 44 44 if(ssize){ -
issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp
r13220 r13622 3 3 4 4 #include "./MeshProfileIntersectionx.h" 5 5 6 6 void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment){ 7 7 … … 14 14 double beta1,beta2; 15 15 double gamma1,gamma2; 16 16 17 17 int edge1,edge2,edge3; 18 18 … … 21 21 double xfinal[2],yfinal[2]; 22 22 23 24 23 /*edge 1: */ 25 24 xel[0]=xnodes[0]; yel[0]=ynodes[0]; xel[1]=xnodes[1]; yel[1]=ynodes[1]; … … 41 40 } 42 41 else if( ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum)) ){ 43 42 44 43 /*segment interscts 2 opposite edges of our triangle, at 2 segment coordinates, pick up the lowest (coord1) and highest (coord2): */ 45 44 if((edge1==IntersectEnum) && (edge2==IntersectEnum)) {coord1=min(alpha1,beta1); coord2=max(alpha1,beta1);} … … 77 76 coord2=1.0; 78 77 } 79 78 80 79 xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]); 81 80 xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]); -
issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp
r8303 r13622 3 3 4 4 #include "./MeshProfileIntersectionx.h" 5 5 6 6 void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes){ 7 7 -
issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp
r13220 r13622 32 32 /*Loop through all contours: */ 33 33 for (i=0;i<numcontours;i++){ 34 34 35 35 /*retrieve contour info: */ 36 36 contouri=*(contours+i); … … 41 41 /*determine segmentsi and numsegsi for this contour and the mesh intersection: */ 42 42 MeshSegmentsIntersection(&segmentsi,&numsegsi,index,x,y,nel,nods,xc,yc,numnodes); 43 43 44 44 /*save segmentsi: */ 45 45 allsegments[i]=segmentsi; -
issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp
r13220 r13622 12 12 Segment<double>* segment=NULL; 13 13 int numsegs; 14 14 15 15 /*intermediary: */ 16 16 DataSet* segments_dataset=NULL; … … 35 35 for(i=0;i<numsegs;i++){ 36 36 Segment<double>* segment=(Segment<double>*)segments_dataset->GetObjectByOffset(i); 37 37 38 38 /*x1,y1,x2,y2 then element_id: */ 39 39 *(segments+5*i+0)=segment->x1; -
issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/NodeInElement.cpp
r8303 r13622 19 19 y3=ynodes[2]; 20 20 21 22 21 /*compute determinant: */ 23 22 det=x1*y2-x1*y3-x3*y2-x2*y1+x2*y3+x3*y1; 24 23 25 24 /*area coordinates: */ 26 25 lambda1=((y2-y3)*(x-x3)+(x3-x2)*(y-y3))/det; -
issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp
r12155 r13622 11 11 double alpha=-1; 12 12 double beta=-1; 13 13 14 14 double xA,xB,xC,xD,yA,yB,yC,yD; 15 15 double O2A[2],O2B[2],O1C[2],O1D[2]; … … 29 29 O1C[0]=xC -(xA/2+xB/2); O1C[1]=yC -(yA/2+yB/2); 30 30 O1D[0]=xD -(xA/2+xB/2); O1D[1]=yD -(yA/2+yB/2); 31 32 31 33 32 n1[0]=yA-yB; n1[1]=xB-xA; //normal vector to segA -
issm/trunk-jpl/src/c/modules/MinVelx/MinVelx.cpp
r13590 r13622 11 11 12 12 void MinVelx( IssmDouble* pminvel, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ 13 13 14 14 int i; 15 15 IssmDouble minvel; … … 21 21 Element* element=(Element*)elements->GetObjectByOffset(i); 22 22 element->MinVel(&element_minvel,process_units); //go pick up the minimum velocity in the inputs 23 23 24 24 if(i==0)minvel=element_minvel; //initialize minvel 25 25 else{ -
issm/trunk-jpl/src/c/modules/MinVxx/MinVxx.cpp
r13590 r13622 10 10 11 11 void MinVxx( IssmDouble* pminvx, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ 12 12 13 13 int i; 14 14 IssmDouble minvx; … … 20 20 Element* element=(Element*)elements->GetObjectByOffset(i); 21 21 element->MinVx(&element_minvx,process_units); //go pick up the minimum velocity in the inputs 22 22 23 23 if(i==0)minvx=element_minvx; //initialize minvx 24 24 else{ -
issm/trunk-jpl/src/c/modules/MinVyx/MinVyx.cpp
r13590 r13622 10 10 11 11 void MinVyx( IssmDouble* pminvy, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ 12 12 13 13 int i; 14 14 IssmDouble minvy; … … 20 20 Element* element=(Element*)elements->GetObjectByOffset(i); 21 21 element->MinVy(&element_minvy,process_units); //go pick up the minimum velocity in the inputs 22 22 23 23 if(i==0)minvy=element_minvy; //initialize minvy 24 24 else{ -
issm/trunk-jpl/src/c/modules/MinVzx/MinVzx.cpp
r13590 r13622 10 10 11 11 void MinVzx( IssmDouble* pminvz, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units){ 12 12 13 13 int i; 14 14 IssmDouble minvz; … … 20 20 Element* element=(Element*)elements->GetObjectByOffset(i); 21 21 element->MinVz(&element_minvz,process_units); //go pick up the minimum velocity in the inputs 22 22 23 23 if(i==0)minvz=element_minvz; //initialize minvz 24 24 else{ -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp
r13550 r13622 14 14 15 15 void CreateParametersAutodiff(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type){ 16 16 17 17 int i; 18 18 Parameters *parameters = NULL; … … 26 26 int* indices=NULL; 27 27 int num_indices; 28 28 29 29 IssmDouble* xp=NULL; 30 30 IssmDouble* xp_backup=NULL; 31 31 int num_ind,local_num_ind; 32 32 DataSet* dependent_objects=NULL; 33 33 34 34 /*Get parameters: */ 35 35 parameters=*pparameters; … … 37 37 /*retrieve some parameters: */ 38 38 iomodel->Constant(&autodiff_analysis,AutodiffIsautodiffEnum); 39 39 40 40 if(autodiff_analysis){ 41 41 42 42 #ifdef _HAVE_ADOLC_ 43 43 44 44 /*Copy some parameters from IoModel to parameters dataset: */ 45 45 parameters->AddObject(iomodel->CopyConstantObject(AutodiffKeepEnum)); … … 48 48 iomodel->Constant(&autodiff_driver,AutodiffDriverEnum); 49 49 parameters->AddObject(iomodel->CopyConstantObject(AutodiffDriverEnum)); 50 50 51 51 if(strcmp(autodiff_driver,"fos_forward")==0){ 52 52 parameters->AddObject(iomodel->CopyConstantObject(AutodiffFosForwardIndexEnum)); … … 66 66 dependent_objects=new DataSet(); 67 67 num_dep=0; 68 68 69 69 if(num_dependent_objects){ 70 70 iomodel->FetchData(&names,&dummy,&dummy,AutodiffDependentObjectNamesEnum); … … 114 114 #endif 115 115 } 116 116 117 117 #ifdef _HAVE_ADOLC_ 118 118 /*initialize a placeholder to store solver pointers: {{{*/ … … 129 129 /*}}}*/ 130 130 #endif 131 131 132 132 /*Assign output pointer: */ 133 133 *pparameters=parameters; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp
r12832 r13622 11 11 12 12 int stabilization; 13 13 14 14 /*Fetch parameters: */ 15 15 iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp
r13073 r13622 31 31 /*Create loads if they do not exist yet*/ 32 32 if(!loads) loads = new Loads(); 33 33 34 34 /*Loads only in DG*/ 35 35 if (stabilization==3){ -
issm/trunk-jpl/src/c/modules/ModelProcessorx/BedSlope/CreateConstraintsBedSlope.cpp
r12832 r13622 15 15 /*Output*/ 16 16 Constraints* constraints = NULL; 17 17 18 18 /*Recover pointer: */ 19 19 constraints=*pconstraints; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/BedSlope/CreateNodesBedSlope.cpp
r12832 r13622 31 31 /*Create nodes if they do not exist yet*/ 32 32 if(!nodes) nodes = new Nodes(); 33 33 34 34 /*Continuous Galerkin partition of nodes: */ 35 35 NodesPartitioning(&iomodel->my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel,continuous_galerkin); 36 36 37 37 /*First fetch data: */ 38 38 iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); … … 40 40 for (i=0;i<numberofvertices;i++){ 41 41 if(iomodel->my_vertices[i]){ 42 42 43 43 /*Add node to nodes dataset: */ 44 44 nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,BedSlopeAnalysisEnum)); … … 49 49 /*Clean fetched data: */ 50 50 iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); 51 51 52 52 /*Assign output pointer: */ 53 53 *pnodes=nodes; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/BedSlope/UpdateElementsBedSlope.cpp
r12832 r13622 42 42 iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); 43 43 } 44 44 45 45 /*Free data: */ 46 46 iomodel->DeleteData(1,MeshElementsEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp
r13283 r13622 13 13 14 14 void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type){ 15 15 16 16 int i; 17 17 Parameters *parameters = NULL; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp
r13119 r13622 25 25 bool control_analysis; 26 26 27 28 27 /*Fetch parameters: */ 29 28 iomodel->Constant(&numberofelements,MeshNumberofelementsEnum); … … 65 64 } 66 65 } 67 66 68 67 /*Free data: */ 69 68 iomodel->DeleteData(1+4+6,MeshElementsEnum,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum,BalancethicknessThickeningRateEnum,VxEnum,VyEnum,FrictionCoefficientEnum,MaterialsRheologyBEnum,MaterialsRheologyZEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp
r13277 r13622 15 15 #include "./ModelProcessorx.h" 16 16 17 18 17 void CreateDataSets(Elements** pelements,Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,const int analysis_type,const int nummodels,int analysis_counter){ 19 18 … … 22 21 Materials *materials = NULL; 23 22 Parameters *parameters = NULL; 24 23 25 24 /*Create elements, vertices and materials, independent of analysis_type: */ 26 25 CreateElementsVerticesAndMaterials(pelements, pvertices, pmaterials, iomodel,nummodels); … … 40 39 UpdateElementsDiagnosticHoriz(elements,iomodel,analysis_counter,analysis_type); 41 40 break; 42 41 43 42 case DiagnosticVertAnalysisEnum: 44 43 CreateNodesDiagnosticVert(pnodes, iomodel); … … 47 46 UpdateElementsDiagnosticVert(elements,iomodel,analysis_counter,analysis_type); 48 47 break; 49 48 50 49 case DiagnosticHutterAnalysisEnum: 51 50 CreateNodesDiagnosticHutter(pnodes, iomodel); … … 55 54 break; 56 55 #endif 57 56 58 57 #ifdef _HAVE_HYDROLOGY_ 59 58 case HydrologyAnalysisEnum: … … 72 71 UpdateElementsThermal(elements,iomodel,analysis_counter,analysis_type); 73 72 break; 74 73 75 74 case EnthalpyAnalysisEnum: 76 75 CreateNodesEnthalpy(pnodes, iomodel); … … 79 78 UpdateElementsEnthalpy(elements,iomodel,analysis_counter,analysis_type); 80 79 break; 81 80 82 81 case MeltingAnalysisEnum: 83 82 CreateNodesMelting(pnodes, iomodel); … … 122 121 #endif 123 122 124 125 123 default: 126 124 _error_("analysis_type: " << EnumToStringx(analysis_type) << " not supported yet!"); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
r13129 r13622 41 41 vertices = new Vertices(); 42 42 materials = new Materials(); 43 43 44 44 /*First, partition elements and vertices. Nodes will partitioned on a per analysis_type basis. If partitining already done, ignore: */ 45 45 ElementsAndVerticesPartitioning(&iomodel->my_elements,&iomodel->my_vertices,iomodel); 46 46 47 47 iomodel->FetchData(2,MeshElementsEnum,MeshElementconnectivityEnum); 48 48 #ifdef _HAVE_3D_ … … 50 50 #endif 51 51 if(control_analysis)iomodel->FetchData(3,InversionControlParametersEnum,InversionMinParametersEnum,InversionMaxParametersEnum); 52 52 53 53 /*Create elements*/ 54 54 for (i=0;i<numberofelements;i++){ … … 62 62 } 63 63 } 64 64 65 65 /*Create materials*/ 66 66 switch(materials_type){ … … 84 84 /*Add new constant material property to materials, at the end: */ 85 85 materials->AddObject(new Matpar(numberofelements+1,iomodel));//put it at the end of the materials 86 86 87 87 /*Create vertices: */ 88 88 … … 90 90 iomodel->FetchData(6,MeshElementsEnum,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum); 91 91 CreateNumberNodeToElementConnectivity(iomodel); 92 92 93 93 for (i=0;i<numberofvertices;i++){ 94 94 95 95 /*vertices and nodes (same number, as we are running continuous galerkin formulation): */ 96 96 if(iomodel->my_vertices[i]){ 97 97 98 98 /*Add vertex to vertices dataset: */ 99 99 vertices->AddObject(new Vertex(i+1,i,i,iomodel)); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
r13486 r13622 8 8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 9 9 #endif 10 11 10 12 11 #include "../../Container/Container.h" … … 20 19 21 20 void CreateParameters(Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,int analysis_type,int analysis_counter){ 22 21 23 22 int i,j,m,k; 24 23 int numoutputs; … … 45 44 /*Initialize dataset: */ 46 45 parameters = new Parameters(); 47 46 48 47 /*Copy some constants from iomodel */ 49 48 parameters->AddObject(iomodel->CopyConstantObject(MeshDimensionEnum)); … … 151 150 if(numoutputs)parameters->AddObject(new IntVecParam(SteadystateRequestedOutputsEnum,requestedoutputs,numoutputs)); 152 151 iomodel->DeleteData(requestedoutputs,SteadystateRequestedOutputsEnum); 153 152 154 153 iomodel->FetchData(&requestedoutputs,&numoutputs,NULL,PrognosticRequestedOutputsEnum); 155 154 parameters->AddObject(new IntParam(PrognosticNumRequestedOutputsEnum,numoutputs)); 156 155 if(numoutputs)parameters->AddObject(new IntVecParam(PrognosticRequestedOutputsEnum,requestedoutputs,numoutputs)); 157 156 iomodel->DeleteData(requestedoutputs,PrognosticRequestedOutputsEnum); 158 159 157 160 158 /*Deal with mass flux segments: {{{*/ 161 159 iomodel->FetchData(&qmu_mass_flux_present,QmuMassFluxSegmentsPresentEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp
r13483 r13622 17 17 /*variable declarations: {{{*/ 18 18 int i,j,k; 19 19 20 20 Parameters* parameters = NULL; 21 21 int second_count; 22 22 23 23 int* part=NULL; 24 24 double* dpart=NULL; … … 35 35 char* qmuerrname=NULL; 36 36 char* qmuoutname=NULL; 37 37 38 38 //descriptors: 39 39 char tag[50]; 40 41 40 42 41 int M; 43 42 int m,n; … … 50 49 51 50 /*}}}*/ 52 51 53 52 /*recover parameters : */ 54 53 parameters=*pparameters; … … 92 91 /*Deal with partitioning: {{{*/ 93 92 /*partition vertices in iomodel->qmu_npart parts, unless a partition is already present: */ 94 93 95 94 parameters->AddObject(iomodel->CopyConstantObject(QmuNumberofpartitionsEnum)); 96 95 iomodel->FetchData(&dpart,NULL,NULL,QmuPartitionEnum); … … 107 106 /*}}}*/ 108 107 /*Deal with data needed because of qmu variables: {{{*/ 109 108 110 109 for(i=0;i<numvariabledescriptors;i++){ 111 110 … … 113 112 /*Ok, we are dealing with a variable that is distributed over nodes. Recover the name of the variable (ex: scaled_Thickness): */ 114 113 sscanf(variabledescriptors[i],"scaled_%s",tag); 115 114 116 115 /*Recover data: */ 117 116 iomodel->FetchData(&dakota_parameter,&nrows,&ncols,StringToEnumx(tag)); … … 124 123 parameters->AddObject(new DoubleTransientMatParam(StringToEnumx(tag),dakota_parameter,nrows,ncols)); 125 124 } 126 125 127 126 /*Free ressources:*/ 128 127 xDelete<double>(dakota_parameter); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
r13283 r13622 64 64 /*Create constraints if they do not exist yet*/ 65 65 if(!constraints) constraints = new Constraints(); 66 66 67 67 /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ 68 68 if(!ismacayealpattyn & !isstokes & !isl1l2){ … … 70 70 return; 71 71 } 72 72 73 73 /*Constraints: fetch data: */ 74 74 iomodel->FetchData(&spcvx,&Mx,&Nx,DiagnosticSpcvxEnum); … … 324 324 } 325 325 } 326 326 327 327 /*Free data: */ 328 328 iomodel->DeleteData(spcvx,DiagnosticSpcvxEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp
r13283 r13622 58 58 return; 59 59 } 60 60 61 61 /*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the nodes 62 62 * referenced by a certain load must belong to the cluster node): */ … … 70 70 /*First load data:*/ 71 71 for (i=0;i<numberofpressureloads;i++){ 72 72 73 73 /*Retrieve element to which this icefront belongs: */ 74 74 if (dim==2) segment_width=4; … … 78 78 /*Now, if this element is not in the partition, pass: */ 79 79 if(!iomodel->my_elements[element]) continue; 80 80 81 81 /*Do not create ice front if Hutter or Stokes elements*/ 82 82 if (reCast<int,IssmDouble>(*(elements_type+element))==HutterApproximationEnum) continue; … … 132 132 iomodel->DeleteData(pressureload,DiagnosticIcefrontEnum); 133 133 134 135 134 /*Create Penpair for penalties: */ 136 135 iomodel->FetchData(&penalties,&numpenalties,NULL,DiagnosticVertexPairingEnum); 137 136 138 137 for(i=0;i<numpenalties;i++){ 139 138 -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp
r13051 r13622 35 35 /*First create nodes*/ 36 36 if(!nodes) nodes = new Nodes(); 37 37 38 38 /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ 39 39 if(!ismacayealpattyn & !isstokes & !isl1l2){ … … 48 48 iomodel->FetchData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBordermacayealEnum,FlowequationBorderstokesEnum, 49 49 MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum); 50 50 51 51 for (i=0;i<numberofvertices;i++){ 52 52 53 53 if(iomodel->my_vertices[i]){ 54 54 55 55 /*Add node to nodes dataset: */ 56 56 nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,DiagnosticHorizAnalysisEnum)); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp
r13129 r13622 85 85 elements->InputDuplicate(VxEnum,InversionVxObsEnum); 86 86 if(dakota_analysis)elements->InputDuplicate(VxEnum,QmuVxEnum); 87 87 88 88 elements->InputDuplicate(VyEnum,VyPicardEnum); 89 89 elements->InputDuplicate(VyEnum,InversionVyObsEnum); 90 90 if(dakota_analysis)elements->InputDuplicate(VyEnum,QmuVyEnum); 91 91 92 92 if(dim==3){ 93 93 elements->InputDuplicate(VzEnum,VzPicardEnum); … … 95 95 if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum); 96 96 } 97 97 98 98 /*Free data: */ 99 99 iomodel->DeleteData(2,MeshElementsEnum,FlowequationElementEquationEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp
r13073 r13622 75 75 /*Free data: */ 76 76 iomodel->DeleteData(3,DiagnosticSpcvxEnum,DiagnosticSpcvyEnum,FlowequationVertexEquationEnum); 77 77 78 78 /*Assign output pointer: */ 79 79 *pconstraints=constraints; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp
r12832 r13622 21 21 /*Create loads if they do not exist yet*/ 22 22 if(!loads) loads = new Loads(); 23 23 24 24 /*Assign output pointer: */ 25 25 *ploads=loads; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHutter/UpdateElementsDiagnosticHutter.cpp
r12832 r13622 19 19 bool ishutter; 20 20 21 22 21 /*Fetch data needed: */ 23 22 iomodel->Constant(&numberofelements,MeshNumberofelementsEnum); … … 38 37 } 39 38 } 40 39 41 40 iomodel->FetchDataToInput(elements,ThicknessEnum); 42 41 iomodel->FetchDataToInput(elements,GeometryHydrostaticRatioEnum); 43 42 44 43 /*Free data: */ 45 44 iomodel->DeleteData(2,MeshElementsEnum,FlowequationElementEquationEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticVert/CreateNodesDiagnosticVert.cpp
r12832 r13622 42 42 /*Continuous Galerkin partition of nodes: */ 43 43 NodesPartitioning(&iomodel->my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel,continuous_galerkin); 44 44 45 45 /*First fetch data: */ 46 46 iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); … … 57 57 /*Clean fetched data: */ 58 58 iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); 59 59 60 60 /*Assign output pointer: */ 61 61 *pnodes=nodes; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticVert/UpdateElementsDiagnosticVert.cpp
r12832 r13622 54 54 /*Free data: */ 55 55 iomodel->DeleteData(1,MeshElementsEnum); 56 56 57 57 } -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp
r13056 r13622 6 6 #include "../../include/include.h" 7 7 #include "../../EnumDefinitions/EnumDefinitions.h" 8 8 9 9 void DistributeNumDofs(DofIndexing* index,int analysis_type,IssmDouble* vertices_type){ 10 10 -
issm/trunk-jpl/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp
r13612 r13622 106 106 107 107 my_elements[i]=true; 108 108 109 109 /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use 110 110 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) … … 114 114 my_vertices[reCast<int>(*(elements+elements_width*i+1))-1]=1; 115 115 my_vertices[reCast<int>(*(elements+elements_width*i+2))-1]=1; 116 116 117 117 if(elements_width==6){ 118 118 my_vertices[reCast<int>(*(elements+elements_width*i+3))-1]=1; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp
r13283 r13622 23 23 IssmDouble heatcapacity; 24 24 IssmDouble referencetemperature; 25 25 26 26 /*Output*/ 27 27 IssmDouble *spcvector = NULL; … … 85 85 /*Create constraints from x,y,z: */ 86 86 for (i=0;i<numberofvertices;i++){ 87 87 88 88 /*keep only this partition's nodes:*/ 89 89 if((iomodel->my_vertices[i])){ -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/CreateLoadsEnthalpy.cpp
r12832 r13622 22 22 if(!loads) loads = new Loads(); 23 23 24 25 24 /*Assign output pointer: */ 26 25 *ploads=loads; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/CreateNodesEnthalpy.cpp
r12832 r13622 25 25 /*Fetch parameters: */ 26 26 iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); 27 27 28 28 /*Recover pointer: */ 29 29 nodes=*pnodes; … … 40 40 for (i=0;i<numberofvertices;i++){ 41 41 if(iomodel->my_vertices[i]){ 42 42 43 43 /*Add node to nodes dataset: */ 44 44 nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,EnthalpyAnalysisEnum)); … … 49 49 /*Clean fetched data: */ 50 50 iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); 51 51 52 52 /*Assign output pointer: */ 53 53 *pnodes=nodes; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp
r12832 r13622 59 59 iomodel->FetchDataToInput(elements,VyEnum); 60 60 iomodel->FetchDataToInput(elements,VzEnum); 61 61 62 62 /*Free data: */ 63 63 iomodel->DeleteData(4,MeshElementsEnum,TemperatureEnum,WaterfractionEnum,PressureEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp
r12832 r13622 23 23 if(!constraints) constraints = new Constraints(); 24 24 IoModelToConstraintsx(constraints,iomodel,HydrologySpcwatercolumnEnum,HydrologyAnalysisEnum); 25 25 26 26 /*Assign output pointer: */ 27 27 *pconstraints=constraints; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp
r12832 r13622 25 25 /*Fetch parameters: */ 26 26 iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); 27 27 28 28 /*Recover pointer: */ 29 29 nodes=*pnodes; … … 46 46 } 47 47 iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); 48 48 49 49 /*Assign output pointer: */ 50 50 *pnodes=nodes; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp
r12832 r13622 17 17 18 18 int numberofelements; 19 19 20 20 /*Fetch data needed: */ 21 21 iomodel->Constant(&numberofelements,MeshNumberofelementsEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Melting/CreateConstraintsMelting.cpp
r12832 r13622 16 16 int i; 17 17 int count; 18 18 19 19 /*Intermediary*/ 20 20 Constraints* constraints = NULL; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Melting/CreateNodesMelting.cpp
r12832 r13622 22 22 /*DataSets: */ 23 23 Nodes* nodes = NULL; 24 24 25 25 /*Fetch parameters: */ 26 26 iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); … … 40 40 41 41 if(iomodel->my_vertices[i]){ 42 42 43 43 /*Add node to nodes dataset: */ 44 44 nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,MeltingAnalysisEnum)); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Melting/UpdateElementsMelting.cpp
r12832 r13622 55 55 iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); 56 56 iomodel->FetchDataToInput(elements,PressureEnum); 57 57 58 58 /*Free data: */ 59 59 iomodel->DeleteData(1,MeshElementsEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp
r13432 r13622 21 21 int i,analysis_type,dim,verbose; 22 22 bool isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy; 23 23 24 24 /*output: */ 25 25 Elements *elements = NULL; … … 31 31 Parameters *parameters = NULL; 32 32 33 34 33 /*Initialize IoModel from input file*/ 35 34 IoModel* iomodel = new IoModel(IOMODEL); … … 67 66 if(solution_type==SteadystateSolutionEnum && analysis_type==MeltingAnalysisEnum && isenthalpy==true) continue; 68 67 if(solution_type==SteadystateSolutionEnum && analysis_type==EnthalpyAnalysisEnum && isenthalpy==false) continue; 69 68 70 69 if(VerboseMProcessor()) _pprintLine_(" creating datasets for analysis " << EnumToStringx(analysis_type)); 71 70 CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,¶meters,iomodel,rootpath,solution_type,analysis_type,nummodels,i); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp
r13056 r13622 22 22 23 23 void NodesPartitioning(bool** pmy_nodes,bool* my_elements, int* my_vertices, IoModel* iomodel, bool continuous){ 24 24 25 25 /*First thing, this is a new partition for a new analysis_type, therefore, to avoid a leak, erase the nodes partition that might come through pmy_nodes: */ 26 26 xDelete<bool>(*pmy_nodes); … … 51 51 } 52 52 53 54 53 void DiscontinuousGalerkinNodesPartitioning(bool** pmy_nodes,bool* my_elements, int* my_vertices, IoModel* iomodel){ 55 54 … … 62 61 * the nodes and the vertices. The vertices are similar to continuous galerkin, but the nodes partitioning involves edges, which mess up sorting of 63 62 * ids. */ 64 63 65 64 int i,j; 66 65 int dim; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp
r12832 r13622 11 11 12 12 int stabilization; 13 13 14 14 /*Fetch parameters: */ 15 15 iomodel->Constant(&stabilization,PrognosticStabilizationEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp
r13283 r13622 25 25 /*DataSet*/ 26 26 Loads* loads = NULL; 27 27 28 28 /*Fetch parameters: */ 29 29 iomodel->Constant(&stabilization,PrognosticStabilizationEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/SortDataSets.cpp
r12832 r13622 14 14 #include "../../EnumDefinitions/EnumDefinitions.h" 15 15 #include "./ModelProcessorx.h" 16 17 16 18 17 void SortDataSets(Elements** pelements,Nodes** pnodes,Vertices** pvertices, Loads** ploads, Materials** pmaterials, Constraints** pconstraints, Parameters** pparameters){ -
issm/trunk-jpl/src/c/modules/ModelProcessorx/SurfaceSlope/CreateConstraintsSurfaceSlope.cpp
r12832 r13622 15 15 /*Output*/ 16 16 Constraints* constraints = NULL; 17 17 18 18 /*Recover pointer: */ 19 19 constraints=*pconstraints; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/SurfaceSlope/CreateNodesSurfaceSlope.cpp
r12832 r13622 31 31 /*Create nodes if they do not exist yet*/ 32 32 if(!nodes) nodes = new Nodes(); 33 33 34 34 /*Continuous Galerkin partition of nodes: */ 35 35 NodesPartitioning(&iomodel->my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel,continuous_galerkin); 36 36 37 37 /*First fetch data: */ 38 38 iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); … … 41 41 42 42 if(iomodel->my_vertices[i]){ 43 43 44 44 /*Add node to nodes dataset: */ 45 45 nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,SurfaceSlopeAnalysisEnum)); … … 50 50 /*Clean fetched data: */ 51 51 iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); 52 52 53 53 /*Assign output pointer: */ 54 54 *pnodes=nodes; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/SurfaceSlope/UpdateElementsSurfaceSlope.cpp
r12832 r13622 37 37 iomodel->FetchDataToInput(elements,BedEnum); 38 38 iomodel->FetchDataToInput(elements,MaskElementonwaterEnum); 39 39 40 40 if (dim==3){ 41 41 iomodel->FetchDataToInput(elements,MeshElementonbedEnum); 42 42 iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); 43 43 } 44 44 45 45 /*Free data: */ 46 46 iomodel->DeleteData(1,MeshElementsEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp
r12832 r13622 18 18 int count; 19 19 int dim; 20 20 21 21 /*Output*/ 22 22 Constraints* constraints = NULL; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp
r13056 r13622 40 40 41 41 for (i=0;i<numberofvertices;i++){ 42 42 43 43 /*keep only this partition's nodes:*/ 44 44 if((iomodel->my_vertices[i]==1)){ -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateNodesThermal.cpp
r12832 r13622 49 49 /*Clean fetched data: */ 50 50 iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); 51 51 52 52 /*Assign output pointer: */ 53 53 *pnodes=nodes; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/UpdateElementsThermal.cpp
r12832 r13622 60 60 iomodel->FetchDataToInput(elements,VyEnum); 61 61 iomodel->FetchDataToInput(elements,VzEnum); 62 62 63 63 if(dakota_analysis){ 64 64 elements->InputDuplicate(TemperatureEnum,QmuTemperatureEnum); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/UpdateCounters.cpp
r12832 r13622 26 26 constraints=*pconstraints; 27 27 28 29 28 if(nodes) iomodel->nodecounter=nodes->NumberOfNodes(); 30 29 else iomodel->nodecounter=0; … … 32 31 if(loads)iomodel->loadcounter=loads->NumberOfLoads(); 33 32 else iomodel->loadcounter=0; 34 33 35 34 if(constraints)iomodel->constraintcounter=constraints->NumberOfConstraints(); 36 35 else iomodel->constraintcounter=0; -
issm/trunk-jpl/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp
r13056 r13622 34 34 double* connectivity=NULL; 35 35 36 37 38 36 /*Allocate connectivity: */ 39 37 connectivity=xNewZeroInit<double>(nods*width); … … 46 44 47 45 for(i=0;i<3;i++){ 48 46 49 47 node=(int)*(elements+n*3+i); //already matlab indexed, elements comes directly from the workspace. 50 48 index=node-1; 51 49 52 50 num_elements=(int)*(connectivity+width*index+maxels); //retrieve number of elements already plugged into the connectivity of this node. 53 51 54 52 already_plugged=0; 55 53 for(j=0;j<num_elements;j++){ … … 64 62 *(connectivity+width*index+num_elements)=element; 65 63 *(connectivity+width*index+maxels)=(double)(num_elements+1); 66 64 67 65 } 68 66 } -
issm/trunk-jpl/src/c/modules/NodesDofx/NodesDofx.cpp
r8816 r13622 15 15 int found=0; 16 16 int i; 17 17 18 18 /*Do we have any nodes for this analysis type? :*/ 19 19 if(nodes->NumberOfNodes(configuration_type)){ -
issm/trunk-jpl/src/c/modules/Orthx/Orthx.cpp
r13216 r13622 6 6 7 7 void Orthx( Vector<IssmDouble>** pnewgradj, Vector<IssmDouble>* gradj, Vector<IssmDouble>* oldgradj){ 8 8 9 9 /*output: */ 10 10 Vector<IssmDouble>* newgradj=NULL; -
issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp
r13612 r13622 15 15 #include "../../io/io.h" 16 16 #include "../../classes/objects/objects.h" 17 17 18 18 void OutputResultsx(Elements* elements, Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Results* results){ 19 19 … … 26 26 char* solutiontypestring = NULL; 27 27 bool dakota_analysis = false; 28 28 29 29 /*retrieve parameters: */ 30 30 parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); 31 31 32 32 /*recover my_rank:*/ 33 33 my_rank=IssmComm::GetRank(); … … 54 54 /*Now, open file for writing, if not already done: */ 55 55 if(!parameters->Exist(OutputFilePointerEnum)){ 56 56 57 57 /*We don't have a file pointer. Retrieve the output file name and open it for writing:*/ 58 58 parameters->FindParam(&outputfilename,OutputfilenameEnum); … … 72 72 } 73 73 xDelete<char>(outputfilename); 74 74 75 75 /*Add file pointer in parameters for further calls to OutputResultsx: */ 76 76 parameters->SetParam(fid,OutputFilePointerEnum); -
issm/trunk-jpl/src/c/modules/ParsePetscOptionsx/ParsePetscOptionsx.cpp
r13612 r13622 59 59 numanalyses=0; 60 60 while ( fgets(line, sizeof line, fid) ){ 61 61 62 62 /*skip comments and empty lines: */ 63 63 if ((line[0]=='%') || (line[0]=='\n') || (line[0]==' ') || (line[0]=='\t') || (line[0]=='\r'))continue; 64 64 65 65 /*Get rid of end of line: */ 66 66 line[strlen(line)-1]='\0'; 67 67 68 68 if (line[0]=='+'){ /*this is the analysis line: */ 69 69 analyses[numanalyses]=StringToEnumx(&line[1]); //skip the '+' -
issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp
r13229 r13622 24 24 my_thread=handle->id; 25 25 num_threads=handle->num; 26 26 27 27 /*recover parameters :*/ 28 28 x=gate->x; … … 57 57 if (j==i)continue; 58 58 distance=sqrt(pow(x[i]-x[j],2)+ pow(y[i]-y[j],2)); 59 59 60 60 if(distance<=mindistance){ 61 61 … … 74 74 /*Free ressources:*/ 75 75 xDelete<bool>(already); 76 76 77 77 return NULL; 78 78 } -
issm/trunk-jpl/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp
r12704 r13622 35 35 IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) 36 36 IssmDouble tstar; // monthly mean surface temp 37 37 38 38 IssmDouble* pdds=NULL; 39 39 IssmDouble* pds=NULL; 40 40 Element* element = NULL; 41 41 42 42 pdds=xNew<IssmDouble>(NPDMAX+1); 43 43 pds=xNew<IssmDouble>(NPDCMAX+1); 44 44 45 45 // initialize PDD (creation of a lookup table) 46 46 tstep = 0.1; … … 55 55 56 56 itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5)); 57 57 58 58 if (itm >= NPDMAX){ 59 59 _printLine_("increase NPDMAX in massBalance.cpp"); … … 73 73 } 74 74 pdds[itm+1] = siglim + DT; 75 75 76 76 //*********compute PD(T) : snow/precip fraction. precipitation falling as snow 77 77 tstepc = 0.1; … … 100 100 pds[itm+1] = 0.; 101 101 // *******END initialize PDD 102 102 103 103 for(i=0;i<elements->Size();i++){ 104 104 element=(Element*)elements->GetObjectByOffset(i); … … 108 108 xDelete<IssmDouble>(pdds); 109 109 xDelete<IssmDouble>(pds); 110 110 111 111 } -
issm/trunk-jpl/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp
r3913 r13622 15 15 RecursivePropagation(pool, connectivity,index, flags); 16 16 } 17 18 17 19 18 void RecursivePropagation(double* pool, double* connectivity, int index, double* flags){ -
issm/trunk-jpl/src/c/modules/Reduceloadx/Reduceloadx.cpp
r13216 r13622 49 49 } 50 50 51 52 51 /*Free ressources and return*/ 53 52 xdelete(&y_s0); -
issm/trunk-jpl/src/c/modules/Reducevectorgtofx/Reducevectorgtofx.cpp
r13216 r13622 5 5 6 6 #include "./Reducevectorgtofx.h" 7 7 8 8 void Reducevectorgtofx(Vector<IssmDouble>** puf, Vector<IssmDouble>* ug, Nodes* nodes,Parameters* parameters){ 9 9 -
issm/trunk-jpl/src/c/modules/RequestedDependentsx/RequestedDependentsx.cpp
r13483 r13622 10 10 11 11 void RequestedDependentsx(Results* results,Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters){ 12 13 14 12 15 13 int i; -
issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp
r13216 r13622 29 29 30 30 void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vector<IssmDouble>* vector){ 31 31 32 32 IssmDouble* serial_vector=NULL; 33 33 -
issm/trunk-jpl/src/c/modules/Shp2Kmlx/Shp2Kmlx.cpp
r13056 r13622 15 15 16 16 #ifdef _HAVE_SHAPELIB_ //only works if Shapelib library has been compiled in. 17 17 18 18 double cm,sp; 19 19 … … 33 33 34 34 #ifdef _HAVE_SHAPELIB_ //only works if Shapelib library has been compiled in. 35 35 36 36 int i,j,k,iret=0; 37 37 int lwidth=1; … … 611 611 #endif 612 612 } 613 -
issm/trunk-jpl/src/c/modules/SmbGradientsx/SmbGradientsx.cpp
r12677 r13622 18 18 19 19 int i; 20 20 21 21 Element* element = NULL; 22 22 23 23 for(i=0;i<elements->Size();i++){ 24 24 element=(Element*)elements->GetObjectByOffset(i); -
issm/trunk-jpl/src/c/modules/SmearFunctionx/SmearFunctionx.cpp
r13599 r13622 11 11 12 12 void SmearFunctionx(Vec* psmearedvector,double (*WeightFunction)(double distance,double radius), int SmearedFieldEnum, double radius,Elements* elements,Nodes* nodes, Parameters* parameters){ 13 13 14 14 Element *element = NULL; 15 15 Vec x=NULL; … … 29 29 gsize=nodes->NumberOfDofs(configuration_type,GsetEnum); 30 30 smearedvector=NewVec(gsize,IssmComm::GetComm()); 31 31 32 32 x=NewVec(gsize,IssmComm::GetComm()); 33 33 y=NewVec(gsize,IssmComm::GetComm()); … … 38 38 } 39 39 40 41 42 40 /*Fill smearedvector vector: */ 43 41 for (i=0;i<elements->Size();i++){ … … 46 44 } 47 45 48 49 46 VecAssemblyBegin(smearedvector); 50 47 VecAssemblyEnd(smearedvector); 51 48 52 49 /*Assign output pointers: */ 53 50 *psmearedvector=smearedvector; 54 51 55 52 } 56 -
issm/trunk-jpl/src/c/modules/Solverx/DofTypesToIndexSet.cpp
r13591 r13622 12 12 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 13 13 #endif 14 14 15 15 void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum){ 16 16 … … 62 62 } 63 63 VecRestoreArray(df,&df_local); 64 64 65 65 /*Create indices sets: */ 66 66 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) -
issm/trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp
r13599 r13622 20 20 Vec uf0_vector = NULL; 21 21 Vec df_vector = NULL; 22 22 23 23 if(uf0) uf0_vector = uf0->vector; 24 24 if(df) df_vector = df->vector; … … 148 148 if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD); 149 149 KSPSolve(ksp,pf,uf); 150 150 151 151 /*Check convergence*/ 152 152 KSPGetIterationNumber(ksp,&iteration_number); … … 155 155 /*Free resources:*/ 156 156 KSPFree(&ksp); 157 157 158 158 /*Assign output pointers:*/ 159 159 *puf=uf; -
issm/trunk-jpl/src/c/modules/SpcNodesx/SpcNodesx.cpp
r9298 r13622 13 13 14 14 for(int i=0;i<constraints->Size();i++){ 15 15 16 16 Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i); 17 17 -
issm/trunk-jpl/src/c/modules/SurfaceAreax/SurfaceAreax.cpp
r13590 r13622 12 12 13 13 void SurfaceAreax( IssmDouble* pS, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ 14 14 15 15 /*Intermediary*/ 16 16 Element* element=NULL; … … 20 20 IssmDouble S=0; 21 21 IssmDouble S_sum; 22 22 23 23 /*Compute gradients: */ 24 24 for (i=0;i<elements->Size();i++){ -
issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
r13216 r13622 11 11 12 12 void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){ 13 13 14 14 /*intermediary: */ 15 15 int i,j; … … 19 19 Element *element = NULL; 20 20 Load *load = NULL; 21 21 22 22 /*output: */ 23 23 Matrix<IssmDouble>* Kff = NULL; … … 70 70 df->Assemble(); 71 71 } 72 72 73 73 if(pflag){ 74 74 … … 106 106 } 107 107 108 109 108 if(penalty_pflag){ 110 109 -
issm/trunk-jpl/src/c/modules/TimeAdaptx/TimeAdaptx.cpp
r13590 r13622 24 24 /*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */ 25 25 element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt(); 26 26 27 27 for (i=1;i<elements->Size();i++){ 28 28 element=(Element*)elements->GetObjectByOffset(i); -
issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp
r13227 r13622 76 76 } 77 77 } 78 78 79 79 /*fill in the point attribute list: */ 80 80 in.pointattributelist = xNew<REAL>(in.numberofpoints*in.numberofpointattributes); 81 81 for (i=0;i<in.numberofpoints;i++) in.pointattributelist[i] = 0.0; 82 82 83 83 /*fill in the point marker list: */ 84 84 in.pointmarkerlist = xNew<int>(in.numberofpoints); … … 96 96 in.numberofsegments+=contour->nods-1; 97 97 } 98 98 99 99 in.segmentlist = xNew<int>(in.numberofsegments*2); 100 100 in.segmentmarkerlist = xNewZeroInit<int>(in.numberofsegments); … … 128 128 counter++; 129 129 } 130 130 131 131 /*Build regions: */ 132 132 in.numberofregions = 0; … … 194 194 index_matrix=new SeqMat<IssmPDouble>(index,out.numberoftriangles,3,1.0); 195 195 *pindex=index_matrix; 196 196 197 197 segments_matrix=new SeqMat<IssmPDouble>(segments,out.numberofsegments,3,1.0); 198 198 *psegments=segments_matrix; -
issm/trunk-jpl/src/c/modules/TriaSearchx/TriaSearchx.cpp
r12832 r13622 40 40 //Get current point coordinates 41 41 r.x=x0[i]; r.y=y0[i]; 42 42 43 43 I=Th.R2ToI2(r); 44 44 -
issm/trunk-jpl/src/c/modules/UpdateConstraintsx/UpdateConstraintsx.cpp
r12520 r13622 23 23 /*start module: */ 24 24 if(VerboseModule()) _pprintLine_(" Updating constraints for time: " << time); 25 25 26 26 /*First, update dof constraints in nodes, using constraints: */ 27 27 SpcNodesx(nodes,constraints,parameters,analysis_type); -
issm/trunk-jpl/src/c/modules/UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.cpp
r13216 r13622 11 11 12 12 void UpdateDynamicConstraintsx(Constraints* constraints,Nodes* nodes,Parameters* parameters,Vector<IssmDouble>* yg){ 13 13 14 14 int configuration_type; 15 15 IssmDouble* yg_serial=NULL; -
issm/trunk-jpl/src/c/modules/VecMergex/VecMergex.cpp
r13216 r13622 18 18 /*retrieve parameters: */ 19 19 parameters->FindParam(&configuration_type,ConfigurationTypeEnum); 20 20 21 21 /*serialize uf: */ 22 22 uf_serial=uf->ToMPISerial(); 23 24 23 25 24 /*Do we have any nodes for this configuration? :*/ -
issm/trunk-jpl/src/c/python/io/CheckNumPythonArguments.cpp
r13056 r13622 8 8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 9 9 #endif 10 11 10 12 11 #define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol … … 23 22 /*figure out size of tuple in input: */ 24 23 size=PyTuple_Size(inputs); 25 24 26 25 /*check on requested size: */ 27 26 if (size==0){ -
issm/trunk-jpl/src/c/python/io/FetchPythonData.cpp
r13484 r13622 47 47 48 48 bool boolean; 49 49 50 50 /*check this is indeed a subtype of long type: */ 51 51 if(!PyBool_Check(py_boolean))_error_("expecting a boolean in input!"); … … 56 56 /*simple copy: */ 57 57 *pboolean=boolean; 58 58 59 59 } 60 60 /*}}}*/ … … 74 74 dims=PyArray_DIMS((PyArrayObject*)py_matrix); 75 75 M=dims[0]; N=dims[1]; 76 76 77 77 if (M && N) { 78 78 /*retrieve internal value: */ … … 110 110 dims=PyArray_DIMS((PyArrayObject*)py_matrix); 111 111 M=dims[0]; N=dims[1]; 112 112 113 113 if (M && N) { 114 114 /*retrieve internal value: */ … … 143 143 dims=PyArray_DIMS((PyArrayObject*)py_vector); 144 144 M=dims[0]; 145 145 146 146 if (M) { 147 147 /*retrieve internal value: */ … … 336 336 /*convert from bytes to string: */ 337 337 string=PyBytes_AS_STRING(py_bytes); 338 338 339 339 *pstring=string; 340 340 } -
issm/trunk-jpl/src/c/python/io/WritePythonData.cpp
r13447 r13622 24 24 /*FUNCTION WriteData(PyObject* py_tuple,int index,int integer){{{*/ 25 25 void WriteData(PyObject* py_tuple, int index, int integer){ 26 26 27 27 PyTuple_SetItem(py_tuple, index, PyInt_FromSsize_t(integer)); 28 28 … … 30 30 /*FUNCTION WriteData(PyObject* py_tuple,int index,char* string){{{*/ 31 31 void WriteData(PyObject* py_tuple, int index, char* string){ 32 32 33 33 PyTuple_SetItem(py_tuple, index, PyUnicode_FromString(string)); 34 34 … … 48 48 /*FUNCTION WriteData(PyObject* py_tuple,int index){{{*/ 49 49 void WriteData(PyObject* py_tuple, int index){ 50 50 51 51 PyTuple_SetItem(py_tuple, index, Py_None); 52 52 … … 118 118 /*FUNCTION WriteData(PyObject* py_tuple,int index,SeqMat<double>* matrix){{{*/ 119 119 void WriteData(PyObject* py_tuple,int index,SeqMat<double>* matrix){ 120 120 121 121 int M,N; 122 122 double* buffer=NULL; 123 123 npy_intp dims[2]={0,0}; 124 124 PyObject* array=NULL; 125 125 126 126 buffer=matrix->ToSerial(); 127 127 matrix->GetSize(&M,&N); … … 129 129 dims[1]=(npy_intp)N; 130 130 array=PyArray_SimpleNewFromData(2,dims,NPY_DOUBLE,buffer); 131 131 132 132 PyTuple_SetItem(py_tuple, index, array); 133 133 … … 135 135 /*FUNCTION WriteData(PyObject* py_tuple,int index,SeqVec<double>* vector){{{*/ 136 136 void WriteData(PyObject* py_tuple,int index,SeqVec<double>* vector){ 137 137 138 138 int M; 139 139 double* buffer=NULL; 140 140 npy_intp dim=10; 141 141 PyObject* array=NULL; 142 142 143 143 buffer=vector->ToMPISerial(); 144 144 vector->GetSize(&M); 145 145 dim=(npy_intp)M; 146 146 array=PyArray_SimpleNewFromData(1,&dim,NPY_DOUBLE,buffer); 147 147 148 148 PyTuple_SetItem(py_tuple, index, array); 149 149 } -
issm/trunk-jpl/src/c/shared/Alloc/alloc.cpp
r13251 r13622 11 11 * the allocation routines described here do not work. 12 12 */ 13 14 13 15 14 #ifdef HAVE_CONFIG_H … … 28 27 29 28 void xdelete(Matrix<IssmDouble>** pv){ 30 29 31 30 if (pv && *pv){ 32 31 /*There is no mxDelete in the Matlab API -> using delete trips up Matlab. So we -
issm/trunk-jpl/src/c/shared/Elements/GetGlobalDofList.cpp
r12437 r13622 9 9 int* ndof_list=NULL; 10 10 int *doflist = NULL; 11 12 11 13 12 if(numnodes){ -
issm/trunk-jpl/src/c/shared/Elements/GetLocalDofList.cpp
r12437 r13622 15 15 ndof_list=xNew<int>(numnodes); 16 16 ngdof_list_cumulative=xNew<int>(numnodes); 17 18 17 19 18 /*Get number of dofs per node, and total for this given set*/ -
issm/trunk-jpl/src/c/shared/Elements/Paterson.cpp
r12475 r13622 10 10 11 11 IssmDouble Paterson(IssmDouble temperature){ 12 12 13 13 /*output: */ 14 14 IssmDouble B; … … 29 29 // % fittedmodel=fit(Temp,B,'cubicspline'); 30 30 // % B=fittedmodel(temperature); 31 32 31 33 32 if(T<=-45.0){ … … 70 69 return B; 71 70 } 72 73 74 -
issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp
r12988 r13622 11 11 12 12 int iqj,imonth, j; 13 13 14 14 IssmDouble saccu; // yearly surface accumulation 15 15 IssmDouble smelt; // yearly melt … … 19 19 IssmDouble runoff; //meltwater only, does not include rain 20 20 IssmDouble sconv; //rhow_rain/rhoi / 12 months 21 21 22 22 IssmDouble density; 23 23 IssmDouble lapser=6.5/1000., sealev=0.; // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper … … 27 27 IssmDouble st; // elevation between altitude of the temp record and current altitude 28 28 IssmDouble sp; // elevation between altitude of the prec record and current altitude 29 29 30 30 // PDD and PD constants and variables 31 31 IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm … … 35 35 IssmDouble DT = 0.02; 36 36 IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow 37 37 38 38 IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist. 39 39 IssmDouble qm = 0.; // snow part of the precipitation … … 42 42 IssmDouble pdd = 0.; 43 43 IssmDouble frzndd = 0.; 44 44 45 45 IssmDouble tstar; // monthly mean surface temp 46 46 IssmDouble Tsum= 0.; // average summer (JJA) temperature 47 47 IssmDouble Tsurf = 0.; // average annual temperature 48 49 48 50 49 IssmDouble deltm=1./12.; 51 50 int ismon[12]={11,0,1,2,3,4,5,6,7,8,9,10}; 52 51 53 52 IssmDouble snwm; // snow that could have been melted in a year. 54 53 IssmDouble snwmf; // ablation factor for snow per positive degree day. 55 54 IssmDouble smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002). 56 55 57 56 IssmDouble dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr 58 57 IssmDouble supice,supcap,diffndd; 59 58 IssmDouble fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice 60 59 IssmDouble pddtj, hmx2; 61 60 62 61 sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months 63 62 … … 68 67 siglim0c = siglimc/DT + 0.5; 69 68 PDup = siglimc+PDCUT; 70 69 71 70 // seasonal loop 72 71 for (iqj = 0; iqj < 12; iqj++){ 73 72 imonth = ismon[iqj]; 74 73 75 74 st=(s-s0t)/1000.; 76 75 tstar = monthlytemperatures[imonth] - lapser *max(st,sealev); 77 76 Tsurf = tstar*deltm+Tsurf; 78 77 79 78 /*********compute PD ****************/ 80 79 if (tstar < PDup){ … … 83 82 else { 84 83 pd = 0.;} 85 84 86 85 /******exp des/elev precip reduction*******/ 87 86 sp=(s-s0p)/1000.; // deselev effect is wrt chng in topo 88 87 if (sp>0.0){q = exp(-desfac*sp);} 89 88 else {q = 1.0;} 90 89 91 90 qmt= qmt + monthlyprec[imonth]*sconv; //*sconv to convert in m of ice equivalent per month 92 91 qmpt= q*monthlyprec[imonth]*sconv; 93 92 qmp= qmp + qmpt; 94 93 qm= qm + qmpt*pd; 95 94 96 95 /*********compute PDD************/ 97 96 // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of … … 105 104 else{frzndd = frzndd - tstar*deltm; } 106 105 } // end of seasonal loop 107 106 108 107 //****************************************************************** 109 108 saccu = qm; 110 109 prect = qmp; // total precipitation during 1 year taking into account des. ef. 111 110 Tsum=Tsum/3; 112 111 113 112 /***** determine PDD factors *****/ 114 113 if(Tsum< -1.) { … … 126 125 snwmf=0.95*snwmf; 127 126 smf=0.95*smf; 128 127 129 128 /***** compute PDD ablation and refreezing *****/ 130 129 pddt = pdd *365; 131 130 snwm = snwmf*pddt; // snow that could have been melted in a year 132 131 hmx2 = min(h,dfrz); // refreeze active layer max depth: dfrz 133 132 134 133 if(snwm < saccu) { 135 134 water=prect-saccu + snwm; //water=rain + snowmelt … … 147 146 // hold the meltwater around for refreezing? And melt-time will have 148 147 // low seasonal frzndd 149 148 150 149 // Superimposed ice : Pfeffer et al. 1991, Tarasov 2002 151 150 152 151 supice= min(hmx2*CovrLm*frzndd+2.2*(saccu-snwm), water); // superimposed ice 153 152 supcap=min(2.2*(saccu-snwm),water); … … 172 171 // does supice make sense when H< 0.1m? then d=thermoactive ice layer //// 173 172 // < 0.1 174 173 175 174 // make more sense to just use residual pdd-ndd except that pdd 176 175 // residual not clear yet … … 197 196 if(Tsurf<0) { 198 197 Tsurf= min(Tsurf+fsupT*diffndd , 0.);} 199 198 200 199 B = -smelt+saccu; 201 200 B = B/yts; 202 201 pddtj=pddt; 203 202 204 203 return B; 205 204 } -
issm/trunk-jpl/src/c/shared/Exceptions/Exceptions.cpp
r13612 r13622 29 29 30 30 ErrorException::~ErrorException() throw(){ 31 31 32 32 int num_procs; 33 33 34 34 /*recover num_procs:*/ 35 35 num_procs=IssmComm::GetSize(); … … 45 45 46 46 void ErrorException::Report() const{ 47 47 48 48 int my_rank; 49 49 int num_procs; 50 50 51 51 /*recover my_rank and num_procs:*/ 52 52 my_rank=IssmComm::GetRank(); 53 53 num_procs=IssmComm::GetSize(); 54 55 56 54 57 55 if (function_name=="" || file_line==0){ //WINDOWS -
issm/trunk-jpl/src/c/shared/Exceptions/exprintf.cpp
r12500 r13622 5 5 * ErrorException(exprintf("%s%i\n","test failed for id:",id)); 6 6 */ 7 8 7 9 8 #include <stdarg.h> -
issm/trunk-jpl/src/c/shared/Exp/DomainOutlineWrite.cpp
r13353 r13622 12 12 13 13 int DomainOutlineWrite(int nprof,int* profnvertices,double** pprofx,double** pprofy,bool* closed,char* domainname){ 14 14 15 15 /*Error management: */ 16 16 int noerr=1; … … 33 33 fprintf(fid,"%s %s\n","##","Icon:0"); 34 34 fprintf(fid,"%s %s %s %s\n","#","Points","Count","Value"); 35 35 36 36 /*Write number of profile vertices: */ 37 37 fprintf(fid,"%u %s\n",profnvertices[counter] ,"1."); 38 38 39 39 /*Write next line: */ 40 40 fprintf(fid,"%s %s %s %s %s\n","#","X","pos","Y","pos"); -
issm/trunk-jpl/src/c/shared/Exp/IsInPoly.cpp
r13220 r13622 18 18 int i, j, c = 0; 19 19 double n1, n2, normp, scalar; 20 20 21 21 /*first test, are they colinear? if yes, is the point in the middle of the segment*/ 22 22 if (edgevalue != 2 ){ -
issm/trunk-jpl/src/c/shared/Exp/IsInPolySerial.cpp
r13367 r13622 5 5 #include <math.h> 6 6 7 8 7 #include "./exp.h" 9 10 8 11 9 int IsInPolySerial(double* in,double* xc,double* yc,int numvertices,double* x,double* y,int nods, int edgevalue){ -
issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp
r13056 r13622 40 40 maxiter=optpars->maxiter; 41 41 cm_jump=optpars->cm_jump; 42 42 43 43 /*initialize counter and get response at the boundaries*/ 44 44 iter=0; -
issm/trunk-jpl/src/c/shared/Numerics/IsInputConverged.cpp
r13056 r13622 20 20 /*output: */ 21 21 IssmDouble eps; 22 22 23 23 /*intermediary: */ 24 24 IssmDouble *newvalues = NULL; -
issm/trunk-jpl/src/c/shared/Numerics/OptimalSearch.cpp
r13056 r13622 38 38 x2 =optpars->xmax; 39 39 maxiter =optpars->maxiter; 40 40 41 41 //get the value of the function at the first boundary 42 42 fx1= (*f)(x1,optargs); -
issm/trunk-jpl/src/c/shared/Numerics/OptionsFromAnalysis.cpp
r13056 r13622 16 16 17 17 char* OptionsFromAnalysis(Parameters* parameters,int analysis_type){ 18 18 19 19 /*output: */ 20 20 char* outstring=NULL; -
issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp
r13521 r13622 40 40 } 41 41 42 43 42 IssmDouble UnitConversionScaleFactor(int type_enum){ 44 43 45 44 IssmDouble yts=365.0*24.0*3600.0; 46 45 47 46 IssmDouble scale; 48 47 switch(type_enum){ -
issm/trunk-jpl/src/c/shared/Numerics/cross.cpp
r12476 r13622 20 20 21 21 } 22 -
issm/trunk-jpl/src/c/shared/Sorting/binary_search.cpp
r13579 r13622 47 47 } 48 48 } 49 49 50 50 //did we find the target? 51 51 if (*mid == target) { … … 60 60 /*Assign output pointers:*/ 61 61 *poffset=offset; 62 62 63 63 /*Return result: */ 64 64 return found; -
issm/trunk-jpl/src/c/shared/String/DescriptorIndex.cpp
r13056 r13622 18 18 19 19 int DescriptorIndex(char* root, int* pindex,char* descriptor){ //We assume root has already been allocated, and we just have to copy into it. 20 20 21 21 char * pch=NULL; 22 22 -
issm/trunk-jpl/src/c/shared/Threads/LaunchThread.cpp
r13056 r13622 31 31 pthread_t *threads = NULL; 32 32 pthread_handle *handles = NULL; 33 33 34 34 /*dynamically allocate: */ 35 35 threads=xNew<pthread_t>(num_threads); … … 52 52 } 53 53 } 54 54 55 55 /*Free ressources:*/ 56 56 xDelete<pthread_t>(threads); … … 62 62 handle.id=0; 63 63 handle.num=1; 64 64 65 65 function((void*)&handle); 66 66 #endif -
issm/trunk-jpl/src/c/shared/Threads/PartitionRange.cpp
r9320 r13622 16 16 int i0=-1; 17 17 int i1=-1; 18 18 19 19 int step; 20 20 int i; 21 22 21 23 22 /*distribute elements across threads :*/ … … 29 28 } 30 29 31 32 30 /*Assign output pointers:*/ 33 31 *pi0=i0; -
issm/trunk-jpl/src/c/shared/TriMesh/AssociateSegmentToElement.cpp
r8301 r13622 6 6 7 7 int AssociateSegmentToElement(double** psegments,int nseg, double* index,int nel){ 8 8 9 9 /*Error management: */ 10 10 int i; … … 30 30 return noerr; 31 31 } 32 33 -
issm/trunk-jpl/src/c/shared/TriMesh/OrderSegments.cpp
r11202 r13622 6 6 7 7 int OrderSegments(double** psegments,int nseg, double* index,int nel){ 8 8 9 9 /*Error management: */ 10 10 int i; … … 18 18 /*element indices: */ 19 19 int el; 20 21 20 22 21 /*Recover segments: */ -
issm/trunk-jpl/src/c/shared/TriMesh/SplitMeshForRifts.cpp
r13248 r13622 12 12 x and y of size nodsx1 13 13 segments of size nsegsx3*/ 14 14 15 15 /*Error management: */ 16 16 int noerr=1; 17 17 18 18 int i,j,k,l; 19 19 int node; … … 23 23 int* riftsegments=NULL; 24 24 int* flags=NULL; 25 25 26 26 int NumGridElementListOnOneSideOfRift; 27 27 int* GridElementListOnOneSideOfRift=NULL; … … 55 55 for (i=0;i<nriftsegs;i++){ 56 56 for (j=0;j<2;j++){ 57 57 58 58 node=riftsegments[4*i+j+2]; 59 59 if(flags[node-1]){ … … 66 66 67 67 if(IsGridOnRift(riftsegments,nriftsegs,node)){ 68 68 69 69 DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,node,index,nel); 70 70 71 71 /*Summary: we have for node, a list of elements 72 72 * (GridElementListOnOneSideOfRift, of size -
issm/trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp
r13363 r13622 17 17 /*Does this node belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift, 18 18 *if it belongs to 2 elements, it is on the tip of a rift, or it has already been split across the rift (see below).*/ 19 19 20 20 int i; 21 21 int j; … … 119 119 /*FUNCTION RiftSegmentsFromSegments{{{*/ 120 120 void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments){ 121 121 122 122 int i,counter; 123 123 int el,el2; 124 124 125 125 int nriftsegs; 126 126 int* riftsegments=NULL; … … 179 179 180 180 xDelete<int>(riftsegments_uncompressed); 181 181 182 182 /*Assign output pointers: */ 183 183 *priftsegments=riftsegments; … … 344 344 *psegmentmarkerlist=segmentmarkerlist; 345 345 *pnsegs=nsegs; 346 346 347 347 return noerr; 348 348 }/*}}}*/ … … 384 384 double *segmentmarkerlist = NULL; 385 385 int numsegs; 386 386 387 387 /*output: */ 388 388 int new_numsegs; … … 460 460 /*FUNCTION PairRiftElements{{{*/ 461 461 int PairRiftElements(int** priftsnumpairs, double*** priftspairs,int numrifts,int* riftsnumsegments, double** riftssegments,double* x,double* y){ 462 463 462 464 463 int noerr=1; … … 693 692 int i; 694 693 int noerr=1; 695 694 696 695 /*output: */ 697 696 int riftflag=0; … … 717 716 /*FUNCTION OrderRifts{{{*/ 718 717 int OrderRifts(double** priftstips, double** riftssegments,double** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels){ 719 718 720 719 int noerr=1; 721 720 int i,j,k,counter; … … 747 746 riftpairs=riftspairs[i]; 748 747 numsegs=riftsnumsegments[i]; 749 748 750 749 /*Allocate copy of riftsegments and riftpairs, 751 750 *as well as ordering vector: */ … … 793 792 } 794 793 } 795 794 796 795 if (node4==node2){ 797 796 /*node2 is a tip*/ … … 818 817 node1=(int)*(riftsegments+3*j+0); 819 818 node2=(int)*(riftsegments+3*j+1); 820 819 821 820 if ((node1==node) || (node2==node)){ 822 821 /*Ok, this segment is connected to node, plug its index into order, unless we already plugged it before: */ … … 851 850 *(riftpairs_copy+2*j+1)=*(riftpairs+2*order[j]+1); 852 851 } 853 852 854 853 for (j=0;j<numsegs;j++){ 855 854 *(riftsegments+3*j+0)=*(riftsegments_copy+3*j+0); … … 873 872 int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,double** riftssegments, 874 873 int* riftsnumsegs,double** riftspairs,double* riftstips,double* x,double* y){ 875 876 874 877 875 int noerr=1; … … 906 904 /*allocate riftpenaltypairs, and riftnumpenaltypairs: */ 907 905 if((numsegs/2-1)!=0)riftpenaltypairs=xNewZeroInit<double>((numsegs/2-1)*RIFTPENALTYPAIRSWIDTH); 908 906 909 907 /*Go through only one flank of the rifts, not counting the tips: */ 910 908 counter=0; … … 997 995 *(riftpenaltypairs+j*7+5)=*(riftpenaltypairs+j*7+5)/magnitude; 998 996 } 999 997 1000 998 riftspenaltypairs[i]=riftpenaltypairs; 1001 999 riftsnumpenaltypairs[i]=(numsegs/2-1); … … 1035 1033 y=*py; 1036 1034 nods=*pnods; 1037 1038 1035 1039 1036 for (i=0;i<num_seg;i++){ -
issm/trunk-jpl/src/c/shared/Wrapper/ModuleBoot.cpp
r13231 r13622 13 13 14 14 int ModuleBoot(void){ 15 15 16 16 /*Some test for MPI_Init crash with mpich2 1.4 on larsen, just ignore*/ 17 17 #ifdef _HAVE_PETSC_ -
issm/trunk-jpl/src/c/solutions/diagnostic_core.cpp
r13621 r13622 95 95 } 96 96 97 98 97 if(save_results){ 99 98 if(VerboseSolution()) _pprintLine_(" saving results"); -
issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp
r13325 r13622 23 23 Vector<IssmDouble>* df = NULL; 24 24 Vector<IssmDouble>* ys = NULL; 25 25 26 26 Loads* loads=NULL; 27 27 bool converged; -
issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp
r13216 r13622 42 42 femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum); 43 43 UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters); 44 44 45 45 count=1; 46 46 converged=false; … … 56 56 femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); 57 57 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); 58 58 59 59 //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) 60 60 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz); … … 77 77 femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum); 78 78 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); 79 79 80 80 /*solve: */ 81 81 SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert, &df_vert,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); -
issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp
r13216 r13622 35 35 bool lowmem=0; 36 36 int configuration_type; 37 38 37 39 38 /*Recover parameters: */ … … 76 75 77 76 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); 78 77 79 78 if(converged)break; 80 79 } -
issm/trunk-jpl/src/c/toolkits/mpi/patches/DetermineLocalSize.cpp
r13612 r13622 26 26 MPI_Comm_size(comm,&num_procs); 27 27 28 29 28 /*We are not bound by any library, just use what seems most logical*/ 30 29 num_local_rows=xNew<int>(num_procs); … … 36 35 num_local_rows[i]=(int)floor((double)global_size/(double)num_procs); 37 36 } 38 37 39 38 /*There may be some rows left. Distribute evenly.*/ 40 39 row_rest=global_size - num_procs*(int)floor((double)global_size/(double)num_procs); … … 43 42 } 44 43 local_size=num_local_rows[my_rank]; 45 44 46 45 /*free ressources: */ 47 46 xDelete<int>(num_local_rows); -
issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp
r13596 r13622 46 46 int* idxm=NULL; 47 47 int* idxn=NULL; 48 48 49 49 if(M)idxm=xNew<int>(M); 50 50 if(N)idxn=xNew<int>(N); … … 52 52 for(i=0;i<M;i++)idxm[i]=i; 53 53 for(i=0;i<N;i++)idxn[i]=i; 54 55 54 56 55 this->matrix=NewMat(M,N,sparsity,IssmComm::GetComm()); … … 66 65 /*FUNCTION PetscMat::PetscMat(int M,int N, int connectivity, int numberofdofspernode){{{*/ 67 66 PetscMat::PetscMat(int M,int N, int connectivity,int numberofdofspernode){ 68 67 69 68 this->matrix=NewMat(M,N,connectivity,numberofdofspernode,IssmComm::GetComm()); 70 69 … … 99 98 IssmDouble PetscMat::Norm(NormMode mode){ 100 99 101 102 100 IssmDouble norm=0; 103 101 _assert_(this->matrix); 104 102 MatNorm(this->matrix,ISSMToPetscNormMode(mode),&norm); 105 103 106 104 return norm; 107 105 -
issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp
r13599 r13622 29 29 /*FUNCTION PetscVec::PetscVec(int M,bool fromlocalsize){{{*/ 30 30 PetscVec::PetscVec(int M,bool fromlocalsize){ 31 31 32 32 this->vector=NewVec(M,IssmComm::GetComm(),fromlocalsize); 33 33 … … 79 79 /*FUNCTION PetscVec::Assemble{{{*/ 80 80 void PetscVec::Assemble(void){ 81 81 82 82 _assert_(this->vector); 83 83 VecAssemblyBegin(this->vector); … … 88 88 /*FUNCTION PetscVec::SetValues{{{*/ 89 89 void PetscVec::SetValues(int ssize, int* list, IssmDouble* values, InsMode mode){ 90 90 91 91 _assert_(this->vector); 92 92 VecSetValues(this->vector,ssize,list,values,ISSMToPetscInsertMode(mode)); … … 128 128 /*FUNCTION PetscVec::Duplicate{{{*/ 129 129 PetscVec* PetscVec::Duplicate(void){ 130 130 131 131 PetscVec* output=NULL; 132 132 _assert_(this->vector); … … 165 165 /*FUNCTION PetscVec::ToMPISerial{{{*/ 166 166 IssmDouble* PetscVec::ToMPISerial(void){ 167 167 168 168 IssmDouble* vec_serial=NULL; 169 169 VecToMPISerial(&vec_serial, this->vector,IssmComm::GetComm()); -
issm/trunk-jpl/src/c/toolkits/petsc/patches/GetOwnershipBoundariesFromRange.cpp
r13612 r13622 3 3 * lower row and upper row from a matrix a cpu owns. 4 4 */ 5 6 5 7 6 #ifdef HAVE_CONFIG_H … … 19 18 int my_rank; 20 19 int num_procs; 21 20 22 21 /*recover my_rank and num_procs:*/ 23 22 MPI_Comm_size(comm,&num_procs); … … 26 25 /*output: */ 27 26 int lower_row,upper_row; 28 27 29 28 /*intermediary :*/ 30 29 int i; -
issm/trunk-jpl/src/c/toolkits/petsc/patches/ISSMToPetscInsertMode.cpp
r13056 r13622 2 2 * \brief: convert InsertMode from ISSM to Petsc 3 3 */ 4 5 4 6 5 #ifdef HAVE_CONFIG_H … … 19 18 #include "../../../shared/shared.h" 20 19 21 22 20 InsertMode ISSMToPetscInsertMode(InsMode mode){ 23 21 -
issm/trunk-jpl/src/c/toolkits/petsc/patches/ISSMToPetscMatrixType.cpp
r13056 r13622 2 2 * \brief: convert MatrixType from ISSM to Petsc 3 3 */ 4 5 4 6 5 #ifdef HAVE_CONFIG_H … … 19 18 #include "../../../shared/shared.h" 20 19 21 22 20 MatType ISSMToPetscMatrixType(MatrixType type){ 23 21 -
issm/trunk-jpl/src/c/toolkits/petsc/patches/ISSMToPetscNormMode.cpp
r13056 r13622 2 2 * \brief: convert NormMode from ISSM to Petsc 3 3 */ 4 5 4 6 5 #ifdef HAVE_CONFIG_H … … 19 18 #include "../../../shared/shared.h" 20 19 21 22 20 NormType ISSMToPetscNormMode(NormMode mode){ 23 21 -
issm/trunk-jpl/src/c/toolkits/petsc/patches/KSPFree.cpp
r9826 r13622 8 8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 9 9 #endif 10 11 10 12 11 /*Petsc includes: */ … … 26 25 27 26 } 28 -
issm/trunk-jpl/src/c/toolkits/petsc/patches/MatFree.cpp
r9826 r13622 8 8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 9 9 #endif 10 11 10 12 11 /*Petsc includes: */ … … 26 25 27 26 } 28 -
issm/trunk-jpl/src/c/toolkits/petsc/patches/MatMultPatch.cpp
r13612 r13622 43 43 44 44 int MatMultCompatible(Mat A,Vec x,COMM comm){ 45 45 46 46 /*error management*/ 47 47 48 48 int local_m,local_n; 49 49 int lower_row,upper_row,range; … … 51 51 int sumresult; 52 52 int num_procs; 53 53 54 54 /*recover num_procs:*/ 55 55 MPI_Comm_size(comm,&num_procs); … … 57 57 MatGetLocalSize(A,&local_m,&local_n);; 58 58 VecGetLocalSize(x,&range);; 59 59 60 60 if (local_n!=range)result=0; 61 61 62 62 /*synchronize result: */ 63 63 MPI_Reduce (&result,&sumresult,1,MPI_INT,MPI_SUM,0,comm ); … … 81 81 /*output: */ 82 82 Vec outvector=NULL; 83 83 84 84 /*Create outvector with local size m*/ 85 85 VecCreate(comm,&outvector); ; -
issm/trunk-jpl/src/c/toolkits/petsc/patches/MatToSerial.cpp
r13612 r13622 9 9 #endif 10 10 11 12 11 #include "../petscincludes.h" 13 12 #include "../../../shared/shared.h" … … 18 17 int my_rank; 19 18 int num_procs; 20 19 21 20 /*Petsc variables*/ 22 21 PetscInt lower_row,upper_row; … … 28 27 double* local_matrix=NULL; /*matrix local to each node used for temporary holding matrix values*/ 29 28 int buffer[3]; 30 29 31 30 /*recover my_rank and num_procs:*/ 32 31 MPI_Comm_size(comm,&num_procs); … … 35 34 /*Output*/ 36 35 double* outmatrix=NULL; 37 36 38 37 /*get matrix size: */ 39 38 MatGetSize(matrix,&M,&N); … … 43 42 upper_row--; 44 43 range=upper_row-lower_row+1; 45 44 46 45 /*Local and global allocation*/ 47 46 if (my_rank==0)outmatrix=xNew<double>(M*N); 48 47 49 48 if (range){ 50 49 local_matrix=xNew<double>(N*range); 51 50 idxm=xNew<int>(range); 52 51 idxn=xNew<int>(N); 53 52 54 53 for (i=0;i<N;i++){ 55 54 *(idxn+i)=i; … … 58 57 *(idxm+i)=lower_row+i; 59 58 } 60 59 61 60 MatGetValues(matrix,range,idxm,N,idxn,local_matrix); 62 61 } … … 64 63 /*Now each node holds its local_matrix containing range rows. 65 64 * We send these rows to the matrix on node 0*/ 66 65 67 66 for (i=1;i<num_procs;i++){ 68 67 if (my_rank==i){ … … 82 81 memcpy(outmatrix,local_matrix,N*range*sizeof(double)); 83 82 } 84 83 85 84 /*Assign output pointer: */ 86 85 *poutmatrix=outmatrix; -
issm/trunk-jpl/src/c/toolkits/petsc/patches/NewMat.cpp
r13602 r13622 34 34 m=DetermineLocalSize(M,comm); 35 35 n=DetermineLocalSize(N,comm); 36 36 37 37 nnz=(int)((double)M*(double)N*sparsity); //number of non zeros. 38 38 d_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2 … … 62 62 m=DetermineLocalSize(M,comm); 63 63 n=DetermineLocalSize(N,comm); 64 64 65 65 nnz=(int)((double)M*(double)N*sparsity); //number of non zeros. 66 66 d_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2 … … 112 112 /*preallocation according to type: */ 113 113 MatGetType(outmatrix,&type); 114 114 115 115 #if _PETSC_MAJOR_ == 2 116 116 if((strcmp(type,"mpiaij")==0) || (strcmp(type,"aijmumps")==0)){ -
issm/trunk-jpl/src/c/toolkits/petsc/patches/NewVec.cpp
r13602 r13622 32 32 local_size=DetermineLocalSize(size,comm); 33 33 } 34 34 35 35 VecCreate(comm,&vector); 36 36 37 37 VecSetSizes(vector,local_size,PETSC_DECIDE); 38 38 VecSetFromOptions(vector); -
issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp
r12431 r13622 2 2 * \brief: convert a sparse or dense Petsc matrix into a matlab matrix 3 3 */ 4 5 4 6 5 #ifdef HAVE_CONFIG_H … … 19 18 /*Petsc includes: */ 20 19 #include "../../../shared/shared.h" 21 22 20 23 21 void PetscMatrixToDoubleMatrix(double** pmatrix, int* prows, int* pcols,Mat petsc_matrix){ -
issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp
r11466 r13622 53 53 } 54 54 55 56 55 #if _PETSC_MAJOR_ >= 3 57 56 PetscOptionsGetString(PETSC_NULL,"-pc_factor_mat_solver_package",&option[0],100,&flag); … … 61 60 #endif 62 61 63 64 62 PetscOptionsGetString(PETSC_NULL,"-issm_option_solver",&option[0],100,&flag); 65 63 if (strcmp(option,"stokes")==0){ -
issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp
r13571 r13622 21 21 void PetscOptionsInsertMultipleString(char* options_string){ 22 22 23 24 23 /*The list of options is going to be pairs of the type "-option option_value"*/ 25 24 #if _PETSC_MAJOR_ == 2 … … 35 34 int first_token=1; 36 35 37 38 36 PetscTokenCreate(options_string,' ',&token); 39 37 for (;;){ 40 41 38 42 39 /*Read next tokens*/ … … 45 42 } 46 43 PetscTokenFind(token,&second); 47 44 48 45 if (!first){ 49 46 /*We are at the end of options*/ -
issm/trunk-jpl/src/c/toolkits/petsc/patches/VecDuplicatePatch.cpp
r9320 r13622 8 8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 9 9 #endif 10 11 10 12 11 /*Petsc includes: */ -
issm/trunk-jpl/src/c/toolkits/petsc/patches/VecFree.cpp
r9826 r13622 8 8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 9 9 #endif 10 11 10 12 11 /*Petsc includes: */ … … 25 24 26 25 } 27 -
issm/trunk-jpl/src/c/toolkits/petsc/patches/VecMerge.cpp
r13056 r13622 22 22 23 23 int i; 24 24 25 25 /*Petsc matrix*/ 26 26 int lower_row,upper_row,range; … … 57 57 VecSetValues(A,range,idxm,values,INSERT_VALUES); 58 58 } 59 59 60 60 /*Assemble vector*/ 61 61 VecAssemblyBegin(A); -
issm/trunk-jpl/src/c/toolkits/petsc/patches/VecToMPISerial.cpp
r13612 r13622 9 9 #endif 10 10 11 12 11 #include "../petscincludes.h" 13 12 #include "../../../shared/shared.h" 14 13 15 14 int VecToMPISerial(double** pgathered_vector, Vec vector,COMM comm){ 16 15 17 16 int i; 18 17 int num_procs; … … 25 24 int * idxn=NULL; 26 25 int buffer[3]; 27 26 28 27 /*intermediary results*/ 29 28 double* local_vector=NULL; … … 31 30 /*input*/ 32 31 int vector_size; 33 32 34 33 /*Output*/ 35 34 double* gathered_vector=NULL; //Global vector holding the final assembled vector on all nodes. 36 35 37 36 /*recover my_rank and num_procs*/ 38 37 MPI_Comm_size(comm,&num_procs); … … 47 46 /*Allocate gathered vector on all nodes .*/ 48 47 gathered_vector=xNew<double>(vector_size); 49 48 50 49 /*Allocate local vectors*/ 51 50 VecGetOwnershipRange(vector,&lower_row,&upper_row); … … 89 88 /*Assign output pointers: */ 90 89 *pgathered_vector=gathered_vector; 91 90 92 91 /*free ressources: */ 93 92 xDelete<int>(idxn); 94 93 xDelete<double>(local_vector); 95 94 96 95 return 1; 97 96 } -
issm/trunk-jpl/src/c/toolkits/plapack/patches/CyclicalFactorization.cpp
r12431 r13622 3 3 */ 4 4 #include <math.h> 5 6 5 7 6 #include "../../../shared/shared.h" … … 10 9 11 10 int CyclicalFactorization(int* pnprows,int* pnpcols,int num_procs){ 12 11 13 12 int nprows,npcols; 14 13 int last_diff; … … 65 64 */ 66 65 int SmallestPrimeFactor(int* output,int input){ 67 66 68 67 int found=0; 69 68 int i; -
issm/trunk-jpl/src/c/toolkits/plapack/patches/PlapackInvertMatrix.cpp
r13598 r13622 14 14 15 15 void PlapackInvertMatrixLocalCleanup(PLA_Obj* pa,PLA_Template* ptempl,double** parrayA,int** pidxnA,MPI_Comm* pcomm_2d); 16 16 17 17 int PlapackInvertMatrix(Mat* A,Mat* inv_A,int status,int con,COMM comm){ 18 18 /*inv_A does not yet exist, inv_A was just allocated, that's all*/ … … 26 26 int lower_row,upper_row,range,this_range,this_lower_row; 27 27 MatType type; 28 28 29 29 /*Plapack: */ 30 30 MPI_Datatype datatype; … … 42 42 int* idxnA=NULL; 43 43 int d_nz,o_nz; 44 44 45 45 /*Feedback to client*/ 46 46 int computation_status; … … 92 92 /* Set the datatype */ 93 93 datatype = MPI_DOUBLE; 94 94 95 95 /* Copy A into a*/ 96 96 PLA_Matrix_create(datatype,mA,nA,templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&a); … … 127 127 xDelete<double>(arrayA); 128 128 xDelete<int>(idxnA); 129 129 130 130 /*Finalize PLAPACK*/ 131 131 PLA_Finalize(); -
issm/trunk-jpl/src/c/toolkits/plapack/patches/PlapackToPetsc.cpp
r13598 r13622 8 8 int PlapackToPetsc(Mat* A,int local_mA,int local_nA,int mA,int nA,MatType type,PLA_Obj a,PLA_Template templ,int nprows,int npcols,int nb,COMM comm){ 9 9 10 11 10 int i; 12 11 … … 22 21 double* local_buffer=NULL; 23 22 24 25 23 /*Create matrix A (right now, we just have an allocated pointer A*/ 26 24 if (strcasecmp_eq(type,MATMPIAIJ)){ … … 33 31 MatCreateMPIDense(comm,local_mA,local_nA, mA,nA,PETSC_NULL,A); 34 32 } 35 33 36 34 MatGetOwnershipRange(*A,&lower_row,&upper_row); 37 35 upper_row--; 38 36 range=upper_row-lower_row+1; 39 37 40 38 /*Build the Plapack row and column indices corresponding to the local_buffer stored in a. 41 39 We need those indices to directly plug local_buffer into the Petsc matrix A. We do not … … 43 41 problem size becomes big. We rely therefore on MatAssembly from Petsc to gather the plapack 44 42 matrix into a Petsc Matrix.*/ 45 43 46 44 /*Vector physically based block cyclic distribution: */ 47 45 row_nodes=xNew<int>(mA); … … 82 80 /*Get local buffer: */ 83 81 PLA_Obj_local_buffer(a,(void**)&local_buffer); 84 82 85 83 /*Insert into invA matrix. Use col oriented insertion, for Plapack is column oriented*/ 86 84 MatSetOption(*A,MAT_COLUMN_ORIENTED);
Note:
See TracChangeset
for help on using the changeset viewer.