30 #ifndef FE_JAC_FILL_FUNCS_HPP 31 #define FE_JAC_FILL_FUNCS_HPP 52 #ifdef PACKAGE_BUGREPORT 53 #undef PACKAGE_BUGREPORT 58 #ifdef PACKAGE_TARNAME 59 #undef PACKAGE_TARNAME 61 #ifdef PACKAGE_VERSION 62 #undef PACKAGE_VERSION 68 #define NUMBER_DIRECTIONS 100 69 #include "adolc/adouble.h" 70 #include "adolc/drivers/drivers.h" 71 #include "adolc/interfaces.h" 90 #define ad_GRAD_MAX 130 99 static const unsigned int nqp = 2;
100 static const unsigned int nnode = 2;
110 for (
unsigned int i=0; i<
nqp; i++) {
115 jac[i] = 0.5*mesh_spacing;
118 phi[i][0] = 0.5*(1.0 - xi[i]);
119 phi[i][1] = 0.5*(1.0 + xi[i]);
126 template <
class FadType>
129 const std::vector<double>& x,
130 std::vector<FadType>& x_fad) {
131 for (
unsigned int node=0; node<e.
nnode; node++)
132 for (
unsigned int eqn=0; eqn<neqn; eqn++)
133 x_fad[node*neqn+eqn] =
FadType(e.
nnode*neqn, node*neqn+eqn,
134 x[e.
gid[node]*neqn+eqn]);
140 const std::vector<T>& x,
145 for (
unsigned int qp=0; qp<e.
nqp; qp++) {
146 for (
unsigned int eqn=0; eqn<neqn; eqn++) {
147 u[qp*neqn+eqn] = 0.0;
148 du[qp*neqn+eqn] = 0.0;
149 for (
unsigned int node=0; node<e.
nnode; node++) {
150 u[qp*neqn+eqn] += x[node*neqn+eqn] * e.
phi[qp][node];
151 du[qp*neqn+eqn] += x[node*neqn+eqn] * e.
dphi[qp][node];
157 std::vector<T> s(e.
nqp);
158 for (
unsigned int qp=0; qp<e.
nqp; qp++) {
160 for (
unsigned int eqn=0; eqn<neqn; eqn++)
161 s[qp] += u[qp*neqn+eqn]*u[qp*neqn+eqn];
165 for (
unsigned int node=0; node<e.
nnode; node++) {
166 for (
unsigned int eqn=0; eqn<neqn; eqn++) {
167 unsigned int row = node*neqn+eqn;
169 for (
unsigned int qp=0; qp<e.
nqp; qp++) {
170 double c1 = e.
w[qp]*e.
jac[qp];
171 double c2 = -e.
dphi[qp][node]/(e.
jac[qp]*e.
jac[qp]);
173 c1*(c2*du[qp*neqn+eqn] + e.
phi[qp][node]*s[qp]*
exp(u[qp*neqn+eqn]));
179 template <
class FadType>
182 const std::vector<FadType>& f_fad,
183 std::vector<double>& f,
184 std::vector< std::vector<double> >& jac) {
185 for (
unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
186 f[e.
gid[0]*neqn+eqn_row] += f_fad[0*neqn+eqn_row].val();
187 f[e.
gid[1]*neqn+eqn_row] += f_fad[1*neqn+eqn_row].val();
188 for (
unsigned int node_col=0; node_col<e.
nnode; node_col++) {
189 for (
unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
190 unsigned int col = node_col*neqn+eqn_col;
191 unsigned int next_col = (node_col+1)*neqn+eqn_col;
192 jac[e.
gid[0]*neqn+eqn_row][next_col] +=
193 f_fad[0*neqn+eqn_row].fastAccessDx(col);
194 jac[e.
gid[1]*neqn+eqn_row][col] +=
195 f_fad[1*neqn+eqn_row].fastAccessDx(col);
201 template <
class FadType>
203 double mesh_spacing) {
207 std::vector<double> x(num_nodes*num_eqns),
f(num_nodes*num_eqns);
208 std::vector< std::vector<double> > jac(num_nodes*num_eqns);
209 for (
unsigned int node_row=0; node_row<num_nodes; node_row++) {
210 for (
unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) {
211 unsigned int row = node_row*num_eqns + eqn_row;
212 x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
214 jac[row] = std::vector<double>((e.
nnode+1)*num_eqns);
215 for (
unsigned int node_col=0; node_col<e.
nnode+1; node_col++) {
216 for (
unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) {
217 unsigned int col = node_col*num_eqns + eqn_col;
226 std::vector<FadType> x_fad(e.
nnode*num_eqns), f_fad(e.
nnode*num_eqns);
227 std::vector<FadType> u(e.
nqp*num_eqns), du(e.
nqp*num_eqns);
228 for (
unsigned int i=0; i<num_nodes-1; i++) {
255 double mesh_spacing);
256 double residual_fill(
unsigned int num_nodes,
unsigned int num_eqns,
257 double mesh_spacing);
260 #ifndef ADOLC_TAPELESS 261 void adolc_init_fill(
bool retape,
265 const std::vector<double>& x,
266 std::vector<double>& x_local,
267 std::vector<adouble>& x_ad);
269 void adolc_process_fill(
bool retape,
273 std::vector<double>& x_local,
274 std::vector<adouble>& f_ad,
275 std::vector<double>& f,
276 std::vector<double>& f_local,
277 std::vector< std::vector<double> >& jac,
281 double adolc_jac_fill(
unsigned int num_nodes,
unsigned int num_eqns,
282 double mesh_spacing);
284 double adolc_retape_jac_fill(
unsigned int num_nodes,
unsigned int num_eqns,
285 double mesh_spacing);
289 void adolc_tapeless_init_fill(
const ElemData& e,
291 const std::vector<double>& x,
292 std::vector<adtl::adouble>& x_ad);
294 void adolc_tapeless_process_fill(
const ElemData& e,
296 const std::vector<adtl::adouble>& f_ad,
297 std::vector<double>& f,
298 std::vector< std::vector<double> >& jac);
300 double adolc_tapeless_jac_fill(
unsigned int num_nodes,
unsigned int num_eqns,
301 double mesh_spacing);
307 void adic_init_fill(
const ElemData& e,
309 const std::vector<double>& x,
310 std::vector<DERIV_TYPE>& x_fad);
312 void adic_process_fill(
const ElemData& e,
314 const std::vector<DERIV_TYPE>& f_fad,
315 std::vector<double>& f,
316 std::vector< std::vector<double> >& jac);
317 double adic_jac_fill(
unsigned int num_nodes,
unsigned int num_eqns,
318 double mesh_spacing);
319 inline void AD_Init(
int arg0) {
320 ad_AD_GradInit(arg0);
double residual_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
void template_element_fill(const ElemData &e, unsigned int neqn, const std::vector< T > &x, std::vector< T > &u, std::vector< T > &du, std::vector< T > &f)
void adic_element_fill(ElemData *e, unsigned int neqn, const DERIV_TYPE *x, DERIV_TYPE *u, DERIV_TYPE *du, DERIV_TYPE *f)
Sacado::Fad::DFad< double > FadType
SimpleFad< ValueT > sqrt(const SimpleFad< ValueT > &a)
ElemData(double mesh_spacing)
void fad_init_fill(const ElemData &e, unsigned int neqn, const std::vector< double > &x, std::vector< FadType > &x_fad)
SimpleFad< ValueT > exp(const SimpleFad< ValueT > &a)
double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
void start(bool reset=false)
double totalElapsedTime(bool readCurrentTime=false) const
void fad_process_fill(const ElemData &e, unsigned int neqn, const std::vector< FadType > &f_fad, std::vector< double > &f, std::vector< std::vector< double > > &jac)
double fad_jac_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)