![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 22360
- Timestamp:
- 01/18/18 08:46:05 (7 years ago)
- Location:
- issm/trunk-jpl/src/c/classes
- Files:
-
- 2 edited
-
FemModel.cpp (modified) (11 diffs)
-
FemModel.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/FemModel.cpp
r22349 r22360 91 91 #if !defined(_HAVE_ADOLC_) 92 92 if(amr_frequency){ 93 /*Verifications. AMR supports SSA, P1 and horizontal 2D domain*/ 94 bool isSSA; 95 int domaintype,element_type,analysis_counter=-1; 96 this->parameters->FindParam(&isSSA,FlowequationIsSSAEnum); 97 this->parameters->FindParam(&domaintype,DomainTypeEnum); 98 for(int i=0;i<this->nummodels;i++) if(this->analysis_type_list[i]==StressbalanceAnalysisEnum){analysis_counter=i;break;} 99 if(analysis_counter==-1) _error_("Could not find alias for analysis_type StressbalanceAnalysisEnum in list of FemModel analyses\n"); 100 for(int i=0;i<this->elements->Size();i++){ 101 Element* element = xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); 102 element_type = element->element_type_list[analysis_counter]; 103 if(element_type!=P1Enum) _error_("Element type "<<EnumToStringx(element_type)<<" not supported with AMR yet!\n"); 104 } 105 if(!isSSA) _error_("Flow equation not supported with AMR yet!\n "); 106 if(domaintype!=Domain2DhorizontalEnum) _error_("Domain "<<EnumToStringx(domaintype)<<" not supported with AMR yet!\n"); 107 93 108 this->parameters->FindParam(&amrtype,AmrTypeEnum); 94 109 switch(amrtype){ … … 2456 2471 default: _error_("not implemented yet"); 2457 2472 } 2458 2473 2459 2474 /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/ 2460 2475 this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices); … … 2477 2492 /*Creating nodes and constraints*/ 2478 2493 /*Just SSA (2D) and P1 in this version*/ 2479 Nodes* new_nodes =new Nodes();2480 Constraints* new_constraints =new Constraints();2481 2482 int nodecounter =0;2494 Nodes* new_nodes = new Nodes(); 2495 Constraints* new_constraints = new Constraints(); 2496 2497 int nodecounter =0; 2483 2498 int constraintcounter=0; 2484 2499 for(int i=0;i<this->nummodels;i++){//create nodes for each analysis in analysis_type_list … … 2487 2502 2488 2503 /*As the domain is 2D, it is not necessary to create nodes for this analysis*/ 2489 /*itapopo must verify if domain is not 3D. Only 2D in this version!*/ 2490 if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; 2504 if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; 2491 2505 2492 2506 this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes); 2493 if(analysis_enum==StressbalanceAnalysisEnum) this->CreateConstraints(newnumberofvertices,newnumberofelements,nodecounter,constraintcounter,newx,newy,my_vertices,new_constraints);2507 this->CreateConstraints(new_vertices,nodecounter,constraintcounter,analysis_enum,new_constraints); 2494 2508 this->UpdateElements(newnumberofelements,newelementslist,my_elements,nodecounter,i,new_elements); 2495 2509 … … 2660 2674 } 2661 2675 /*}}}*/ 2662 void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/ 2663 2664 int maxinputs = MaximumNumberOfDefinitionsEnum; 2676 void FemModel::GetInputs(int* pnumP0inputs,IssmDouble** pP0inputs,int** pP0input_enums,int** pP0input_interp,int* pnumP1inputs,IssmDouble** pP1inputs,int** pP1input_enums,int** pP1input_interp){/*{{{*/ 2677 2678 int maxinputs = MaximumNumberOfDefinitionsEnum; 2679 int numberofvertices = this->vertices->NumberOfVertices(); 2680 int numberofelements = this->elements->NumberOfElements(); 2681 int elementswidth = this->GetElementsWidth(); 2682 int numP0inputs = -1; 2683 IssmDouble* P0inputs = NULL; 2684 Vector<IssmDouble>* vP0inputs = NULL; 2685 int* P0input_enums = NULL; 2686 int* P0input_interp = NULL; 2687 int numP1inputs = -1; 2688 IssmDouble* P1inputs = NULL; 2689 Vector<IssmDouble>* vP1inputs = NULL; 2690 int* P1input_enums = NULL; 2691 int* P1input_interp = NULL; 2692 Vector<IssmDouble>* input_interpolations = NULL; 2693 IssmDouble* input_interpolations_serial = NULL; 2694 int* pos = NULL; 2695 IssmDouble value = 0; 2665 2696 2666 2697 /*Figure out how many inputs we have and their respective interpolation*/ 2667 Vector<IssmDouble>*input_interpolations=new Vector<IssmDouble>(maxinputs);2698 input_interpolations=new Vector<IssmDouble>(maxinputs); 2668 2699 if(this->elements->Size()){ 2669 2700 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(0)); 2670 2701 element->GetInputsInterpolations(input_interpolations); 2671 2702 } 2703 2704 /*Assemble and serialize*/ 2672 2705 input_interpolations->Assemble(); 2673 2674 /*Serialize and set output*/ 2675 IssmDouble* input_interpolations_serial = input_interpolations->ToMPISerial(); 2676 delete input_interpolations; 2706 input_interpolations_serial = input_interpolations->ToMPISerial(); 2677 2707 2678 2708 /*Count and get enums of all inputs in old mesh*/ 2679 int numP0inputs;2680 int numP1inputs;2681 int *P0input_enums = NULL;2682 int *P1input_enums = NULL;2683 int *P0input_interp = NULL;2684 int *P1input_interp = NULL;2685 2709 for(int step=0;step<2;step++){ 2686 2710 if(step){ 2687 2711 P0input_enums = xNew<int>(numP0inputs); 2712 P0input_interp = xNew<int>(numP0inputs); 2688 2713 P1input_enums = xNew<int>(numP1inputs); 2689 P0input_interp = xNew<int>(numP0inputs);2690 2714 P1input_interp = xNew<int>(numP1inputs); 2691 2715 } … … 2720 2744 } 2721 2745 } 2722 2723 /*========== Deal with P0 inputs ==========*/ 2724 int numelementsold = this->elements->NumberOfElements(); 2725 int numelementsnew = newfemmodel_elements->NumberOfElements(); 2726 int numverticesold = this->vertices->NumberOfVertices(); 2727 int numverticesnew = newfemmodel_vertices->NumberOfVertices(); 2728 IssmDouble* P0inputsold = xNew<IssmDouble>(numelementsold*numP0inputs); 2729 IssmDouble* P0inputsnew = NULL; 2730 IssmDouble* vector = NULL; 2731 2732 for(int i=0;i<numP0inputs;i++){ 2733 GetVectorFromInputsx(&vector,this,P0input_enums[i],ElementSIdEnum); 2734 2735 /*Copy vector in matrix*/ 2736 for(int j=0;j<numelementsold;j++) P0inputsold[j*numP0inputs+i] = vector[j]; 2737 xDelete<IssmDouble>(vector); 2738 } 2739 2740 /*========== Deal with P1 inputs ==========*/ 2741 IssmDouble* P1inputsold = xNew<IssmDouble>(numverticesold*numP1inputs); 2742 IssmDouble* P1inputsnew = NULL; 2743 vector = NULL; 2744 2745 for(int i=0;i<numP1inputs;i++){ 2746 GetVectorFromInputsx(&vector,this,P1input_enums[i],VertexSIdEnum); 2747 2748 /*Copy vector in matrix*/ 2749 for(int j=0;j<numverticesold;j++) P1inputsold[j*numP1inputs+i] = vector[j]; 2750 xDelete<IssmDouble>(vector); 2751 } 2752 2753 /*Old mesh coordinates*/ 2754 IssmDouble *Xold = NULL; 2755 IssmDouble *Yold = NULL; 2756 IssmDouble *Zold = NULL; 2757 int *Indexold = NULL; 2758 IssmDouble *Xnew = NULL; 2759 IssmDouble *Ynew = NULL; 2760 IssmDouble *Znew = NULL; 2761 IssmDouble* XC_new = NULL; 2762 IssmDouble* YC_new = NULL; 2763 int *Indexnew = NULL; 2764 2765 /*Get the old mesh*/ 2766 this->GetMesh(this->vertices,this->elements,&Xold,&Yold,&Zold,&Indexold); 2767 2768 /*Get the new mesh*/ 2769 this->GetMesh(newfemmodel_vertices,newfemmodel_elements,&Xnew,&Ynew,&Znew,&Indexnew); 2770 2771 /*Calculate the center points xc and xy*/ 2772 int elementswidth = this->GetElementsWidth(); //just tria in this version 2773 /*new mesh*/ 2774 XC_new=xNewZeroInit<IssmDouble>(numelementsnew); 2775 YC_new=xNewZeroInit<IssmDouble>(numelementsnew); 2776 for(int i=0;i<numelementsnew;i++){ 2746 2747 /*Get P0 and P1 inputs over the elements*/ 2748 pos = xNew<int>(elementswidth); 2749 vP0inputs= new Vector<IssmDouble>(numberofelements*numP0inputs); 2750 vP1inputs= new Vector<IssmDouble>(numberofvertices*numP1inputs); 2751 for(int i=0;i<this->elements->Size();i++){ 2752 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); 2753 2754 /*Get P0 inputs*/ 2755 for(int j=0;j<numP0inputs;j++){ 2756 TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j])); 2757 input->GetInputAverage(&value); 2758 pos[0]=element->Sid()*numP0inputs+j; 2759 /*Insert input in the vector*/ 2760 vP0inputs->SetValues(1,pos,&value,INS_VAL); 2761 } 2762 2763 /*Get P1 inputs*/ 2764 for(int j=0;j<numP1inputs;j++){ 2765 TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P1input_enums[j])); 2766 pos[0]=element->vertices[0]->Sid()*numP1inputs+j; 2767 pos[1]=element->vertices[1]->Sid()*numP1inputs+j; 2768 pos[2]=element->vertices[2]->Sid()*numP1inputs+j; 2769 /*Insert input in the vector*/ 2770 vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL); 2771 } 2772 } 2773 2774 /*Assemble and serialize*/ 2775 vP0inputs->Assemble(); 2776 vP1inputs->Assemble(); 2777 P0inputs=vP0inputs->ToMPISerial(); 2778 P1inputs=vP1inputs->ToMPISerial(); 2779 2780 /*Assign pointers*/ 2781 *pnumP0inputs = numP0inputs; 2782 *pP0inputs = P0inputs; 2783 *pP0input_enums = P0input_enums; 2784 *pP0input_interp = P0input_interp; 2785 *pnumP1inputs = numP1inputs; 2786 *pP1inputs = P1inputs; 2787 *pP1input_enums = P1input_enums; 2788 *pP1input_interp = P1input_interp; 2789 2790 /*Cleanup*/ 2791 delete input_interpolations; 2792 delete vP0inputs; 2793 delete vP1inputs; 2794 xDelete<IssmDouble>(input_interpolations_serial); 2795 xDelete<int>(pos); 2796 } 2797 /*}}}*/ 2798 void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/ 2799 2800 int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh 2801 int newnumberofelements = newfemmodel_elements->Size(); //just on the new partition 2802 int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh 2803 int newnumberofvertices = newfemmodel_vertices->Size(); //just on the new partition 2804 int elementswidth = this->GetElementsWidth(); //just tria in this version 2805 int numP0inputs = -1; 2806 IssmDouble* P0inputs = NULL; //global, entire old mesh 2807 IssmDouble* newP0inputs = NULL; //just on the new partition 2808 int* P0input_enums = NULL; 2809 int* P0input_interp = NULL; 2810 int numP1inputs = -1; 2811 IssmDouble* P1inputs = NULL; //global, entire old mesh 2812 IssmDouble* newP1inputs = NULL; //just on the new partition 2813 int* P1input_enums = NULL; 2814 int* P1input_interp = NULL; 2815 IssmDouble* values = NULL; 2816 IssmDouble* vector = NULL; 2817 IssmDouble* x = NULL;//global, entire old mesh 2818 IssmDouble* y = NULL;//global, entire old mesh 2819 IssmDouble* z = NULL;//global, entire old mesh 2820 int* elementslist = NULL;//global, entire old mesh 2821 IssmDouble* newx = NULL;//just on the new partition 2822 IssmDouble* newy = NULL;//just on the new partition 2823 IssmDouble* newz = NULL;//just on the new partition 2824 IssmDouble* newxc = NULL;//just on the new partition 2825 IssmDouble* newyc = NULL;//just on the new partition 2826 int* newelementslist = NULL;//just on the new partition 2827 int* sidtoindex = NULL;//global vertices sid to partition index 2828 2829 /*Get old P0 and P1 inputs (entire mesh)*/ 2830 this->GetInputs(&numP0inputs,&P0inputs,&P0input_enums,&P0input_interp,&numP1inputs,&P1inputs,&P1input_enums,&P1input_interp); 2831 2832 /*Get the old mesh (global, entire mesh)*/ 2833 this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist); 2834 2835 /*Get the new mesh (just on the new partition)*/ 2836 this->GetMeshOnPartition(newfemmodel_vertices,newfemmodel_elements,&newx,&newy,&newz,&newelementslist,&sidtoindex); 2837 2838 /*Calculate the center points xc and xy (new mesh, new partition)*/ 2839 newxc=xNewZeroInit<IssmDouble>(newnumberofelements); 2840 newyc=xNewZeroInit<IssmDouble>(newnumberofelements); 2841 for(int i=0;i<newnumberofelements;i++){ 2777 2842 for(int j=0;j<elementswidth;j++){ 2778 int vid = Indexnew[i*elementswidth+j]-1;//Transform to C indexing 2779 XC_new[i]+=Xnew[vid]; 2780 YC_new[i]+=Ynew[vid]; 2781 } 2782 XC_new[i]=XC_new[i]/3.; 2783 YC_new[i]=YC_new[i]/3.; 2784 } 2785 2786 /*Interplate P0 inputs in the new mesh*/ 2787 InterpFromMeshToMesh2dx(&P0inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold, 2788 P0inputsold,numelementsold,numP0inputs, 2789 XC_new,YC_new,numelementsnew,NULL); 2790 2791 /*Interpolate P1 inputs in the new mesh*/ 2792 InterpFromMeshToMesh2dx(&P1inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold, 2793 P1inputsold,numverticesold,numP1inputs, 2794 Xnew,Ynew,numverticesnew,NULL); 2795 2796 /*Insert P0 inputs into the new elements.*/ 2797 vector=NULL; 2798 for(int i=0;i<numP0inputs;i++){ 2799 2800 /*Get P0 input vector from the interpolated matrix*/ 2801 vector=xNew<IssmDouble>(numelementsnew); 2802 for(int j=0;j<numelementsnew;j++) vector[j]=P0inputsnew[j*numP0inputs+i];//vector has values for all elements (serial) 2803 2804 /*Update elements from inputs: */ 2805 for(int j=0;j<newfemmodel_elements->Size();j++){ 2806 Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(j)); 2807 switch(P0input_interp[i]){ 2843 int vid = newelementslist[i*elementswidth+j]-1;//Transform to C indexing 2844 newxc[i]+=newx[vid]/elementswidth; 2845 newyc[i]+=newy[vid]/elementswidth; 2846 } 2847 } 2848 2849 /*Interplate P0 inputs in the new mesh (just on the new partition)*/ 2850 InterpFromMeshToMesh2dx(&newP0inputs,elementslist,x,y,numberofvertices,numberofelements, 2851 P0inputs,numberofelements,numP0inputs, 2852 newxc,newyc,newnumberofelements,NULL); 2853 2854 /*Interpolate P1 inputs in the new mesh (just on the new partition)*/ 2855 InterpFromMeshToMesh2dx(&newP1inputs,elementslist,x,y,numberofvertices,numberofelements, 2856 P1inputs,numberofvertices,numP1inputs, 2857 newx,newy,newnumberofvertices,NULL); 2858 2859 /*Insert P0 and P1 inputs into the new elements (just on the new partition)*/ 2860 values=xNew<IssmDouble>(elementswidth); 2861 for(int i=0;i<newfemmodel_elements->Size();i++){//just on the new partition 2862 Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(i)); 2863 /*newP0inputs is just on the new partition*/ 2864 for(int j=0;j<numP0inputs;j++){ 2865 switch(P0input_interp[j]){ 2808 2866 case P0Enum: 2809 2867 case DoubleInputEnum: 2810 element->AddInput(new DoubleInput(P0input_enums[ i],vector[element->sid]));//sid because newfemmodel has just a partitioning2868 element->AddInput(new DoubleInput(P0input_enums[j],newP0inputs[i*numP0inputs+j])); 2811 2869 break; 2812 2870 case IntInputEnum: 2813 element->AddInput(new IntInput(P0input_enums[ i],reCast<int>(vector[element->sid])));//sid because newfemmodel has just a partitioning2871 element->AddInput(new IntInput(P0input_enums[j],reCast<int>(newP0inputs[i*numP0inputs+j]))); 2814 2872 break; 2815 2873 case BoolInputEnum: 2816 element->AddInput(new BoolInput(P0input_enums[ i],reCast<bool>(vector[element->sid])));//sid because newfemmodel has just a partitioning2874 element->AddInput(new BoolInput(P0input_enums[j],reCast<bool>(newP0inputs[i*numP0inputs+j]))); 2817 2875 break; 2818 2876 default: 2819 _error_(EnumToStringx(P0input_enums[ i])<<" Not supported yet");2877 _error_(EnumToStringx(P0input_enums[j])<<" Not supported yet"); 2820 2878 } 2821 2879 } 2822 2823 xDelete<IssmDouble>(vector); 2824 } 2825 2826 /*Insert P1 inputs into the new elements.*/ 2827 vector=NULL; 2828 for(int i=0;i<numP1inputs;i++){ 2829 2830 /*Get P1 input vector from the interpolated matrix*/ 2831 vector=xNew<IssmDouble>(numverticesnew); 2832 for(int j=0;j<numverticesnew;j++) vector[j]=P1inputsnew[j*numP1inputs+i];//vector has all vertices (serial) 2833 2834 /*Update elements from inputs: */ 2835 //InputUpdateFromVectorx(newfemmodel,vector,P1input_enums[i],VertexSIdEnum);//VertexSId because vector is serial in SId indexing 2836 for(int j=0;j<newfemmodel_elements->Size();j++){ 2837 Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(j)); 2838 element->InputUpdateFromVector(vector,P1input_enums[i],VertexSIdEnum); 2839 } 2840 2841 xDelete<IssmDouble>(vector); 2842 } 2843 2880 /*newP1inputs is just on the new partition*/ 2881 for(int j=0;j<numP1inputs;j++){ 2882 values[0]=newP1inputs[sidtoindex[element->vertices[0]->Sid()]*numP1inputs+j]; 2883 values[1]=newP1inputs[sidtoindex[element->vertices[1]->Sid()]*numP1inputs+j]; 2884 values[2]=newP1inputs[sidtoindex[element->vertices[2]->Sid()]*numP1inputs+j]; 2885 element->inputs->AddInput(new TriaInput(P1input_enums[j],values,P1Enum)); 2886 } 2887 } 2888 2844 2889 /*Cleanup*/ 2845 xDelete<IssmDouble>(input_interpolations_serial); 2846 xDelete<IssmDouble>(P0inputsold); 2847 xDelete<IssmDouble>(P0inputsnew); 2890 xDelete<IssmDouble>(P0inputs); 2891 xDelete<IssmDouble>(newP0inputs); 2848 2892 xDelete<int>(P0input_enums); 2849 2893 xDelete<int>(P0input_interp); 2850 xDelete<IssmDouble>(P1inputs old);2851 xDelete<IssmDouble>( P1inputsnew);2894 xDelete<IssmDouble>(P1inputs); 2895 xDelete<IssmDouble>(newP1inputs); 2852 2896 xDelete<int>(P1input_enums); 2853 2897 xDelete<int>(P1input_interp); 2854 xDelete<IssmDouble>(Xold); 2855 xDelete<IssmDouble>(Yold); 2856 xDelete<IssmDouble>(Zold); 2857 xDelete<int>(Indexold); 2858 xDelete<IssmDouble>(Xnew); 2859 xDelete<IssmDouble>(Ynew); 2860 xDelete<IssmDouble>(Znew); 2861 xDelete<IssmDouble>(XC_new); 2862 xDelete<IssmDouble>(YC_new); 2863 xDelete<int>(Indexnew); 2898 xDelete<IssmDouble>(x); 2899 xDelete<IssmDouble>(y); 2900 xDelete<IssmDouble>(z); 2901 xDelete<int>(elementslist); 2902 xDelete<IssmDouble>(newx); 2903 xDelete<IssmDouble>(newy); 2904 xDelete<IssmDouble>(newz); 2905 xDelete<IssmDouble>(newxc); 2906 xDelete<IssmDouble>(newyc); 2907 xDelete<int>(newelementslist); 2908 xDelete<int>(sidtoindex); 2909 xDelete<IssmDouble>(values); 2864 2910 } 2865 2911 /*}}}*/ … … 3109 3155 } 3110 3156 } 3111 return;3112 3157 } 3113 3158 /*}}}*/ … … 3117 3162 if(!femmodel_elements) _error_("GetMesh: elements are NULL."); 3118 3163 3119 int numberofvertices, numberofelements; 3120 int elementswidth = this->GetElementsWidth(); // just 2D mesh in this version (just tria elements) 3121 IssmDouble *x = NULL; 3122 IssmDouble *y = NULL; 3123 IssmDouble *z = NULL; 3124 int* elementslist = NULL; 3125 3164 int numberofvertices = femmodel_vertices->NumberOfVertices(); 3165 int numberofelements = femmodel_elements->NumberOfElements(); 3166 int elementswidth = this->GetElementsWidth(); // just 2D mesh in this version (just tria elements) 3167 IssmDouble* x = NULL; 3168 IssmDouble* y = NULL; 3169 IssmDouble* z = NULL; 3170 int* elementslist = NULL; 3171 int* elem_vertices = NULL; 3172 IssmDouble *id1 = NULL; 3173 IssmDouble *id2 = NULL; 3174 IssmDouble *id3 = NULL; 3175 3126 3176 /*Get vertices coordinates*/ 3127 VertexCoordinatesx(&x, &y, &z, femmodel_vertices,false) ; 3128 3129 numberofvertices = femmodel_vertices->NumberOfVertices(); 3130 numberofelements = femmodel_elements->NumberOfElements(); 3177 VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ; 3131 3178 3132 3179 /*Get element vertices*/ 3133 int* elem_vertices=xNew<int>(elementswidth);3180 elem_vertices = xNew<int>(elementswidth); 3134 3181 Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements); 3135 3182 Vector<IssmDouble>* vid2= new Vector<IssmDouble>(numberofelements); … … 3151 3198 3152 3199 /*Serialize*/ 3153 IssmDouble *id1 = vid1->ToMPISerial();3154 IssmDouble *id2 = vid2->ToMPISerial();3155 IssmDouble *id3 = vid3->ToMPISerial();3200 id1 = vid1->ToMPISerial(); 3201 id2 = vid2->ToMPISerial(); 3202 id3 = vid3->ToMPISerial(); 3156 3203 3157 3204 /*Construct elements list*/ … … 3180 3227 } 3181 3228 /*}}}*/ 3182 void FemModel::CreateConstraints(int newnumberofvertices,int newnumberofelements,int nodecounter,int constraintcounter,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints){/*{{{*/ 3183 3184 /*itapopo ATTENTION: JUST SPCVX AND SPCVY TO TEST!!!*/ 3185 /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED!!!*/ 3186 3187 /*Get x and y of the mesh i-1*/ 3188 int numberofvertices = this->vertices->NumberOfVertices(); 3189 int numberofelements = this->elements->NumberOfElements(); 3190 IssmDouble *x = NULL; 3191 IssmDouble *y = NULL; 3192 IssmDouble *z = NULL; 3193 int *elementslist = NULL; 3194 3195 /*elementslist is in Matlab indexing*/ 3229 void FemModel::GetMeshOnPartition(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist,int** psidtoindex){/*{{{*/ 3230 3231 if(!femmodel_vertices) _error_("GetMesh: vertices are NULL."); 3232 if(!femmodel_elements) _error_("GetMesh: elements are NULL."); 3233 3234 int numberofvertices = femmodel_vertices->Size(); //number of vertices of this partition 3235 int numbertotalofvertices = femmodel_vertices->NumberOfVertices(); //number total of vertices (entire mesh) 3236 int numberofelements = femmodel_elements->Size(); //number of elements of this partition 3237 int elementswidth = this->GetElementsWidth(); //just 2D mesh in this version (just tria elements) 3238 IssmDouble* x = NULL; 3239 IssmDouble* y = NULL; 3240 IssmDouble* z = NULL; 3241 int* elementslist = NULL; 3242 int* sidtoindex = NULL; 3243 int* elem_vertices = NULL; 3244 3245 /*Get vertices coordinates of this partition*/ 3246 sidtoindex = xNewZeroInit<int>(numbertotalofvertices);//entire mesh, all vertices 3247 x = xNew<IssmDouble>(numberofvertices);//just this partition 3248 y = xNew<IssmDouble>(numberofvertices);//just this partitio; 3249 z = xNew<IssmDouble>(numberofvertices);//just this partitio; 3250 3251 /*Go through in this partition (vertices)*/ 3252 for(int i=0;i<numberofvertices;i++){//just this partition 3253 Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i); 3254 /*Attention: no spherical coordinates*/ 3255 x[i]=vertex->GetX(); 3256 y[i]=vertex->GetY(); 3257 z[i]=vertex->GetZ(); 3258 /*Keep the index and sid pair*/ 3259 sidtoindex[vertex->Sid()]=i; 3260 } 3261 3262 /*Go through in this partition (elements) and build the element list*/ 3263 elem_vertices= xNew<int>(elementswidth); 3264 elementslist = xNew<int>(numberofelements*elementswidth); 3265 if(numberofelements*elementswidth<0) _error_("numberofelements negative."); 3266 3267 for(int i=0;i<numberofelements;i++){//just this partition 3268 Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i)); 3269 element->GetVerticesSidList(elem_vertices); 3270 elementslist[elementswidth*i+0] = sidtoindex[elem_vertices[0]]+1; //InterpMesh wants Matlab indexing 3271 elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing 3272 elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing 3273 } 3274 3275 /*Assign pointers*/ 3276 *px = x; 3277 *py = y; 3278 *pz = z; 3279 *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type. 3280 *psidtoindex = sidtoindex; //it is ncessary to insert inputs 3281 3282 /*Cleanup*/ 3283 xDelete<int>(elem_vertices); 3284 } 3285 /*}}}*/ 3286 void FemModel::CreateConstraints(Vertices* newfemmodel_vertices,int nodecounter,int constraintcounter,int analysis_enum,Constraints* newfemmodel_constraints){/*{{{*/ 3287 3288 /*ATTENTION: JUST SPCVX AND SPCVY*/ 3289 /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/ 3290 if(analysis_enum!=StressbalanceAnalysisEnum) return; 3291 3292 int numberofnodes_analysistype= this->nodes->NumberOfNodes(analysis_enum); 3293 int dofpernode = 2; //vx and vy 3294 int numberofcols = dofpernode*2; //to keep dofs and flags in the vspc vector 3295 int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh 3296 int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh 3297 int newnumberofvertices = newfemmodel_vertices->Size(); //local, just the new partition 3298 int count = 0; 3299 IssmDouble* x = NULL; //global, entire old mesh 3300 IssmDouble* y = NULL; //global, entire old mesh 3301 IssmDouble* z = NULL; //global, entire old mesh 3302 int* elementslist = NULL; //global, entire old mesh 3303 IssmDouble* spc = NULL; //global, entire old mesh 3304 IssmDouble* newx = NULL; //local, just new partition 3305 IssmDouble* newy = NULL; //local, just new partition 3306 IssmDouble* newspc = NULL; //local, just new partition 3307 IssmDouble eps = 1.e-8; 3308 Vector<IssmDouble>* vspc = new Vector<IssmDouble>(numberofnodes_analysistype*numberofcols); 3309 3310 /*Get old mesh (global, entire mesh). Elementslist comes in Matlab indexing*/ 3196 3311 this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist); 3197 3198 /*Get spcvx and spcvy for mesh i-1*/3199 IssmDouble *spcvx = NULL;3200 IssmDouble *spcvy = NULL;3201 IssmDouble *spcvxflag = NULL;3202 IssmDouble *spcvyflag = NULL;3203 int numberofnodes_analysistype = this->nodes->NumberOfNodes(StressbalanceAnalysisEnum);3204 Vector<IssmDouble>* vspcvx = new Vector<IssmDouble>(numberofnodes_analysistype);3205 Vector<IssmDouble>* vspcvy = new Vector<IssmDouble>(numberofnodes_analysistype);3206 Vector<IssmDouble>* vspcvxflag = new Vector<IssmDouble>(numberofnodes_analysistype);3207 Vector<IssmDouble>* vspcvyflag = new Vector<IssmDouble>(numberofnodes_analysistype); 3208 3312 3313 /*Get vertices coordinates of the new partition*/ 3314 newx=xNew<IssmDouble>(newnumberofvertices);//just the new partition 3315 newy=xNew<IssmDouble>(newnumberofvertices);//just the new partition 3316 for(int i=0;i<newnumberofvertices;i++){//just the new partition 3317 Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i); 3318 /*Attention: no spherical coordinates*/ 3319 newx[i]=vertex->GetX(); 3320 newy[i]=vertex->GetY(); 3321 } 3322 3323 /*Get spcvx and spcvy of old mesh*/ 3209 3324 for(int i=0;i<this->constraints->Size();i++){ 3210 SpcStatic* spc = xDynamicCast<SpcStatic*>(this->constraints->GetObjectByOffset(i)); 3211 int dof = spc->GetDof(); 3212 int node = spc->GetNodeId(); 3213 IssmDouble spcvalue = spc->GetValue(); 3325 3326 Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i); 3327 if(!constraint->InAnalysis(analysis_enum)) _error_("AMR create constraints for "<<EnumToStringx(analysis_enum)<<" not supported yet!\n"); 3328 3329 SpcStatic* spcstatic = xDynamicCast<SpcStatic*>(constraint); 3330 int dof = spcstatic->GetDof(); 3331 int node = spcstatic->GetNodeId(); 3332 IssmDouble spcvalue = spcstatic->GetValue(); 3214 3333 int nodeindex = node-1; 3215 3334 3335 /*vx and vx flag insertion*/ 3216 3336 if(dof==0) {//vx 3217 vspc vx->SetValue(nodeindex,spcvalue,INS_VAL);3218 vspc vxflag->SetValue(nodeindex,1,INS_VAL);3337 vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL); //vx 3338 vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag 3219 3339 } 3340 /*vy and vy flag insertion*/ 3220 3341 if(dof==1){//vy 3221 vspcvy->SetValue(nodeindex,spcvalue,INS_VAL); 3222 vspcvyflag->SetValue(nodeindex,1,INS_VAL); 3223 } 3224 } 3225 3226 /*Assemble*/ 3227 vspcvx->Assemble(); 3228 vspcvy->Assemble(); 3229 vspcvxflag->Assemble(); 3230 vspcvyflag->Assemble(); 3231 3232 /*Serialize*/ 3233 spcvx = vspcvx->ToMPISerial(); 3234 spcvy = vspcvy->ToMPISerial(); 3235 spcvxflag = vspcvxflag->ToMPISerial(); 3236 spcvyflag = vspcvyflag->ToMPISerial(); 3237 3238 IssmDouble *newspcvx = NULL; 3239 IssmDouble *newspcvy = NULL; 3240 IssmDouble *newspcvxflag = NULL; 3241 IssmDouble *newspcvyflag = NULL; 3242 int nods_data = numberofvertices; 3243 int nels_data = numberofelements; 3244 int M_data = numberofvertices; 3245 int N_data = 1; 3246 int N_interp = newnumberofvertices; 3247 3248 /*Interpolate spcvx and spcvy in the new mesh*/ 3249 InterpFromMeshToMesh2dx(&newspcvx,elementslist,x,y,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,NULL); 3250 InterpFromMeshToMesh2dx(&newspcvy,elementslist,x,y,nods_data,nels_data,spcvy,M_data,N_data,newx,newy,N_interp,NULL); 3251 InterpFromMeshToMesh2dx(&newspcvxflag,elementslist,x,y,nods_data,nels_data,spcvxflag,M_data,N_data,newx,newy,N_interp,NULL); 3252 InterpFromMeshToMesh2dx(&newspcvyflag,elementslist,x,y,nods_data,nels_data,spcvyflag,M_data,N_data,newx,newy,N_interp,NULL); 3253 3254 int count = 0; 3255 IssmDouble eps = 1.e-8; 3256 3342 vspc->SetValue(nodeindex*numberofcols+1,spcvalue,INS_VAL); //vy 3343 vspc->SetValue(nodeindex*numberofcols+dofpernode+1,1,INS_VAL);//vyflag 3344 } 3345 } 3346 3347 /*Assemble and serialize*/ 3348 vspc->Assemble(); 3349 spc=vspc->ToMPISerial(); 3350 3351 /*Interpolate spc values and flags in the new partition*/ 3352 InterpFromMeshToMesh2dx(&newspc,elementslist,x,y,numberofvertices,numberofelements, 3353 spc,numberofvertices,numberofcols, 3354 newx,newy,newnumberofvertices,NULL); 3355 3257 3356 /*Now, insert the interpolated constraints in the data set (constraints)*/ 3258 for(int i=0;i<newnumberofvertices;i++){3259 if(my_vertices[i]){3260 /*spcvx*/3261 if(!xIsNan<IssmDouble>(newspcvx[i]) && newspcvxflag[i]>(1-eps)){3262 constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,newspcvx[i],StressbalanceAnalysisEnum));3263 //add count'th spc, on node i+1, setting dof 1 to vx.3264 count++;3265 }3266 }3267 }3268 3357 count=0; 3269 for(int i=0;i<newnumberofvertices;i++){ 3270 if(my_vertices[i]){ 3271 /*spcvy*/ 3272 if(!xIsNan<IssmDouble>(newspcvy[i]) && newspcvyflag[i]>(1-eps) ){ 3273 constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,newspcvy[i],StressbalanceAnalysisEnum)); 3274 //add count'th spc, on node i+1, setting dof 1 to vx. 3275 count++; 3276 } 3358 for(int i=0;i<newnumberofvertices;i++){//just in the new partition 3359 Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i); 3360 /*spcvx*/ 3361 if(!xIsNan<IssmDouble>(newspc[i*numberofcols]) && newspc[i*numberofcols+dofpernode]>(1-eps)){ 3362 newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,0,newspc[i*numberofcols],analysis_enum)); 3363 //add count'th spc, on node i+1, setting dof 1 to vx. 3364 count++; 3365 } 3366 } 3367 count=0; 3368 for(int i=0;i<newnumberofvertices;i++){//just in the new partition 3369 Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i); 3370 /*spcvy*/ 3371 if(!xIsNan<IssmDouble>(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){ 3372 newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); 3373 //add count'th spc, on node i+1, setting dof 1 to vx. 3374 count++; 3277 3375 } 3278 3376 } … … 3283 3381 xDelete<IssmDouble>(z); 3284 3382 xDelete<int>(elementslist); 3285 xDelete<IssmDouble>(spcvx); 3286 xDelete<IssmDouble>(spcvy); 3287 xDelete<IssmDouble>(spcvxflag); 3288 xDelete<IssmDouble>(spcvyflag); 3289 xDelete<IssmDouble>(newspcvx); 3290 xDelete<IssmDouble>(newspcvy); 3291 xDelete<IssmDouble>(newspcvxflag); 3292 xDelete<IssmDouble>(newspcvyflag); 3293 delete vspcvx; 3294 delete vspcvy; 3295 delete vspcvxflag; 3296 delete vspcvyflag; 3383 xDelete<IssmDouble>(spc); 3384 xDelete<IssmDouble>(newspc); 3385 xDelete<IssmDouble>(newx); 3386 xDelete<IssmDouble>(newy); 3387 delete vspc; 3297 3388 } 3298 3389 /*}}}*/ -
issm/trunk-jpl/src/c/classes/FemModel.h
r22349 r22360 167 167 void AdjustBaseThicknessAndMask(void); 168 168 void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist); 169 void GetMeshOnPartition(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist,int** psidtoindex); 169 170 void GetGroundediceLevelSet(IssmDouble** pmasklevelset); 170 171 void CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices); … … 172 173 void CreateMaterials(int newnumberofelements,bool* my_elements,Materials* materials); 173 174 void CreateNodes(int newnumberofvertices,int* my_vertices,int nodecounter,int analysis_enum,Nodes* nodes); 174 void CreateConstraints(int newnumberofvertices,int newnumberofelements,int nodecounter,int constraintcounter,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints); 175 void CreateConstraints(Vertices* newfemmodel_vertices,int nodecounter,int constraintcounter,int analysis_enum,Constraints* newfemmodel_constraints); 176 void GetInputs(int* pnumP0inputs,IssmDouble** pP0inputs,int** pP0input_enums,int** pP0input_interp,int* pnumP1inputs,IssmDouble** pP1inputs,int** pP1input_enums,int** pP1input_interp); 175 177 void InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements); 176 178 void UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements);
Note:
See TracChangeset
for help on using the changeset viewer.
