Index: ../trunk-jpl/src/c/classes/Loads/Channel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23959) +++ ../trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23960) @@ -27,101 +27,39 @@ /*}}}*/ Channel::Channel(int channel_id,int i,int index,IoModel* iomodel){/*{{{*/ - /* Intermediary */ - int j; - int pos1,pos2,pos3,pos4; - int num_nodes; + this->id=channel_id; + this->parameters = NULL; + this->element = NULL; + this->nodes = NULL; - /*channel constructor data: */ - int channel_elem_ids[2]; - int channel_vertex_ids[2]; - int channel_node_ids[4]; - int channel_type; - - /*Get edge*/ + /*Get edge info*/ int i1 = iomodel->faces[4*index+0]; int i2 = iomodel->faces[4*index+1]; int e1 = iomodel->faces[4*index+2]; 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; - } + /*Set Element hook (4th column may be -1 for boundary edges)*/ + this->helement = new Hook(&e1,1); - /*1: Get vertices ids*/ + /*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; + this->hvertices =new Hook(&channel_vertex_ids[0],2); - /*2: Get node ids*/ - if (channel_type==InternalEnum){ + /*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); - /*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; -} -/*}}}*/ +}/*}}}*/ Channel::~Channel(){/*{{{*/ this->parameters=NULL; delete helement; delete hnodes; delete hvertices; -} -/*}}}*/ +}/*}}}*/ /*Object virtual functions definitions:*/ Object* Channel::copy() {/*{{{*/ @@ -233,8 +171,13 @@ int analysis_type; 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*/ if(Ke){ @@ -250,9 +193,15 @@ ElementVector* pe=NULL; 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*/ if(pe){ pe->AddToGlobal(pf); @@ -371,3 +320,31 @@ *po_nz=o_nz; } /*}}}*/ + +/*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: ../trunk-jpl/src/c/classes/Loads/Channel.h =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Channel.h (revision 23959) +++ ../trunk-jpl/src/c/classes/Loads/Channel.h (revision 23960) @@ -71,6 +71,8 @@ /*}}}*/ /*Channel management:{{{*/ void GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]); + ElementVector* CreatePVectorHydrologyGlaDS(void); + ElementMatrix* CreateKMatrixHydrologyGlaDS(void); /*}}}*/ }; Index: ../trunk-jpl/src/c/classes/Loads/Moulin.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23959) +++ ../trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23960) @@ -158,20 +158,19 @@ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); 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){ pe->AddToGlobal(pf); Index: ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23959) +++ ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 23960) @@ -25,7 +25,25 @@ /*Now, do we really want GlaDS?*/ 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){/*{{{*/