Changeset 3588
- Timestamp:
- 04/21/10 09:27:00 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 7 added
- 19 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp
r3446 r3588 18 18 extern int my_rank; 19 19 20 //_printf_(" Configuring elements...\n");20 _printf_(" Configuring elements...\n"); 21 21 elements->Configure(elements,loads,nodes,vertices,materials,parameters); 22 //_printf_(" Configuring loads...\n");22 _printf_(" Configuring loads...\n"); 23 23 loads->Configure(elements,loads,nodes,vertices,materials,parameters); 24 //_printf_(" Configuring nodes...\n");24 _printf_(" Configuring nodes...\n"); 25 25 nodes->Configure(elements,loads,nodes,vertices,materials,parameters); 26 //_printf_(" Configuring parameters...\n");26 _printf_(" Configuring parameters...\n"); 27 27 parameters->Configure(elements,loads, nodes,vertices, materials,parameters); 28 28 -
issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
r3567 r3588 48 48 Prognostic2AnalysisEnum, 49 49 BalancedthicknessAnalysisEnum, 50 Balancedthickness2AnalysisEnum, 50 51 BalancedvelocitiesAnalysisEnum, 51 52 //melting -
issm/trunk/src/c/Makefile.am
r3529 r3588 245 245 ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\ 246 246 ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\ 247 ./ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp\ 248 ./ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp\ 249 ./ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp\ 250 ./ModelProcessorx/Balancedthickness2/CreateParametersBalancedthickness2.cpp\ 247 251 ./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\ 248 252 ./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\ … … 645 649 ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\ 646 650 ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\ 651 ./ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp\ 652 ./ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp\ 653 ./ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp\ 654 ./ModelProcessorx/Balancedthickness2/CreateParametersBalancedthickness2.cpp\ 647 655 ./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\ 648 656 ./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\ … … 761 769 ./parallel/prognostic2_core.cpp\ 762 770 ./parallel/balancedthickness_core.cpp\ 771 ./parallel/balancedthickness2_core.cpp\ 763 772 ./parallel/balancedvelocities_core.cpp\ 764 773 ./parallel/slopecompute_core.cpp\ … … 835 844 bin_PROGRAMS = 836 845 else 837 bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe prognostic2.exe balancedthickness.exe balanced velocities.exe transient.exe steadystate.exe slopecompute.exe846 bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe prognostic2.exe balancedthickness.exe balancedthickness2.exe balancedvelocities.exe transient.exe steadystate.exe slopecompute.exe 838 847 endif 839 848 … … 858 867 balancedthickness_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 859 868 869 balancedthickness2_exe_SOURCES = parallel/balancedthickness2.cpp 870 balancedthickness2_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 871 860 872 balancedvelocities_exe_SOURCES = parallel/balancedvelocities.cpp 861 873 balancedvelocities_exe_CXXFLAGS= -fPIC -D_PARALLEL_ -
issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp
r3582 r3588 11 11 12 12 void CreateConstraintsBalancedthickness(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){ 13 13 14 int i; 15 int count=0; 14 16 DataSet* constraints = NULL; 15 17 16 18 /*Create constraints: */ 17 19 constraints = new DataSet(ConstraintsEnum); 20 21 /*Fetch data: */ 22 IoModelFetchData(&iomodel->spcthickness,NULL,NULL,iomodel_handle,"spcthickness"); 23 24 count=1; //matlab indexing 25 /*Create spcs from x,y,z, as well as the spc values on those spcs: */ 26 for (i=0;i<iomodel->numberofvertices;i++){ 27 if(iomodel->my_vertices[i]){ 28 29 if ((int)iomodel->spcthickness[2*i]){ 30 31 constraints->AddObject(new Spc(count,i+1,1,*(iomodel->spcthickness+2*i+1)));//we enforce first translation degree of freedom, for temperature 32 count++; 33 } 34 } 35 } 36 37 /*Free data: */ 38 xfree((void**)&iomodel->spcthickness); 39 40 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 41 * datasets, it will not be redone: */ 42 constraints->Presort(); 43 44 cleanup_and_return: 18 45 19 46 /*Assign output pointer: */ -
issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp
r3582 r3588 8 8 #include "../../objects/objects.h" 9 9 #include "../../shared/shared.h" 10 #include "../../include/macros.h"11 10 #include "../../include/typedefs.h" 12 11 #include "../../MeshPartitionx/MeshPartitionx.h" 13 12 #include "../IoModel.h" 14 13 15 void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){14 void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ 16 15 17 16 /*Intermediary*/ 18 int i,j; 19 int vertex_index; 20 int node_index; 17 int i,j,k,n; 21 18 22 19 /*DataSets: */ 23 DataSet* elements = NULL;24 DataSet* nodes = NULL;25 DataSet* vertices = NULL;26 DataSet* materials = NULL;20 DataSet* elements = NULL; 21 DataSet* nodes = NULL; 22 DataSet* vertices = NULL; 23 DataSet* materials = NULL; 27 24 28 25 /*First create the elements, nodes and material properties: */ … … 35 32 Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); 36 33 37 /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */38 34 /*2d mesh: */ 39 35 if (strcmp(iomodel->meshtype,"2d")==0){ … … 58 54 } 59 55 }//for (i=0;i<numberofelements;i++) 60 56 61 57 /*Free data : */ 58 xfree((void**)&iomodel->elements); 62 59 xfree((void**)&iomodel->thickness); 63 60 xfree((void**)&iomodel->surface); … … 68 65 } 69 66 else{ // if (strcmp(meshtype,"2d")==0) 70 ISSMERROR("not implemented yet"); 67 68 /*Fetch data needed: */ 69 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); 70 IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); 71 IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); 72 IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); 73 IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); 74 IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); 75 IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); 76 IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); 77 78 for (i=0;i<iomodel->numberofelements;i++){ 79 if(iomodel->my_elements[i]){ 80 /*Create and add penta element to elements dataset: */ 81 elements->AddObject(new Penta(i,iomodel)); 82 83 /*Create and add material property to materials dataset: */ 84 materials->AddObject(new Matice(i,iomodel,6)); 85 } 86 }//for (i=0;i<numberofelements;i++) 87 88 /*Free data: */ 89 xfree((void**)&iomodel->elements); 90 xfree((void**)&iomodel->thickness); 91 xfree((void**)&iomodel->surface); 92 xfree((void**)&iomodel->bed); 93 xfree((void**)&iomodel->elementoniceshelf); 94 xfree((void**)&iomodel->elementonbed); 95 xfree((void**)&iomodel->elementonsurface); 96 xfree((void**)&iomodel->elementonwater); 97 71 98 } //if (strcmp(meshtype,"2d")==0) 72 99 73 /*Add new constrant material property t go materials, at the end: */100 /*Add new constrant material property to materials, at the end: */ 74 101 materials->AddObject(new Matpar(iomodel)); 75 102 76 /*Create nodes and vertices: */ 103 /*First fetch data: */ 104 if (strcmp(iomodel->meshtype,"3d")==0){ 105 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); 106 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); 107 } 77 108 IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); 78 109 IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); … … 82 113 IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed"); 83 114 IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface"); 84 IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");85 115 IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); 86 116 IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); 87 if (strcmp(iomodel->meshtype,"3d")==0){88 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");89 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");90 }91 117 92 /*Build Vertices dataset*/93 118 for (i=0;i<iomodel->numberofvertices;i++){ 94 119 … … 99 124 vertices->AddObject(new Vertex(i,iomodel)); 100 125 101 }102 }126 /*Add node to nodes dataset: */ 127 nodes->AddObject(new Node(i,iomodel)); 103 128 104 /*Build Nodes dataset -> 3 for each element (Discontinuous Galerkin)*/105 for (i=0;i<iomodel->numberofelements;i++){106 for (j=0;j<3;j++){107 108 if(iomodel->my_nodes[3*i+j]){109 110 //Get index of the vertex on which the current node is located111 vertex_index=(int)*(iomodel->elements+3*i+j)-1; //(Matlab to C indexing)112 ISSMASSERT(vertex_index>=0 && vertex_index<iomodel->numberofvertices);113 114 //Compute Node index (id-1)115 node_index=3*i+j;116 117 /*Add node to nodes dataset: */118 nodes->AddObject(new Node(vertex_index,node_index,iomodel));119 120 }121 129 } 122 130 } … … 131 139 xfree((void**)&iomodel->gridonbed); 132 140 xfree((void**)&iomodel->gridonsurface); 133 xfree((void**)&iomodel->gridonhutter);134 141 xfree((void**)&iomodel->uppernodes); 135 142 xfree((void**)&iomodel->gridonicesheet); … … 139 146 * datasets, it will not be redone: */ 140 147 elements->Presort(); 148 vertices->Presort(); 141 149 nodes->Presort(); 142 vertices->Presort();143 150 materials->Presort(); 144 151 … … 150 157 *pvertices=vertices; 151 158 *pmaterials=materials; 159 152 160 } -
issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp
r3582 r3588 8 8 #include "../../shared/shared.h" 9 9 #include "../../include/macros.h" 10 #include "../../include/typedefs.h"11 10 #include "../IoModel.h" 11 12 12 13 13 void CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){ 14 14 15 /*Intermediary*/ 16 int i,j; 17 int i1,i2; 18 int pos1,pos2; 19 double e1,e2; 20 21 /*Output*/ 22 DataSet* loads=NULL; 23 24 /*numericalflux intermediary data: */ 25 char numericalflux_type[NUMERICALFLUXSTRING]; 26 int numericalflux_id; 27 int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES]; 28 int numericalflux_elem_id; 29 double numericalflux_h[MAX_NUMERICALFLUX_NODES]; 15 DataSet* loads = NULL; 30 16 31 17 /*Create loads: */ 32 18 loads = new DataSet(LoadsEnum); 33 34 /*Get edges and elements*/ 35 IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges"); 36 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); 37 IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); 38 39 /*First load data:*/ 40 for (i=0;i<iomodel->numberofedges;i++){ 41 42 /*Get left and right elements*/ 43 e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2] 44 e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2] 45 46 /*Now, if this element is not in the partition, pass: */ 47 if(!iomodel->my_elements[(int)e1]) continue; 48 49 /*Create load*/ 50 numericalflux_id=i+1; //Matlab indexing 51 numericalflux_elem_id=(int)e1+1;//id is in matlab index 52 53 /*1: Get vertices ids*/ 54 i1=(int)iomodel->edges[4*i+0]; 55 i2=(int)iomodel->edges[4*i+1]; 56 57 if (!isnan(e2)){ 58 strcpy(numericalflux_type,"internal"); 59 60 /*Now, we must get the nodes of the 4 nodes located on the edge*/ 61 62 /*2: Get the column where these ids are located in the index*/ 63 pos1=pos2=UNDEF; 64 for(j=0;j<3;j++){ 65 if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1; 66 if (iomodel->elements[3*(int)e2+j]==i1) pos2=j+1; 67 } 68 ISSMASSERT(pos1!=UNDEF && pos2!=UNDEF); 69 70 /*3: We have the id of the elements and the position of the vertices in the index 71 * we can compute their dofs!*/ 72 numericalflux_node_ids[0]=3*(int)e1+pos1; //ex: 1 2 3 73 numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; //ex: 2 3 1 74 numericalflux_node_ids[2]=3*(int)e2+pos2; //ex: 1 2 3 75 numericalflux_node_ids[3]=3*(int)e2+((pos2+1)%3)+1; //ex: 3 1 2 76 77 numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1] -1]; 78 numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1]; 79 numericalflux_h[2]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[2]-1]-1]; 80 numericalflux_h[3]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[3]-1]-1]; 81 } 82 else{ 83 strcpy(numericalflux_type,"boundary"); 84 85 /*2: Get the column where these ids are located in the index*/ 86 pos1==UNDEF; 87 for(j=0;j<3;j++){ 88 if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1; 89 } 90 ISSMASSERT(pos1!=UNDEF); 91 92 /*3: We have the id of the elements and the position of the vertices in the index 93 * we can compute their dofs!*/ 94 numericalflux_node_ids[0]=3*(int)e1+pos1; 95 numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; 96 97 numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1]-1]; 98 numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1]; 99 numericalflux_h[2]=UNDEF; 100 numericalflux_h[3]=UNDEF; 101 } 102 103 loads->AddObject(new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_id,numericalflux_h)); 104 } 105 106 /*Free data: */ 107 xfree((void**)&iomodel->edges); 108 xfree((void**)&iomodel->elements); 109 xfree((void**)&iomodel->thickness); 110 19 20 111 21 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 112 22 * datasets, it will not be redone: */ 113 23 loads->Presort(); 24 25 cleanup_and_return: 114 26 115 27 /*Assign output pointer: */ … … 117 29 118 30 } 31 32 -
issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp
r3582 r3588 1 /*!\file: CreateParametersPrognostic.cpp 2 * \brief driver for creating parameters dataset, for prognostic analysis. 1 /* \brief driver for creating parameters dataset, for prognostic analysis. 3 2 */ 4 3 … … 16 15 int count; 17 16 int i; 18 int dim; 19 double* vx_g=NULL; 20 double* vy_g=NULL; 17 double* u_g=NULL; 21 18 22 19 /*recover parameters : */ … … 28 25 IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); 29 26 IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); 27 IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz"); 30 28 31 vx_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double)); 32 vy_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double)); 29 u_g=(double*)xcalloc(iomodel->numberofvertices*3,sizeof(double)); 33 30 34 if(iomodel->vx) for(i=0;i<iomodel->numberofvertices;i++)vx_g[i]=iomodel->vx[i]/iomodel->yts; 35 if(iomodel->vy) for(i=0;i<iomodel->numberofvertices;i++)vy_g[i]=iomodel->vy[i]/iomodel->yts; 31 if(iomodel->vx) for(i=0;i<iomodel->numberofvertices;i++)u_g[3*i+0]=iomodel->vx[i]/iomodel->yts; 32 if(iomodel->vy) for(i=0;i<iomodel->numberofvertices;i++)u_g[3*i+1]=iomodel->vy[i]/iomodel->yts; 33 if(iomodel->vz) for(i=0;i<iomodel->numberofvertices;i++)u_g[3*i+2]=iomodel->vz[i]/iomodel->yts; 36 34 37 35 count++; 38 param= new Param(count," vx_g",DOUBLEVEC);39 param->SetDoubleVec( vx_g,iomodel->numberofvertices,1);36 param= new Param(count,"u_g",DOUBLEVEC); 37 param->SetDoubleVec(u_g,3*iomodel->numberofvertices,3); 40 38 parameters->AddObject(param); 41 count++; 42 param= new Param(count,"vy_g",DOUBLEVEC); 43 param->SetDoubleVec(vy_g,iomodel->numberofvertices,1); 44 parameters->AddObject(param); 39 45 40 46 41 xfree((void**)&iomodel->vx); 47 42 xfree((void**)&iomodel->vy); 48 xfree((void**)&vx_g); 49 xfree((void**)&vy_g); 50 51 /*Get thickness: */ 52 IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); 53 54 count++; 55 param= new Param(count,"h_g",DOUBLEVEC); 56 if(iomodel->thickness) param->SetDoubleVec(iomodel->thickness,iomodel->numberofvertices,1); 57 else param->SetDoubleVec(iomodel->thickness,0,0); 58 parameters->AddObject(param); 59 60 /*Free thickness: */ 61 xfree((void**)&iomodel->thickness); 62 63 /*Get dhdt: */ 64 IoModelFetchData(&iomodel->dhdt,NULL,NULL,iomodel_handle,"dhdt"); 65 if(iomodel->dhdt) for(i=0;i<iomodel->numberofvertices;i++)iomodel->dhdt[i]=iomodel->dhdt[i]/iomodel->yts; 66 67 count++; 68 param= new Param(count,"dhdt_g",DOUBLEVEC); 69 if(iomodel->dhdt) param->SetDoubleVec(iomodel->dhdt,iomodel->numberofvertices,1); 70 else param->SetDoubleVec(iomodel->dhdt,0,1); 71 parameters->AddObject(param); 72 73 /*Free dhdt: */ 74 xfree((void**)&iomodel->dhdt); 43 xfree((void**)&iomodel->vz); 44 xfree((void**)&u_g); 75 45 76 46 /*Get melting: */ -
issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp
r3570 r3588 100 100 101 101 } 102 else if (iomodel->analysis_type==Balancedthickness2AnalysisEnum){ 103 104 CreateElementsNodesAndMaterialsBalancedthickness2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); 105 CreateConstraintsBalancedthickness2(pconstraints,iomodel,iomodel_handle); 106 CreateLoadsBalancedthickness2(ploads,iomodel,iomodel_handle); 107 CreateParametersBalancedthickness2(pparameters,iomodel,iomodel_handle); 108 109 } 102 110 else if (iomodel->analysis_type==BalancedvelocitiesAnalysisEnum){ 103 111 -
issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp
r3582 r3588 247 247 if ( 248 248 iomodel->analysis_type==Prognostic2AnalysisEnum || 249 iomodel->analysis_type==Balancedthickness AnalysisEnum249 iomodel->analysis_type==Balancedthickness2AnalysisEnum 250 250 ) 251 251 param->SetDouble(3*iomodel->numberofelements); -
issm/trunk/src/c/ModelProcessorx/IoModel.h
r3582 r3588 263 263 void CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); 264 264 265 /*balancedthickness2:*/ 266 void CreateElementsNodesAndMaterialsBalancedthickness2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); 267 void CreateConstraintsBalancedthickness2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); 268 void CreateLoadsBalancedthickness2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); 269 void CreateParametersBalancedthickness2(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); 270 265 271 /*balancedvelocities:*/ 266 272 void CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); -
issm/trunk/src/c/ModelProcessorx/Partitioning.cpp
r3582 r3588 22 22 void Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ 23 23 24 if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness AnalysisEnum)24 if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum) 25 25 DiscontinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle); 26 26 else -
issm/trunk/src/c/objects/Numericalflux.cpp
r3582 r3588 225 225 if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!"); 226 226 } 227 else if (analysis_type==Balancedthickness AnalysisEnum){227 else if (analysis_type==Balancedthickness2AnalysisEnum){ 228 228 /*No transient term is involved*/ 229 229 dt=1; … … 256 256 GetParameterValue(&vy, &vy_list[0],gauss_coord); 257 257 UdotN=vx*normal[0]+vy*normal[1]; 258 if (fabs(UdotN)<1.0e-9) printf("Edge number %i has a flux very small (u.n = %g ), which could lead to unaccurate results\n",id,UdotN); 259 //if (fabs(UdotN)>0 && fabs(UdotN)< 1.0e-7) UdotN= 1.0e-7; 260 //if (fabs(UdotN)<0 && fabs(UdotN)>-1.0e-7) UdotN=-1.0e-7; 258 261 259 262 /*Get L and B: */ … … 341 344 if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!"); 342 345 } 343 else if (analysis_type==Balancedthickness AnalysisEnum){346 else if (analysis_type==Balancedthickness2AnalysisEnum){ 344 347 /*No transient term is involved*/ 345 348 dt=1; … … 483 486 if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!"); 484 487 } 485 else if (analysis_type==Balancedthickness AnalysisEnum){488 else if (analysis_type==Balancedthickness2AnalysisEnum){ 486 489 /*No transient term is involved*/ 487 490 dt=1; -
issm/trunk/src/c/objects/Tria.cpp
r3582 r3588 77 77 /*hooks: */ 78 78 //go recover node ids, needed to initialize the node hook. 79 if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness AnalysisEnum){79 if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum){ 80 80 /*Discontinuous Galerkin*/ 81 81 tria_node_ids[0]=3*i+1; … … 598 598 CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type); 599 599 } 600 else if (analysis_type==Balancedthickness2AnalysisEnum){ 601 602 CreateKMatrixBalancedthickness2( Kgg,inputs,analysis_type,sub_analysis_type); 603 } 600 604 else if (analysis_type==BalancedvelocitiesAnalysisEnum){ 601 605 … … 611 615 /*FUNCTION Tria::CreateKMatrixBalancedthickness {{{1*/ 612 616 void Tria::CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ 617 618 /* local declarations */ 619 int i,j; 620 621 /* node data: */ 622 const int numgrids=3; 623 const int NDOF1=1; 624 const int numdof=NDOF1*numgrids; 625 double xyz_list[numgrids][3]; 626 int doflist[numdof]; 627 int numberofdofspernode; 628 629 /* gaussian points: */ 630 int num_gauss,ig; 631 double* first_gauss_area_coord = NULL; 632 double* second_gauss_area_coord = NULL; 633 double* third_gauss_area_coord = NULL; 634 double* gauss_weights = NULL; 635 double gauss_weight; 636 double gauss_l1l2l3[3]; 637 638 /* matrices: */ 639 double L[numgrids]; 640 double B[2][numgrids]; 641 double Bprime[2][numgrids]; 642 double DL[2][2]={0.0}; 643 double DLprime[2][2]={0.0}; 644 double DL_scalar; 645 double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix 646 double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. 647 double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. 648 double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. 649 650 double Jdettria; 651 652 /*input parameters for structural analysis (diagnostic): */ 653 double vxvy_list[numgrids][2]={0.0}; 654 double vx_list[numgrids]={0.0}; 655 double vy_list[numgrids]={0.0}; 656 double dvx[2]; 657 double dvy[2]; 658 double vx,vy; 659 double dvxdx,dvydy; 660 double v_gauss[2]={0.0}; 661 double K[2][2]={0.0}; 662 double KDL[2][2]={0.0}; 663 int dofs[2]={0,1}; 664 int found=0; 665 666 /*dynamic objects pointed to by hooks: */ 667 Node** nodes=NULL; 668 Matpar* matpar=NULL; 669 Matice* matice=NULL; 670 Numpar* numpar=NULL; 671 672 ParameterInputs* inputs=NULL; 673 674 /*recover pointers: */ 675 inputs=(ParameterInputs*)vinputs; 676 677 /*recover objects from hooks: */ 678 nodes=(Node**)hnodes.deliverp(); 679 matpar=(Matpar*)hmatpar.delivers(); 680 matice=(Matice*)hmatice.delivers(); 681 numpar=(Numpar*)hnumpar.delivers(); 682 683 /*recover extra inputs from users, at current convergence iteration: */ 684 found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); 685 if(!found)ISSMERROR(" could not find velocity_average in inputs!"); 686 687 for(i=0;i<numgrids;i++){ 688 vx_list[i]=vxvy_list[i][0]; 689 vy_list[i]=vxvy_list[i][1]; 690 } 691 692 /* Get node coordinates and dof list: */ 693 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 694 GetDofList(&doflist[0],&numberofdofspernode); 695 696 //Create Artificial diffusivity once for all if requested 697 if(numpar->artdiff){ 698 //Get the Jacobian determinant 699 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; 700 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); 701 702 //Build K matrix (artificial diffusivity matrix) 703 v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); 704 v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); 705 706 K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); 707 K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); 708 } 709 710 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 711 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); 712 713 /* Start looping on the number of gaussian points: */ 714 for (ig=0; ig<num_gauss; ig++){ 715 /*Pick up the gaussian point: */ 716 gauss_weight=*(gauss_weights+ig); 717 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 718 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 719 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 720 721 /* Get Jacobian determinant: */ 722 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); 723 724 /*Get B and B prime matrix: */ 725 GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); 726 GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); 727 728 //Get vx, vy and their derivatives at gauss point 729 GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); 730 GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); 731 732 GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3); 733 GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3); 734 735 dvxdx=dvx[0]; 736 dvydy=dvy[1]; 737 738 DL_scalar=gauss_weight*Jdettria; 739 740 //Create DL and DLprime matrix 741 DL[0][0]=DL_scalar*dvxdx; 742 DL[1][1]=DL_scalar*dvydy; 743 744 DLprime[0][0]=DL_scalar*vx; 745 DLprime[1][1]=DL_scalar*vy; 746 747 //Do the triple product tL*D*L. 748 //Ke_gg_thickness=B'*DLprime*Bprime; 749 750 TripleMultiply( &B[0][0],2,numdof,1, 751 &DL[0][0],2,2,0, 752 &B[0][0],2,numdof,0, 753 &Ke_gg_thickness1[0][0],0); 754 755 TripleMultiply( &B[0][0],2,numdof,1, 756 &DLprime[0][0],2,2,0, 757 &Bprime[0][0],2,numdof,0, 758 &Ke_gg_thickness2[0][0],0); 759 760 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ 761 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j]; 762 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j]; 763 764 if(numpar->artdiff){ 765 766 /* Compute artificial diffusivity */ 767 KDL[0][0]=DL_scalar*K[0][0]; 768 KDL[1][1]=DL_scalar*K[1][1]; 769 770 TripleMultiply( &Bprime[0][0],2,numdof,1, 771 &KDL[0][0],2,2,0, 772 &Bprime[0][0],2,numdof,0, 773 &Ke_gg_gaussian[0][0],0); 774 775 /* Add artificial diffusivity matrix */ 776 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 777 778 } 779 780 #ifdef _DEBUGELEMENTS_ 781 if(my_rank==RANK && id==ELID){ 782 printf(" B:\n"); 783 for(i=0;i<3;i++){ 784 for(j=0;j<numdof;j++){ 785 printf("%g ",B[i][j]); 786 } 787 printf("\n"); 788 } 789 printf(" Bprime:\n"); 790 for(i=0;i<3;i++){ 791 for(j=0;j<numdof;j++){ 792 printf("%g ",Bprime[i][j]); 793 } 794 printf("\n"); 795 } 796 } 797 #endif 798 } // for (ig=0; ig<num_gauss; ig++) 799 800 /*Add Ke_gg to global matrix Kgg: */ 801 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); 802 803 #ifdef _DEBUGELEMENTS_ 804 if(my_rank==RANK && id==ELID){ 805 printf(" Ke_gg erms:\n"); 806 for( i=0; i<numdof; i++){ 807 for (j=0;j<numdof;j++){ 808 printf("%g ",Ke_gg[i][j]); 809 } 810 printf("\n"); 811 } 812 printf(" Ke_gg row_indices:\n"); 813 for( i=0; i<numdof; i++){ 814 printf("%i ",doflist[i]); 815 } 816 817 } 818 #endif 819 820 cleanup_and_return: 821 xfree((void**)&first_gauss_area_coord); 822 xfree((void**)&second_gauss_area_coord); 823 xfree((void**)&third_gauss_area_coord); 824 xfree((void**)&gauss_weights); 825 826 } 827 /*}}}*/ 828 /*FUNCTION Tria::CreateKMatrixBalancedthickness2 {{{1*/ 829 void Tria::CreateKMatrixBalancedthickness2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ 613 830 614 831 /* local declarations */ … … 2109 2326 CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type); 2110 2327 } 2328 else if (analysis_type==Balancedthickness2AnalysisEnum){ 2329 2330 CreatePVectorBalancedthickness2( pg,inputs,analysis_type,sub_analysis_type); 2331 } 2111 2332 else if (analysis_type==BalancedvelocitiesAnalysisEnum){ 2112 2333 … … 2121 2342 /*FUNCTION Tria::CreatePVectorBalancedthickness {{{1*/ 2122 2343 void Tria::CreatePVectorBalancedthickness(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){ 2344 2345 2346 /* local declarations */ 2347 int i,j; 2348 2349 /* node data: */ 2350 const int numgrids=3; 2351 const int NDOF1=1; 2352 const int numdof=NDOF1*numgrids; 2353 double xyz_list[numgrids][3]; 2354 int doflist[numdof]; 2355 int numberofdofspernode; 2356 2357 /* gaussian points: */ 2358 int num_gauss,ig; 2359 double* first_gauss_area_coord = NULL; 2360 double* second_gauss_area_coord = NULL; 2361 double* third_gauss_area_coord = NULL; 2362 double* gauss_weights = NULL; 2363 double gauss_weight; 2364 double gauss_l1l2l3[3]; 2365 2366 /* matrix */ 2367 double pe_g[numgrids]={0.0}; 2368 double L[numgrids]; 2369 double Jdettria; 2370 2371 /*input parameters for structural analysis (diagnostic): */ 2372 double accumulation_list[numgrids]={0.0}; 2373 double accumulation_g; 2374 double melting_list[numgrids]={0.0}; 2375 double melting_g; 2376 double thickness_list[numgrids]={0.0}; 2377 double thickness_g; 2378 int dofs[1]={0}; 2379 int found=0; 2380 2381 /*dynamic objects pointed to by hooks: */ 2382 Node** nodes=NULL; 2383 Matpar* matpar=NULL; 2384 Matice* matice=NULL; 2385 Numpar* numpar=NULL; 2386 2387 ParameterInputs* inputs=NULL; 2388 2389 /*recover pointers: */ 2390 inputs=(ParameterInputs*)vinputs; 2391 2392 /*recover objects from hooks: */ 2393 nodes=(Node**)hnodes.deliverp(); 2394 matpar=(Matpar*)hmatpar.delivers(); 2395 matice=(Matice*)hmatice.delivers(); 2396 numpar=(Numpar*)hnumpar.delivers(); 2397 2398 /*recover extra inputs from users, at current convergence iteration: */ 2399 found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes); 2400 if(!found)ISSMERROR(" could not find accumulation in inputs!"); 2401 found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); 2402 if(!found)ISSMERROR(" could not find melting in inputs!"); 2403 2404 /* Get node coordinates and dof list: */ 2405 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2406 GetDofList(&doflist[0],&numberofdofspernode); 2407 2408 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 2409 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); 2410 2411 /* Start looping on the number of gaussian points: */ 2412 for (ig=0; ig<num_gauss; ig++){ 2413 /*Pick up the gaussian point: */ 2414 gauss_weight=*(gauss_weights+ig); 2415 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 2416 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 2417 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 2418 2419 /* Get Jacobian determinant: */ 2420 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); 2421 2422 /*Get L matrix: */ 2423 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); 2424 2425 /* Get accumulation, melting and thickness at gauss point */ 2426 GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3); 2427 GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3); 2428 2429 /* Add value into pe_g: */ 2430 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(accumulation_g-melting_g)*L[i]; 2431 2432 } // for (ig=0; ig<num_gauss; ig++) 2433 2434 /*Add pe_g to global matrix Kgg: */ 2435 VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); 2436 2437 cleanup_and_return: 2438 xfree((void**)&first_gauss_area_coord); 2439 xfree((void**)&second_gauss_area_coord); 2440 xfree((void**)&third_gauss_area_coord); 2441 xfree((void**)&gauss_weights); 2442 2443 } 2444 /*}}}*/ 2445 /*FUNCTION Tria::CreatePVectorBalancedthickness2 {{{1*/ 2446 void Tria::CreatePVectorBalancedthickness2(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){ 2123 2447 2124 2448 -
issm/trunk/src/c/objects/Tria.h
r3529 r3588 117 117 void CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); 118 118 void CreatePVectorBalancedthickness(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type); 119 void CreateKMatrixBalancedthickness2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); 120 void CreatePVectorBalancedthickness2(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type); 119 121 void CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); 120 122 void CreatePVectorBalancedvelocities(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type); -
issm/trunk/src/c/parallel/ProcessResults.cpp
r3567 r3588 135 135 if(analysis_type==BalancedthicknessAnalysisEnum){ 136 136 fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum); 137 } 138 if(analysis_type==Balancedthickness2AnalysisEnum){ 139 fem_p=model->GetFormulation(Balancedthickness2AnalysisEnum); 137 140 } 138 141 if(analysis_type==BalancedvelocitiesAnalysisEnum){ -
issm/trunk/src/c/parallel/balancedthickness.cpp
r3582 r3588 25 25 Model* model=NULL; 26 26 27 double* vx_g=NULL; 28 double* vy_g=NULL; 29 double* m_g=NULL; 30 double* a_g=NULL; 31 double* h_g=NULL; 32 double* dhdt_g=NULL; 27 Vec u_g=NULL; 28 double* u_g_serial=NULL; 29 double* melting_g=NULL; 30 double* accumulation_g=NULL; 33 31 double dt; 34 32 double yts; … … 62 60 MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 63 61 64 _printf_("recover input file name and output file name:\n");62 _printf_("recover , input file name and output file name:\n"); 65 63 inputfilename=argv[2]; 66 64 outputfilename=argv[3]; … … 86 84 _printf_("initialize inputs:\n"); 87 85 88 model->FindParam(&vx_g,NULL,NULL,"vx_g",BalancedthicknessAnalysisEnum); 89 model->FindParam(&vy_g,NULL,NULL,"vy_g",BalancedthicknessAnalysisEnum); 90 model->FindParam(&m_g,NULL,NULL,"m_g",BalancedthicknessAnalysisEnum); 91 model->FindParam(&a_g,NULL,NULL,"a_g",BalancedthicknessAnalysisEnum); 92 model->FindParam(&h_g,NULL,NULL,"h_g",BalancedthicknessAnalysisEnum); 93 model->FindParam(&dhdt_g,NULL,NULL,"dhdt_g",BalancedthicknessAnalysisEnum); 86 model->FindParam(&u_g_serial,NULL,NULL,"u_g",BalancedthicknessAnalysisEnum); 87 model->FindParam(&melting_g,NULL,NULL,"m_g",BalancedthicknessAnalysisEnum); 88 model->FindParam(&accumulation_g,NULL,NULL,"a_g",BalancedthicknessAnalysisEnum); 94 89 model->FindParam(&numberofnodes,"numberofnodes"); 95 90 96 91 inputs=new ParameterInputs; 97 inputs->Add("vx",vx_g,1,numberofnodes); 98 inputs->Add("vy",vy_g,1,numberofnodes); 99 inputs->Add("melting",m_g,1,numberofnodes); 100 inputs->Add("accumulation",a_g,1,numberofnodes); 101 inputs->Add("thickness",h_g,1,numberofnodes); 102 inputs->Add("dhdt",dhdt_g,1,numberofnodes); 92 inputs->Add("velocity",u_g_serial,3,numberofnodes); 93 inputs->Add("melting",melting_g,1,numberofnodes); 94 inputs->Add("accumulation",accumulation_g,1,numberofnodes); 103 95 104 96 _printf_("initialize results:\n"); … … 146 138 147 139 /*Free ressources:*/ 148 xfree((void**)&vx_g);149 xfree((void**)&vy_g);150 xfree((void**)&m_g);151 xfree((void**)&a_g);152 xfree((void**)&dhdt_g);153 xfree((void**)&h_g);154 140 delete processedresults; 155 141 delete results; 156 delete model;157 142 delete inputs; 158 143 -
issm/trunk/src/c/parallel/balancedthickness_core.cpp
r3582 r3588 19 19 20 20 /*intermediary: */ 21 Vec vx_g=NULL; 22 Vec vy_g=NULL; 21 Vec u_g=NULL; 23 22 24 23 /*solutions: */ … … 29 28 int numberofdofspernode; 30 29 int numberofnodes; 31 int numberofvertices; 32 int dofs[1]={1}; 30 int dofs[2]={1,1}; 33 31 34 32 /*fem balancedthickness model: */ … … 36 34 37 35 36 /*recover fem model: */ 38 37 fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum); 39 38 … … 41 40 model->FindParam(&verbose,"verbose"); 42 41 model->FindParam(&numberofnodes,"numberofnodes"); 43 model->FindParam(&numberofvertices,"numberofvertices");44 42 model->FindParam(&numberofdofspernode,"numberofdofspernode"); 45 43 46 44 _printf_("depth averaging velocity...\n"); 47 vx_g=inputs->Get("vx",&dofs[0],1); 48 vy_g=inputs->Get("vy",&dofs[0],1); 49 /* NOT WORKING YET.... 50 FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vx"); 51 FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vy"); 52 */ 53 inputs->Add("vx_average",vx_g,1,numberofvertices); 54 inputs->Add("vy_average",vy_g,1,numberofvertices); 45 u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity 46 FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity"); 47 inputs->Add("velocity_average",u_g,2,numberofnodes); 55 48 56 49 _printf_("call computational core:\n"); 57 50 diagnostic_core_linear(&h_g,fem_p,inputs,BalancedthicknessAnalysisEnum,NoneAnalysisEnum); 58 51 59 _printf_("Averaging over vertices:\n"); 60 FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters); 61 62 // _printf_("extrude computed thickness on all layers:\n"); 63 // FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); 52 _printf_("extrude computed thickness on all layers:\n"); 53 FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); 64 54 65 55 /*Plug results into output dataset: */ … … 68 58 69 59 /*Free ressources:*/ 70 VecFree(&vx_g); 71 VecFree(&vy_g); 60 VecFree(&u_g); 72 61 VecFree(&h_g); 73 62 } -
issm/trunk/src/c/parallel/parallel.h
r3373 r3588 19 19 void prognostic2_core(DataSet* results,Model* model, ParameterInputs* inputs); 20 20 void balancedthickness_core(DataSet* results,Model* model, ParameterInputs* inputs); 21 void balancedthickness2_core(DataSet* results,Model* model, ParameterInputs* inputs); 21 22 void balancedvelocities_core(DataSet* results,Model* model, ParameterInputs* inputs); 22 23 void slopecompute_core(DataSet* results,Model* model, ParameterInputs* inputs); -
issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp
r3570 r3588 59 59 numdofs=1; 60 60 } 61 else if (analysis_type==Balancedthickness2AnalysisEnum){ 62 numdofs=1; 63 } 61 64 else if (analysis_type==BalancedvelocitiesAnalysisEnum){ 62 65 numdofs=1;
Note:
See TracChangeset
for help on using the changeset viewer.