53 #ifndef ROL_POISSONCONTROL_HPP 54 #define ROL_POISSONCONTROL_HPP 71 typedef typename vector::size_type
uint;
76 Teuchos::RCP<const vector>
getVector(
const V& x ) {
77 using Teuchos::dyn_cast;
82 using Teuchos::dyn_cast;
97 Real h = 1.0/((Real)n+1.0);
98 for (uint i=0; i<n; i++) {
100 (*Mzp)[i] = h/6.0*(4.0*(*zp)[i] + (*zp)[i+1]);
102 else if ( i == n-1 ) {
103 (*Mzp)[i] = h/6.0*((*zp)[i-1] + 4.0*(*zp)[i]);
106 (*Mzp)[i] = h/6.0*((*zp)[i-1] + 4.0*(*zp)[i] + (*zp)[i+1]);
119 Real h = 1.0/((Real)n+1.0);
120 SV b( Teuchos::rcp(
new vector(n,0.0) ) );
129 (*up)[0] = (*bp)[0]/d;
130 for ( uint i = 1; i < n; i++ ) {
131 m = 1.0/(d - o*c[i-1]);
133 (*up)[i] = ( (*bp)[i] - o*(*up)[i-1] )*m;
135 for ( uint i = n-1; i > 0; i-- ) {
136 (*up)[i-1] = (*up)[i-1] - c[i-1]*(*up)[i];
141 Real val = 1.0/3.0*std::pow(x,4.0) - 2.0/3.0*std::pow(x,3.0) + 1.0/3.0*x + 8.0*
alpha_;
151 Real h = 1.0/((Real)n+1.0);
153 SV u( rcp(
new vector(n,0.0) ) );
162 for (uint i=0; i<n; i++) {
163 res = alpha_*(*zp)[i];
165 res *= h/6.0*(4.0*(*zp)[i] + (*zp)[i+1]);
168 res += h/6.0*(4.0*res1 + res2)*res1;
170 else if ( i == n-1 ) {
171 res *= h/6.0*((*zp)[i-1] + 4.0*(*zp)[i]);
174 res += h/6.0*(res1 + 4.0*res2)*res2;
177 res *= h/6.0*((*zp)[i-1] + 4.0*(*zp)[i] + (*zp)[i+1]);
181 res += h/6.0*(res1 + 4.0*res2 + res3)*res2;
196 Real h = 1.0/((Real)n+1.0);
199 SV u( rcp(
new vector(n,0.0) ) );
207 for (uint i=0; i<n; i++) {
211 SV p( rcp(
new vector(n,0.0) ) );
218 for (uint i=0; i<n; i++) {
220 res1 = alpha_*(*zp)[i] - (*pp)[i];
221 res2 = alpha_*(*zp)[i+1] - (*pp)[i+1];
222 (*gp)[i] = h/6.0*(4.0*res1 + res2);
224 else if ( i == n-1 ) {
225 res1 = alpha_*(*zp)[i-1] - (*pp)[i-1];
226 res2 = alpha_*(*zp)[i] - (*pp)[i];
227 (*gp)[i] = h/6.0*(res1 + 4.0*res2);
230 res1 = alpha_*(*zp)[i-1] - (*pp)[i-1];
231 res2 = alpha_*(*zp)[i] - (*pp)[i];
232 res3 = alpha_*(*zp)[i+1] - (*pp)[i+1];
233 (*gp)[i] = h/6.0*(res1 + 4.0*res2 + res3);
247 Real h = 1.0/((Real)n+1.0);
250 SV u( rcp(
new vector(n,0.0) ) );
255 SV p( rcp(
new vector(n,0.0) ) );
263 for (uint i=0; i<n; i++) {
265 res1 = alpha_*(*vp)[i] + (*pp)[i];
266 res2 = alpha_*(*vp)[i+1] + (*pp)[i+1];
267 (*hvp)[i] = h/6.0*(4.0*res1 + res2);
269 else if ( i == n-1 ) {
270 res1 = alpha_*(*vp)[i-1] + (*pp)[i-1];
271 res2 = alpha_*(*vp)[i] + (*pp)[i];
272 (*hvp)[i] = h/6.0*(res1 + 4.0*res2);
275 res1 = alpha_*(*vp)[i-1] + (*pp)[i-1];
276 res2 = alpha_*(*vp)[i] + (*pp)[i];
277 res3 = alpha_*(*vp)[i+1] + (*pp)[i+1];
278 (*hvp)[i] = h/6.0*(res1 + 4.0*res2 + res3);
293 Teuchos::RCP<std::vector<Real> > x0p = Teuchos::rcp(
new std::vector<Real>(n,0.0));
294 for (
int i=0; i<n; i++) {
300 Teuchos::RCP<std::vector<Real> > xp = Teuchos::rcp(
new std::vector<Real>(n,0.0));
301 Real h = 1.0/((Real)n+1.0), pt = 0.0;
302 for(
int i = 0; i < n; i++ ) {
304 (*xp)[i] = 4.0*pt*(1.0-pt);
Provides the interface to evaluate objective functions.
Objective_PoissonControl(Real alpha=1.e-4)
Real value(const Vector< Real > &z, Real &tol)
Compute value.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Teuchos::RCP< vector > getVector(V &x)
void getPoissonControl(Teuchos::RCP< Objective< Real > > &obj, Teuchos::RCP< Vector< Real > > &x0, Teuchos::RCP< Vector< Real > > &x)
Defines the linear algebra or vector space interface.
std::vector< Real > vector
void solve_poisson(Vector< Real > &u, const Vector< Real > &z)
void gradient(Vector< Real > &g, const Vector< Real > &z, Real &tol)
Compute gradient.
Poisson distributed control.
void apply_mass(Vector< Real > &Mz, const Vector< Real > &z)
Teuchos::RCP< const vector > getVector(const V &x)
Real evaluate_target(Real x)