Changeset 6029
- Timestamp:
- 09/24/10 13:23:32 (14 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified issm/trunk/src/c/objects/Elements/Tria.cpp ¶
r5989 r6029 718 718 Ke=CreateKMatrixPrognostic(); 719 719 break; 720 case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:720 case BalancedthicknessAnalysisEnum: 721 721 Ke=CreateKMatrixBalancedthickness(); 722 break; 723 case AdjointBalancedthicknessAnalysisEnum: 724 Ke=CreateKMatrixAdjointBalancedthickness(); 722 725 break; 723 726 case BalancedvelocitiesAnalysisEnum: … … 2384 2387 2385 2388 /*Tria specific routines: */ 2389 /*FUNCTION Tria::CreateKMatrixAdjointBalancedthickness {{{1*/ 2390 ElementMatrix* Tria::CreateKMatrixAdjointBalancedthickness(void){ 2391 2392 ElementMatrix* Ke=NULL; 2393 2394 /*Get Element Matrix of the forward model*/ 2395 switch(GetElementType()){ 2396 case P1Enum: 2397 Ke=CreateKMatrixBalancedthickness_CG(); 2398 break; 2399 case P1DGEnum: 2400 Ke=CreateKMatrixBalancedthickness_DG(); 2401 break; 2402 default: 2403 ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType())); 2404 } 2405 2406 /*Transpose and return Ke*/ 2407 Ke->Transpose(); 2408 return Ke; 2409 } 2410 /*}}}*/ 2386 2411 /*FUNCTION Tria::CreateKMatrixBalancedthickness {{{1*/ 2387 2412 ElementMatrix* Tria::CreateKMatrixBalancedthickness(void){ -
TabularUnified issm/trunk/src/c/objects/Elements/Tria.h ¶
r5933 r6029 121 121 /*}}}*/ 122 122 /*Tria specific routines:{{{1*/ 123 ElementMatrix* CreateKMatrixAdjointBalancedthickness(void); 123 124 ElementMatrix* CreateKMatrixBalancedthickness(void); 124 125 ElementMatrix* CreateKMatrixBalancedthickness_DG(void); -
TabularUnified issm/trunk/src/c/objects/Loads/Numericalflux.cpp ¶
r6026 r6029 509 509 510 510 /*Initialize Element matrix and return if necessary*/ 511 ElementMatrix* Ke = NULL; 511 512 Tria* tria=(Tria*)element; 512 513 if(tria->IsOnWater()) return NULL; 513 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);514 514 515 515 /*Retrieve all inputs and parameters*/ … … 532 532 UdotN=mean_vx*normal[0]+mean_vy*normal[1]; 533 533 if (UdotN<=0){ 534 /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ 535 return NULL; 534 return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ 535 } 536 else{ 537 Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters); 536 538 } 537 539 … … 659 661 660 662 /*Initialize Element matrix and return if necessary*/ 663 ElementMatrix* Ke = NULL; 661 664 Tria* tria=(Tria*)element; 662 665 if(tria->IsOnWater()) return NULL; 663 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);664 666 665 667 /*Retrieve all inputs and parameters*/ … … 681 683 UdotN=mean_vx*normal[0]+mean_vy*normal[1]; 682 684 if (UdotN<=0){ 683 /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ 684 return NULL; 685 return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ 686 } 687 else{ 688 Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters); 685 689 } 686 690 … … 731 735 ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthicknessInternal(void){ 732 736 733 return CreateKMatrixBalancedthicknessInternal(); 737 ElementMatrix* Ke=CreateKMatrixBalancedthicknessInternal(); 738 if (Ke) Ke->Transpose(); 739 return Ke; 734 740 } 735 741 /*}}}*/ … … 737 743 ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthicknessBoundary(void){ 738 744 739 return CreateKMatrixBalancedthicknessBoundary(); 745 ElementMatrix* Ke=CreateKMatrixBalancedthicknessBoundary(); 746 if(Ke) Ke->Transpose(); 747 return Ke; 740 748 } 741 749 /*}}}*/ … … 779 787 780 788 /*Initialize Load Vector and return if necessary*/ 789 ElementVector* pe = NULL; 781 790 Tria* tria=(Tria*)element; 782 791 if(tria->IsOnWater()) return NULL; 783 ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);784 792 785 793 /*Retrieve all inputs and parameters*/ … … 800 808 vyaverage_input->GetParameterValue(&mean_vy,gauss); 801 809 delete gauss; 810 802 811 UdotN=mean_vx*normal[0]+mean_vy*normal[1]; 803 812 if (UdotN>0){ 804 /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/ 805 return NULL; 813 return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/ 814 } 815 else{ 816 ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters); 806 817 } 807 818 … … 868 879 869 880 /*Initialize Load Vector and return if necessary*/ 881 ElementVector* pe = NULL; 870 882 Tria* tria=(Tria*)element; 871 883 if(tria->IsOnWater()) return NULL; 872 ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);873 884 874 885 /*Retrieve all inputs and parameters*/ … … 876 887 Input* vxaverage_input=tria->inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); 877 888 Input* vyaverage_input=tria->inputs->GetInput(VyEnum); ISSMASSERT(vyaverage_input); 878 Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum); 879 880 /*Here, as it is a forcing, we have H=Hobs by default (for control methods)*/ 881 if (!thickness_input){ 882 thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 883 } 889 Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); 884 890 GetNormal(&normal[0],xyz_list); 885 891 … … 895 901 UdotN=mean_vx*normal[0]+mean_vy*normal[1]; 896 902 if (UdotN>0){ 897 /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/ 898 return NULL; 903 return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/ 904 } 905 else{ 906 pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters); 899 907 } 900 908 … … 925 933 ElementVector* Numericalflux::CreatePVectorAdjointBalancedthickness(void){ 926 934 927 int type; 928 inputs->GetParameterValue(&type,TypeEnum); 929 930 switch(type){ 931 case InternalEnum: 932 return CreatePVectorAdjointBalancedthicknessInternal(); 933 case BoundaryEnum: 934 return CreatePVectorAdjointBalancedthicknessBoundary(); 935 default: 936 ISSMERROR("type not supported yet"); 937 } 938 } 939 /*}}}*/ 940 /*FUNCTION Numericalflux::CreatePVectorAdjointBalancedthicknessInternal{{{1*/ 941 ElementVector* Numericalflux::CreatePVectorAdjointBalancedthicknessInternal(void){ 942 943 /*Nothing added to PVector*/ 935 /*No PVector for the Adjoint*/ 944 936 return NULL; 945 946 }947 /*}}}*/948 /*FUNCTION Numericalflux::CreatePVectorAdjointBalancedthicknessBoundary{{{1*/949 ElementVector* Numericalflux::CreatePVectorAdjointBalancedthicknessBoundary(void){950 951 /* constants*/952 const int numdof=NDOF1*NUMVERTICES_BOUNDARY;953 954 /* Intermediaries*/955 int i,j,ig,index1,index2;956 double DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness;957 double xyz_list[NUMVERTICES_BOUNDARY][3];958 double normal[2];959 double L[numdof];960 GaussTria *gauss;961 962 /*Initialize Load Vector and return if necessary*/963 Tria* tria=(Tria*)element;964 if(tria->IsOnWater()) return NULL;965 ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);966 967 /*Retrieve all inputs and parameters*/968 GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY);969 Input* vxaverage_input=tria->inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input);970 Input* vyaverage_input=tria->inputs->GetInput(VyEnum); ISSMASSERT(vyaverage_input);971 Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum); ISSMASSERT(thickness_input);972 GetNormal(&normal[0],xyz_list);973 974 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/975 index1=tria->GetNodeIndex(nodes[0]);976 index2=tria->GetNodeIndex(nodes[1]);977 978 gauss=new GaussTria();979 gauss->GaussEdgeCenter(index1,index2);980 vxaverage_input->GetParameterValue(&mean_vx,gauss);981 vyaverage_input->GetParameterValue(&mean_vy,gauss);982 delete gauss;983 UdotN=mean_vx*normal[0]+mean_vy*normal[1];984 if (UdotN>0){985 /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/986 return NULL;987 }988 989 /* Start looping on the number of gaussian points: */990 gauss=new GaussTria(index1,index2,2);991 for(ig=gauss->begin();ig<gauss->end();ig++){992 993 gauss->GaussPoint(ig);994 995 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);996 997 vxaverage_input->GetParameterValue(&vx,gauss);998 vyaverage_input->GetParameterValue(&vy,gauss);999 thickness_input->GetParameterValue(&thickness,gauss);1000 UdotN=vx*normal[0]+vy*normal[1];1001 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);1002 DL= - gauss->weight*Jdet*UdotN*thickness;1003 1004 for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];1005 }1006 1007 /*Clean up and return*/1008 delete gauss;1009 return pe;1010 937 } 1011 938 /*}}}*/ -
TabularUnified issm/trunk/src/c/objects/Loads/Numericalflux.h ¶
r6026 r6029 90 90 ElementVector* CreatePVectorBalancedthicknessBoundary(void); 91 91 ElementVector* CreatePVectorAdjointBalancedthickness(void); 92 ElementVector* CreatePVectorAdjointBalancedthicknessInternal(void);93 ElementVector* CreatePVectorAdjointBalancedthicknessBoundary(void);94 92 /*}}}*/ 95 93 -
TabularUnified issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp ¶
r6027 r6029 43 43 this->col_sglobaldoflist=NULL; 44 44 45 } 46 /*}}}*/ 47 /*FUNCTION ElementMatrix::ElementMatrix(ElementMatrix* Ke){{{1*/ 48 ElementMatrix::ElementMatrix(ElementMatrix* Ke){ 49 50 if(!Ke) ISSMERROR("Input Element Matrix is a NULL pointer"); 51 this->Init(Ke); 52 return; 45 53 } 46 54 /*}}}*/ … … 319 327 320 328 /*Intermediaries*/ 321 double *values_copy=NULL; 322 int temp; 329 ElementMatrix* Ke_copy=new ElementMatrix(this); 330 331 /*Update sizes*/ 332 this->nrows=Ke_copy->ncols; 333 this->ncols=Ke_copy->nrows; 334 335 /*Transpose values*/ 336 for (int i=0;i<this->nrows;i++) for(int j=0;j<this->ncols;j++) this->values[i*this->ncols+j]=Ke_copy->values[j*Ke_copy->ncols+i]; 323 337 324 338 /*Transpose indices*/ … … 327 341 } 328 342 329 /*Transpose values*/330 values_copy=(double*)xmalloc(this->nrows*this->ncols*sizeof(double));331 memcpy(values_copy,this->values,this->nrows*this->ncols*sizeof(double));332 for (int i=0;i<this->nrows;i++) for(int j=0;j<this->ncols;j++) this->values[j*this->nrows+i]=values_copy[i*this->ncols+j];333 334 /*Update sizes*/335 temp=this->nrows;336 this->nrows=this->ncols;337 this->ncols=temp;338 339 343 /*Clean up and return*/ 340 xfree((void**)&values_copy); 344 delete Ke_copy; 345 return; 341 346 } 342 347 /*}}}*/ -
TabularUnified issm/trunk/src/c/objects/Numerics/ElementMatrix.h ¶
r6027 r6029 52 52 /*ElementMatrix constructors, destructors {{{1*/ 53 53 ElementMatrix(); 54 ElementMatrix(ElementMatrix* Ke); 54 55 ElementMatrix(ElementMatrix* Ke1,ElementMatrix* Ke2); 55 56 ElementMatrix(ElementMatrix* Ke1,ElementMatrix* Ke2,ElementMatrix* Ke3);
Note:
See TracChangeset
for help on using the changeset viewer.