Changeset 22739 for issm/trunk-jpl/src/c/classes/Elements/Element.cpp

Timestamp:

05/03/18 13:31:36 (7 years ago)

Author:

erobo

Message:

CHG: allow for transient ControlInputs

File:

• issm/trunk-jpl/src/c/classes/Elements/Element.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -1387 +1387 @@
         /*Get number of controls*/
         int num_controls;
+        bool isautodiff;
         parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+        parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
+
 
         /*Get number of vertices*/
         int numvertices = this->GetNumberOfVertices();
-
+        if(isautodiff){
+                int* N;
+                int* M;
+                int start = 0;
+                parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
+                parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
+                if(control_index>0) for(int n=0;n<control_index-1;n++) start+=N[n]*M[n];
+
+                for(int n=0;n<N[control_index];n++){
+                        for(int i=0;i<numvertices;i++){
+                                indexing[i+n*numvertices]=this->vertices[i]->Sid() + start + n*M[control_index];
+                        }
+                }
+        }
+        else{
+                int M;
+                parameters->FindParam(&M,ControlInputSizeMEnum);
         /*get gradient indices*/
-        if(onsid){
-                for(int i=0;i<numvertices;i++){
-                        indexing[i]=num_controls*this->vertices[i]->Sid() + control_index;
-                }
-        }
-        else{
-                for(int i=0;i<numvertices;i++){
-                        indexing[i]=num_controls*this->vertices[i]->Pid() + control_index;
-                }
-        }
-
+                if(onsid){
+                        for(int i=0;i<numvertices;i++){
+                                indexing[i]=this->vertices[i]->Sid() + (control_index)*M;
+                        }
+                }
+                else{
+                        for(int i=0;i<numvertices;i++){
+                                indexing[i]=this->vertices[i]->Pid() + (control_index)*M;
+                        }
+                }
+        }
 }
 /*}}}*/
@@ -1619 +1638 @@
 
         /*For now we only support nodal vectors*/
-        if(M!=iomodel->numberofvertices) _error_("not supported");
-        if(N!=1) _error_("not supported");
+        //if(M!=iomodel->numberofvertices) _error_("not supported");
+        //if(N!=1) _error_("not supported");
 
         /*Recover vertices ids needed to initialize inputs*/
@@ -1637 +1656 @@
                 this->AddControlInput(input_enum,values,min_vector,max_vector,P1Enum,id);
         }
+
+        else if(M==iomodel->numberofvertices+1){
+            /*create transient input: */
+            IssmDouble* times = xNew<IssmDouble>(N);
+            for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
+                                /*Create the three transient inputs for the control input*/
+            TransientInput* values_input=new TransientInput(ControlInputValuesEnum,times,N);
+                                TransientInput* mins_input = new TransientInput(ControlInputMinsEnum,times,N);
+                                TransientInput* maxs_input = new TransientInput(ControlInputMaxsEnum,times,N);
+                           TransientInput* grad_input = new TransientInput(ControlInputGradEnum);
+                                for(int t=0;t<N;t++){
+                for(int i=0;i<numvertices;i++){ 
+                                                values[i]=vector[N*(vertexids[i]-1)+t];
+                                                values_min[i] = min_vector[N*(vertexids[i]-1)+t];
+                                                values_max[i] = max_vector[N*(vertexids[i]-1)+t];
+                                         } 
+                                        switch(this->ObjectEnum()){
+                    case TriaEnum:
+                                                                        values_input->AddTimeInput(new TriaInput(ControlInputValuesEnum,values,P1Enum)); 
+                                                                        mins_input->AddTimeInput(new TriaInput(ControlInputMinsEnum,values_min,P1Enum)); 
+                                                                        maxs_input->AddTimeInput(new TriaInput(ControlInputMaxsEnum,values_max,P1Enum));
+                                                                break;
+                    case PentaEnum:
+                                                                printf("-------------- file: Element.cpp line: %i\n",__LINE__); 
+                                                                        values_input->AddTimeInput(new PentaInput(ControlInputValuesEnum,values,P1Enum)); 
+                                                                        mins_input->AddTimeInput(new PentaInput(ControlInputMinsEnum,values_min,P1Enum)); 
+                                                                        maxs_input->AddTimeInput(new PentaInput(ControlInputMaxsEnum,values_max,P1Enum)); 
+                                                                        break;
+                    case TetraEnum:
+                                                                        values_input->AddTimeInput(new TetraInput(ControlInputValuesEnum,values,P1Enum)); 
+                                                                        mins_input->AddTimeInput(new TetraInput(ControlInputMinsEnum,values_min,P1Enum)); 
+                                                                        maxs_input->AddTimeInput(new TetraInput(ControlInputMaxsEnum,values_max,P1Enum)); 
+                                                                        break;
+                    default: _error_("Not implemented yet");
+                }
+            }
+            this->inputs->AddInput(new ControlInput(input_enum,TransientInputEnum,values_input,mins_input,maxs_input,grad_input,P1Enum,id));
+            xDelete<IssmDouble>(times);
+        }
+                  else _error_("not currently supported type of M and N attempted");
 
         /*clean up*/