Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 3428) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 3429) @@ -15,113 +15,63 @@ int i,j; - int count; + int count=0; DataSet* constraints = NULL; Spc* spc = NULL; Rgb* rgb = NULL; + int node1,node2; - /*spc intermediary data: */ - int spc_sid; - int spc_node; - int spc_dof; - double spc_value; - - /*rgb constructor data: */ - int rgb_id; - int rgb_dof; - int rgb_nodeid1; - int rgb_nodeid2; - int grid1,grid2; - - double* spcvelocity=NULL; - double* gridonhutter=NULL; + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ + if (!iomodel->ismacayealpattyn)goto cleanup_and_return; /*Create constraints: */ constraints = new DataSet(ConstraintsEnum()); - /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ - if (!iomodel->ismacayealpattyn)goto cleanup_and_return; + /*Spcs: fetch data: */ + IoModelFetchData(&iomodel->spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity"); + IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter"); - /*Fetch data: */ - IoModelFetchData(&spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity"); - IoModelFetchData(&gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter"); + count=1; //matlab indexing + /*Create spcs from x,y,z, as well as the spc values on those spcs: */ + for (i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ - count=0; + if ((int)iomodel->spcvelocity[6*i+0] | (int)iomodel->gridonhutter[i]) + constraints->AddObject(new Spc(count,i+1,1,*(iomodel->spcvelocity+6*i+3)/iomodel->yts)); //add count'th spc, on node i+1, setting dof 1 to vx. - /*Create spcs from x,y,z, as well as the spc values on those spcs: */ - for (i=0;inumberofnodes;i++){ - #ifdef _PARALLEL_ - /*keep only this partition's nodes:*/ - if((iomodel->my_grids[i]==1)){ - #endif + count++; + + if ((int)iomodel->spcvelocity[6*i+1] | (int)iomodel->gridonhutter[i]) + constraints->AddObject(new Spc(count,i+1,2,*(iomodel->spcvelocity+6*i+4)/iomodel->yts)); //add count'th spc, on node i+1, setting dof 2 to vy - if ((int)spcvelocity[6*i+0] | (int)gridonhutter[i]){ - - /*This grid needs to be spc'd to vx_obs:*/ - - spc_sid=count; - spc_node=i+1; - spc_dof=1; //we enforce first x translation degree of freedom - spc_value=*(spcvelocity+6*i+3)/iomodel->yts; - - spc = new Spc(spc_sid,spc_node,spc_dof,spc_value); - constraints->AddObject(spc); count++; } + } + + /*Free data: */ + xfree((void**)&iomodel->spcvelocity); + xfree((void**)&iomodel->gridonhutter); - if ((int)spcvelocity[6*i+1] | (int)gridonhutter[i]){ - - /*This grid needs to be spc'd to vy_obs:*/ - - spc_sid=count; - spc_node=i+1; - spc_dof=2; //we enforce first y translation degree of freedom - spc_value=*(spcvelocity+6*i+4)/iomodel->yts; - - spc = new Spc(spc_sid,spc_node,spc_dof,spc_value); - constraints->AddObject(spc); - count++; - } - - #ifdef _PARALLEL_ - } //if((my_grids[i]==1)) - #endif - } - - /*Create penalties loads for linking of collapsed formulation grids with non collapsed grids: */ - - /*First fetch penalties: */ - if (strcmp(iomodel->meshtype,"3d")==0){ - IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties"); - } - - count=0; + + /*penalty loads: */ + IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties"); + + count=1; //matlab indexing if(iomodel->numpenalties){ - /*First deal with internal grids: */ for (i=0;inumpenalties;i++){ - if (iomodel->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU - for(j=1;jnumlayers;j++){ - /*We are pairing grids along a vertical profile.*/ - grid1=(int)*(iomodel->penalties+iomodel->numlayers*i); - grid2=(int)*(iomodel->penalties+iomodel->numlayers*i+j); + for(j=1;jnumlayers;j++){ + + /*We are pairing nodes along a vertical profile.*/ + node1=(int)*(iomodel->penalties+iomodel->numlayers*i); + node2=(int)*(iomodel->penalties+iomodel->numlayers*i+j); - rgb_id=count+1; //matlab indexing - rgb_dof=1; - rgb_nodeid1=grid1; - rgb_nodeid2=grid2; - rgb= new Rgb(rgb_id,rgb_nodeid1,rgb_nodeid2,rgb_dof); - constraints->AddObject(rgb); - count++; + constraints->AddObject(new Rgb(count,node1,node2,1)); //add count'th Rgb on dof 1 between node1 and node2 + + count++; + + constraints->AddObject(new Rgb(count,node1,node2,2)); //add count'th Rgb on dof 1 between node1 and node2 - rgb_id=count+1; //matlab indexing - rgb_dof=2; - rgb_nodeid1=grid1; - rgb_nodeid2=grid2; - rgb= new Rgb(rgb_id,rgb_nodeid1,rgb_nodeid2,rgb_dof); - constraints->AddObject(rgb); - count++; - } //for ( i=0; iPresort(); - /*Free data: */ - xfree((void**)&spcvelocity); - xfree((void**)&gridonhutter); - cleanup_and_return: /*Assign output pointer: */ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3428) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3429) @@ -23,4 +23,12 @@ DataSet* materials = NULL; + /*Objects: */ + Node* node = NULL; + Vertex* vertex = NULL; + Tria* tria = NULL; + Penta* penta = NULL; + Matice* matice = NULL; + Matpar* matpar = NULL; + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ if (!iomodel->ismacayealpattyn)goto cleanup_and_return; Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 3428) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 3429) @@ -13,10 +13,4 @@ void CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){ - int i; - int element; - - extern int my_rank; - extern int num_procs; - DataSet* loads = NULL; Icefront* icefront = NULL; @@ -24,50 +18,6 @@ int segment_width; - int element_type; - int i1,i2,i3,i4; - int i0,range; + int element; - - /*icefront intermediary data: */ - char icefront_type[ICEFRONTSTRING]; - int icefront_fill; - int icefront_element_type; - int icefront_sid; - int icefront_eid; - int icefront_mparid; - int icefront_node_ids[MAX_ICEFRONT_GRIDS]; - double icefront_h[MAX_ICEFRONT_GRIDS]; - double icefront_b[MAX_ICEFRONT_GRIDS]; - - /*rifts: */ - char riftfront_type[RIFTFRONTSTRING]; - int riftfront_id; - int riftfront_node_ids[2]; - int riftfront_mparid; - double riftfront_h[2]; - double riftfront_b[2]; - double riftfront_s[2]; - double riftfront_normal[2]; - double riftfront_length; - int riftfront_fill; - double riftfront_friction; - double riftfront_fraction; - double riftfront_fractionincrement; - bool riftfront_shelf; - double riftfront_penalty_offset; - int riftfront_penalty_lock; - bool riftfront_active; - bool riftfront_frozen; - int riftfront_counter; - bool riftfront_prestable; - int el1,el2; - int grid1,grid2; - double normal[2]; - double length; - int fill; - double friction; - double fraction; - double fractionincrement; - /*Create loads: */ loads = new DataSet(LoadsEnum()); @@ -86,70 +36,17 @@ for (i=0;inumberofpressureloads;i++){ - if (strcmp(iomodel->meshtype,"2d")==0){ - segment_width=4; - element_type=TriaEnum(); - } - else{ - segment_width=6; - element_type=PentaEnum(); - } - + /*Retrieve element to which this icefront belongs: */ + if (strcmp(iomodel->meshtype,"2d")==0) segment_width=4; + else segment_width=6; element=(int)(*(iomodel->pressureload+segment_width*i+segment_width-2)-1); //element is in the penultimate column (grid1 grid2 ... elem fill) - #ifdef _PARALLEL_ - if (iomodel->epart[element]!=my_rank){ - /*This load does not belong to this cluster node, as it references an element - *that does not belong to this node's partition. Drop this 'i':*/ - continue; - } + /*Now, if this element is not in the partition, pass: */ + if(!iomodel->my_elements[i])continue; + + /*Do not create ice front if Hutter or Stokes elements*/ + if ((int)*(iomodel->elements_type+2*element+0)==(HutterFormulationEnum() || StokesFormulationEnum()))continue; - #endif - - /*Do not create ice front if Hutter or Stokes elements*/ - if ((int)*(iomodel->elements_type+2*element+0)==(HutterFormulationEnum() || StokesFormulationEnum())){ - continue; - } - - icefront_mparid=iomodel->numberofelements+1; //matlab indexing - icefront_sid=i+1; //matlab indexing - icefront_fill=(int)*(iomodel->pressureload+segment_width*i+segment_width-1); - icefront_eid=(int) *(iomodel->pressureload+segment_width*i+segment_width-2); //matlab indexing - icefront_element_type=element_type; - - i1=(int)*(iomodel->pressureload+segment_width*i+0); - i2=(int)*(iomodel->pressureload+segment_width*i+1); - - icefront_node_ids[0]=i1; - icefront_node_ids[1]=i2; - - icefront_h[0]=iomodel->thickness[i1-1]; - icefront_h[1]=iomodel->thickness[i2-1]; - - icefront_b[0]=iomodel->bed[i1-1]; - icefront_b[1]=iomodel->bed[i2-1]; - - if ((int)*(iomodel->elements_type+2*element+0)==MacAyealFormulationEnum()){ //this is a collapsed 3d element, icefront will be 2d - strcpy(icefront_type,"segment"); - } - else if ((int)*(iomodel->elements_type+2*element+0)==PattynFormulationEnum()){ //this is a real 3d element, icefront will be 3d. - strcpy(icefront_type,"quad"); - i3=(int)*(iomodel->pressureload+segment_width*i+2); - i4=(int)*(iomodel->pressureload+segment_width*i+3); - icefront_node_ids[2]=i3; - icefront_node_ids[3]=i4; - - icefront_h[2]=iomodel->thickness[i3-1]; - icefront_h[3]=iomodel->thickness[i4-1]; - - icefront_b[2]=iomodel->bed[i3-1]; - icefront_b[3]=iomodel->bed[i4-1]; - } - else{ - ISSMERROR(exprintf(" element type %i not supported yet",(int)*(iomodel->elements_type+2*element+0))); - } - - icefront = new Icefront(icefront_type,icefront_fill,icefront_sid,icefront_mparid,icefront_eid,icefront_element_type,icefront_node_ids,icefront_h,icefront_b); - - loads->AddObject(icefront); + /*Create and add load: */ + loads->AddObject(new Icefront(i,iomodel)); } @@ -160,6 +57,5 @@ xfree((void**)&iomodel->bed); - /*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */ - /*First fetch rifts: */ + /*Create Riffront loads for rifts: */ IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo"); IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); @@ -171,68 +67,12 @@ for(i=0;inumrifts;i++){ - el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2); - #ifdef _PARALLEL_ - if (iomodel->epart[el1-1]!=my_rank){ - /*This load does not belong to this cluster node, as it references an element - *that does not belong to this node's partition. Drop this 'j':*/ - continue; + if(iomodel->my_elements[(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)]){ + + loads->AddObject(new Riftfront(i,iomodel)); } - #endif - - /*Ok, retrieve all the data needed to add a penalty between the two grids: */ - el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3); - - grid1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+0); - grid2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+1); - - normal[0]=*(iomodel->riftinfo+RIFTINFOSIZE*i+4); - normal[1]=*(iomodel->riftinfo+RIFTINFOSIZE*i+5); - length=*(iomodel->riftinfo+RIFTINFOSIZE*i+6); - - fill = (int)*(iomodel->riftinfo+RIFTINFOSIZE*i+7); - friction=*(iomodel->riftinfo+RIFTINFOSIZE*i+8); - fraction=*(iomodel->riftinfo+RIFTINFOSIZE*i+9); - fractionincrement=*(iomodel->riftinfo+RIFTINFOSIZE*i+10); - - strcpy(riftfront_type,"2d"); - riftfront_id=i+1; //matlab indexing - riftfront_node_ids[0]=grid1; - riftfront_node_ids[1]=grid2; - riftfront_mparid=iomodel->numberofelements+1; //matlab indexing - - riftfront_h[0]=iomodel->thickness[grid1-1]; - riftfront_h[1]=iomodel->thickness[grid2-1]; - - riftfront_b[0]=iomodel->bed[grid1-1]; - riftfront_b[1]=iomodel->bed[grid2-1]; - - riftfront_s[0]=iomodel->surface[grid1-1]; - riftfront_s[1]=iomodel->surface[grid2-1]; - - riftfront_normal[0]=normal[0]; - riftfront_normal[1]=normal[1]; - riftfront_length=length; - - riftfront_fill=fill; - riftfront_friction=friction; - riftfront_fraction=fraction; - riftfront_fractionincrement=fractionincrement; - riftfront_shelf=(bool)iomodel->gridoniceshelf[grid1-1]; - - riftfront_penalty_offset=iomodel->penalty_offset; - riftfront_penalty_lock=iomodel->penalty_lock; - - riftfront_active=0; - riftfront_frozen=0; - riftfront_counter=0; - riftfront_prestable=0; - - riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_frozen,riftfront_counter,riftfront_prestable,riftfront_shelf); - - loads->AddObject(riftfront); } } - - + + /*free ressources: */ xfree((void**)&iomodel->riftinfo); Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp (revision 3428) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp (revision 3429) @@ -35,13 +35,13 @@ IoModelFetchData(&vz,NULL,NULL,iomodel_handle,"vz"); - ug=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double)); + ug=(double*)xcalloc(iomodel->numberofvertices*3,sizeof(double)); - if(vx) for(i=0;inumberofnodes;i++)ug[3*i+0]=vx[i]/iomodel->yts; - if(vy) for(i=0;inumberofnodes;i++)ug[3*i+1]=vy[i]/iomodel->yts; - if(vz) for(i=0;inumberofnodes;i++)ug[3*i+2]=vz[i]/iomodel->yts; + if(vx) for(i=0;inumberofvertices;i++)ug[3*i+0]=vx[i]/iomodel->yts; + if(vy) for(i=0;inumberofvertices;i++)ug[3*i+1]=vy[i]/iomodel->yts; + if(vz) for(i=0;inumberofvertices;i++)ug[3*i+2]=vz[i]/iomodel->yts; count++; param= new Param(count,"u_g",DOUBLEVEC); - param->SetDoubleVec(ug,3*iomodel->numberofnodes,3); + param->SetDoubleVec(ug,3*iomodel->numberofvertices,3); parameters->AddObject(param);