Panzer  Version of the Day
Panzer_Response_Probe_impl.hpp
Go to the documentation of this file.
1 #ifndef __Panzer_Response_Probe_impl_hpp__
2 #define __Panzer_Response_Probe_impl_hpp__
3 
4 #include "Teuchos_Comm.hpp"
5 #include "Teuchos_CommHelpers.hpp"
6 #include "Teuchos_dyn_cast.hpp"
7 
8 #include "Epetra_LocalMap.h"
9 
10 #include "Sacado_Traits.hpp"
11 
12 namespace panzer {
13 
14 template <typename EvalT>
16 Response_Probe(const std::string & responseName, MPI_Comm comm,
17  const Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> > & linObjFact)
18  : ResponseMESupport_Default<EvalT>(responseName,comm), value(0.0),
19  have_probe(false), linObjFactory_(linObjFact)
20 {
21  if(linObjFactory_!=Teuchos::null) {
22  // requires thyra object factory
23  thyraObjFactory_ = Teuchos::rcp_dynamic_cast<const panzer::ThyraObjFactory<double> >(linObjFactory_,true);
25 
26  // build a ghosted container, with a solution vector
27  ghostedContainer_ = linObjFactory_->buildGhostedLinearObjContainer();
28 
29  // set ghosted container (work space for assembly)
31  }
32 }
33 
34 template <typename EvalT>
37 {
38  value = 0.0;
39  have_probe = false;
40 
41  if(ghostedContainer_!=Teuchos::null) ghostedContainer_->initialize();
42 }
43 
44 template <typename EvalT>
47 {
48  double glbValue = Sacado::ScalarValue<ScalarT>::eval(value);
49 
50  // find the minimum processor who has the probe value
51  int locProc = have_probe ? this->getComm()->getRank() : this->getComm()->getSize();
52  int glbProc = 0;
53  Teuchos::reduceAll(*this->getComm(), Teuchos::REDUCE_MIN, Thyra::Ordinal(1), &locProc, &glbProc);
54 
55  TEUCHOS_ASSERT(glbProc < this->getComm()->getSize());
56 
57  // now broadcast the value from proc glbProc
58  Teuchos::broadcast(*this->getComm(), glbProc, Thyra::Ordinal(1), &glbValue);
59 
60  value = glbValue;
61 
62  // built data in vectors
63  if(this->useEpetra()) {
64  // use epetra
65  this->getEpetraVector()[0] = glbValue;
66  }
67  else {
68  // use thyra
69  TEUCHOS_ASSERT(this->useThyra());
70 
71  this->getThyraVector()[0] = glbValue;
72  }
73 }
74 
75 template < >
78 {
79  using Teuchos::rcp_dynamic_cast;
80 
81  Teuchos::RCP<Thyra::MultiVectorBase<double> > dgdx_unique = getDerivative();
82 
83  uniqueContainer_ = linObjFactory_->buildLinearObjContainer();
84  Teuchos::rcp_dynamic_cast<ThyraObjContainer<double> >(uniqueContainer_)->set_x_th(dgdx_unique->col(0));
85 
86  linObjFactory_->ghostToGlobalContainer(*ghostedContainer_,*uniqueContainer_,LinearObjContainer::X);
87 
88  uniqueContainer_ = Teuchos::null;
89 }
90 
91 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
92 template < >
95 {
96  using Teuchos::rcp_dynamic_cast;
97 
98  Teuchos::RCP<Thyra::MultiVectorBase<double> > dgdx_unique = getDerivative();
99 
100  uniqueContainer_ = linObjFactory_->buildLinearObjContainer();
101  Teuchos::rcp_dynamic_cast<ThyraObjContainer<double> >(uniqueContainer_)->set_x_th(dgdx_unique->col(0));
102 
103  linObjFactory_->ghostToGlobalContainer(*ghostedContainer_,*uniqueContainer_,LinearObjContainer::X);
104 
105  uniqueContainer_ = Teuchos::null;
106 }
107 #endif
108 
109 template < >
112 {
113  const int n = value.size();
114  const int num_deriv = this->numDeriv();
115  TEUCHOS_ASSERT(n == 0 || n == num_deriv);
116  if (n == 0)
117  value.resize(num_deriv);
118 
119  // find the minimum processor who has the probe value
120  if (num_deriv > 0) {
121  int locProc = have_probe ? this->getComm()->getRank() : this->getComm()->getSize();
122  int glbProc = 0;
123  Teuchos::reduceAll(*this->getComm(), Teuchos::REDUCE_MIN, Thyra::Ordinal(1), &locProc, &glbProc);
124 
125  TEUCHOS_ASSERT(glbProc < this->getComm()->getSize());
126 
127  // now broadcast the derivatives from proc glbProc
128  Teuchos::broadcast(*this->getComm(), glbProc, Thyra::Ordinal(num_deriv), &value.fastAccessDx(0));
129  }
130 
131  // copy data in vectors
132  if(this->useEpetra()) {
133  // use epetra
134  Epetra_MultiVector& deriv = this->getEpetraMultiVector();
135  for (int i=0; i<num_deriv; ++i)
136  deriv[i][0] = value.dx(i);
137  }
138  else {
139  // use thyra
140  TEUCHOS_ASSERT(this->useThyra());
141  Thyra::ArrayRCP< Thyra::ArrayRCP<double> > deriv = this->getThyraMultiVector();
142  for (int i=0; i<num_deriv; ++i)
143  deriv[i][0] = value.dx(i);
144  }
145 }
146 
147 // Do nothing unless derivatives are actually required
148 template <typename EvalT>
151 
152 // derivatives are required for
153 template < >
156 {
157  setDerivativeVectorSpace(soln_vs);
158 }
159 
160 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
161 // derivatives are required for
162 template < >
165 {
166  setDerivativeVectorSpace(soln_vs);
167 }
168 #endif
169 
170 // Do nothing unless derivatives are required
171 template <typename EvalT>
173 adjustForDirichletConditions(const GlobalEvaluationData & localBCRows,const GlobalEvaluationData & globalBCRows) { }
174 
175 // Do nothing unless derivatives are required
176 template < >
179 {
180  linObjFactory_->adjustForDirichletConditions(Teuchos::dyn_cast<const LinearObjContainer>(localBCRows),
181  Teuchos::dyn_cast<const LinearObjContainer>(globalBCRows),
182  *ghostedContainer_,true,true);
183 }
184 
185 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
186 // Do nothing unless derivatives are required
187 template < >
189 adjustForDirichletConditions(const GlobalEvaluationData & localBCRows,const GlobalEvaluationData & globalBCRows)
190 {
191  linObjFactory_->adjustForDirichletConditions(Teuchos::dyn_cast<const LinearObjContainer>(localBCRows),
192  Teuchos::dyn_cast<const LinearObjContainer>(globalBCRows),
193  *ghostedContainer_,true,true);
194 }
195 #endif
196 
197 }
198 
199 #endif
virtual void scatterResponse()
This simply does global summation, then shoves the result into a vector.
Teuchos::RCP< LinearObjContainer > ghostedContainer_
Teuchos::RCP< const panzer::LinearObjFactory< panzer::Traits > > linObjFactory_
void adjustForDirichletConditions(const GlobalEvaluationData &localBCRows, const GlobalEvaluationData &globalBCRows)
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
virtual Teuchos::RCP< const Thyra::VectorSpaceBase< ScalarT > > getThyraDomainSpace() const =0
Get the domain space.
Teuchos::RCP< const panzer::ThyraObjFactory< double > > thyraObjFactory_
Teuchos::RCP< const Teuchos::Comm< int > > comm
void setSolnVectorSpace(const Teuchos::RCP< const Thyra::VectorSpaceBase< double > > &soln_vs)
Set solution vector space.
#define TEUCHOS_ASSERT(assertion_test)