Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23959) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23960) @@ -26,5 +26,23 @@ if(hydrology_model!=HydrologyGlaDSEnum) return; - /*No extra loads*/ + /*Add channels?*/ + bool ischannels; + iomodel->FindConstant(&ischannels,"md.hydrology.ischannels"); + if(ischannels){ + /*Get faces (edges in 2d)*/ + CreateFaces(iomodel); + for(int i=0;inumberoffaces;i++){ + /*Get left and right elements*/ + int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] + + /*Now, if this element is not in the partition*/ + if(!iomodel->my_elements[element]) continue; + + /* Add load */ + loads->AddObject(new Channel(i+1,i,i,iomodel)); + } + + /*Free data: */ + } }/*}}}*/ void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23959) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23960) @@ -28,16 +28,10 @@ Channel::Channel(int channel_id,int i,int index,IoModel* iomodel){/*{{{*/ - /* Intermediary */ - int j; - int pos1,pos2,pos3,pos4; - int num_nodes; - - /*channel constructor data: */ - int channel_elem_ids[2]; - int channel_vertex_ids[2]; - int channel_node_ids[4]; - int channel_type; - - /*Get edge*/ + this->id=channel_id; + this->parameters = NULL; + this->element = NULL; + this->nodes = NULL; + + /*Get edge info*/ int i1 = iomodel->faces[4*index+0]; int i2 = iomodel->faces[4*index+1]; @@ -45,75 +39,20 @@ int e2 = iomodel->faces[4*index+3]; - /*First, see wether this is an internal or boundary edge (if e2=-1)*/ - if(e2==-1){ - /* Boundary edge, only one element */ - num_nodes=2; - channel_type=BoundaryEnum; - channel_elem_ids[0]=e1; - } - else{ - /* internal edge: connected to 2 elements */ - num_nodes=4; - channel_type=InternalEnum; - channel_elem_ids[0]=e1; - channel_elem_ids[1]=e2; - } - - /*1: Get vertices ids*/ + /*Set Element hook (4th column may be -1 for boundary edges)*/ + this->helement = new Hook(&e1,1); + + /*Set Vertices hooks (4th column may be -1 for boundary edges)*/ + int channel_vertex_ids[2]; channel_vertex_ids[0]=i1; channel_vertex_ids[1]=i2; - - /*2: Get node ids*/ - if (channel_type==InternalEnum){ - - /*Now, we must get the nodes of the 4 nodes located on the edge*/ - - /*2: Get the column where these ids are located in the index*/ - pos1=pos2=pos3=pos4=UNDEF; - for(j=0;j<3;j++){ - if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1; - if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1; - if(iomodel->elements[3*(e2-1)+j]==i1) pos3=j+1; - if(iomodel->elements[3*(e2-1)+j]==i2) pos4=j+1; - } - _assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF); - - /*3: We have the id of the elements and the position of the vertices in the index - * we can compute their dofs!*/ - channel_node_ids[0]=3*(e1-1)+pos1; - channel_node_ids[1]=3*(e1-1)+pos2; - channel_node_ids[2]=3*(e2-1)+pos3; - channel_node_ids[3]=3*(e2-1)+pos4; - } - else{ - - /*2: Get the column where these ids are located in the index*/ - pos1=pos2=UNDEF; - for(j=0;j<3;j++){ - if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1; - if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1; - } - _assert_(pos1!=UNDEF && pos2!=UNDEF); - - /*3: We have the id of the elements and the position of the vertices in the index - * we can compute their dofs!*/ - channel_node_ids[0]=3*(e1-1)+pos1; - channel_node_ids[1]=3*(e1-1)+pos2; - } - - /*Ok, we have everything to build the object: */ - this->id=channel_id; - - /*Hooks: */ - this->hnodes =new Hook(channel_node_ids,num_nodes); this->hvertices =new Hook(&channel_vertex_ids[0],2); - this->helement =new Hook(channel_elem_ids,1); // take only the first element for now - - //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. - this->parameters=NULL; - this->element=NULL; - this->nodes=NULL; -} -/*}}}*/ + + /*Set Nodes hooks (!! Assumes P1 CG)*/ + int channel_node_ids[2]; + channel_node_ids[0]=i1; + channel_node_ids[1]=i2; + this->hnodes=new Hook(&channel_node_ids[0],2); + +}/*}}}*/ Channel::~Channel(){/*{{{*/ this->parameters=NULL; @@ -121,6 +60,5 @@ delete hnodes; delete hvertices; -} -/*}}}*/ +}/*}}}*/ /*Object virtual functions definitions:*/ @@ -234,6 +172,11 @@ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - if(analysis_type != HydrologyGlaDSAnalysisEnum) _error_("Analysis not supported!!"); - _error_("STOP"); + switch(analysis_type){ + case HydrologyGlaDSAnalysisEnum: + Ke = this->CreateKMatrixHydrologyGlaDS(); + break; + default: + _error_("Don't know why we should be here"); + } /*Add to global matrix*/ @@ -251,6 +194,12 @@ int analysis_type; this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - if(analysis_type != HydrologyGlaDSAnalysisEnum) _error_("Analysis not supported!!"); - _error_("STOP"); + + switch(analysis_type){ + case HydrologyGlaDSAnalysisEnum: + pe = this->CreatePVectorHydrologyGlaDS(); + break; + default: + _error_("Don't know why we should be here"); + } /*Add to global matrix*/ @@ -372,2 +321,30 @@ } /*}}}*/ + +/*Channel specific functions*/ +ElementVector* Channel::CreatePVectorHydrologyGlaDS(void){/*{{{*/ + + _error_("not implemented :( "); + + /*Initialize Element matrix and return if necessary*/ + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; + ElementVector* pe=new ElementVector(nodes,NUMNODES,this->parameters); + + /*Clean up and return*/ + return pe; +} +/*}}}*/ +ElementMatrix* Channel::CreateKMatrixHydrologyGlaDS(void){/*{{{*/ + + _error_("not implemented :( "); + + /*Initialize Element matrix and return if necessary*/ + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; + ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES,this->parameters); + + /*Clean up and return*/ + return Ke; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.h (revision 23959) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.h (revision 23960) @@ -72,4 +72,6 @@ /*Channel management:{{{*/ void GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]); + ElementVector* CreatePVectorHydrologyGlaDS(void); + ElementMatrix* CreateKMatrixHydrologyGlaDS(void); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23959) +++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23960) @@ -159,18 +159,17 @@ switch(analysis_type){ - - case HydrologyShaktiAnalysisEnum: - pe = this->CreatePVectorHydrologyShakti(); - break; - case HydrologyDCInefficientAnalysisEnum: - pe = CreatePVectorHydrologyDCInefficient(); - break; - case HydrologyDCEfficientAnalysisEnum: - pe = CreatePVectorHydrologyDCEfficient(); - break; - default: - _error_("Don't know why we should be here"); - /*No loads applied, do nothing: */ - return; + case HydrologyShaktiAnalysisEnum: + pe = this->CreatePVectorHydrologyShakti(); + break; + case HydrologyDCInefficientAnalysisEnum: + pe = CreatePVectorHydrologyDCInefficient(); + break; + case HydrologyDCEfficientAnalysisEnum: + pe = CreatePVectorHydrologyDCEfficient(); + break; + default: + _error_("Don't know why we should be here"); + /*No loads applied, do nothing: */ + return; } if(pe){