oomph::ScaleCoupledElement< ELEMENT > Class Template Reference

Wrapper class for solid elements to be coupled with discrete solid particles on the surfaces. More...

#include <ScaleCoupledElement.h>

+ Inheritance diagram for oomph::ScaleCoupledElement< ELEMENT >:

Public Member Functions

void clear_coupling_residual ()
 Empties the coupling residual. More...
 
void set_coupling_residual (Vector< Vector< double >> &residual)
 Sets the coupling residual. More...
 
double get_coupling_residual (const unsigned &n, const unsigned &d)
 Returns the coupling residual. More...
 
void clear_coupling_jacobian ()
 Empties the coupling Jacobian. More...
 
void set_coupling_jacobian (Vector< Vector< double >> &jacobian)
 Sets the coupling residual. More...
 
void set_coupling_weight (Vector< double > &weight)
 Sets the coupling weight. More...
 
double get_coupling_weight (const unsigned &n)
 Sets the coupling weight. More...
 
void fill_in_contribution_to_residuals (Vector< double > &residuals) override
 Adds the elements contribution to its residual vector (wrapper) More...
 
void fill_in_contribution_to_jacobian (Vector< double > &residuals, DenseMatrix< double > &jacobian) override
 Add the elements contribution to its residual vector and Jacobian matrix (wrapper) More...
 
void fill_in_generic_contribution_to_residuals_pvd (Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag) override
 \TW I think this just addes the coupling weight More...
 
void add_internal_coupling_forces_to_residuals (Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
 Add the point source contribution to the residual vector. More...
 

Private Attributes

Vector< Vector< double > > coupling_residual
 Contribution of coupling to residual (residuals[local_eqn] += coupling_residual[node][dim];) More...
 
Vector< Vector< double > > coupling_jacobian
 Contribution of coupling to jacobian (jacobian(local_eqn,local_unknown) += coupling_jacobian[node0][node1];) More...
 
Vector< doublecoupling_weight
 Nodal coupling weight (interpolated to compute coupling weight at integration points) More...
 

Detailed Description

template<class ELEMENT>
class oomph::ScaleCoupledElement< ELEMENT >

Wrapper class for solid elements to be coupled with discrete solid particles on the surfaces.

Member Function Documentation

◆ add_internal_coupling_forces_to_residuals()

template<class ELEMENT >
void oomph::ScaleCoupledElement< ELEMENT >::add_internal_coupling_forces_to_residuals ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
const unsigned &  flag 
)
inline

Add the point source contribution to the residual vector.

484  {
485  // No further action if no coupling forces
486  if (coupling_residual.empty()) return;
487 
488  // Find out how many nodes there are
489  const unsigned n_node = this->nnode();
490  const unsigned DIM = this->dim();
491 
492  // Find out how many positional dofs there are
493  unsigned n_position_type = this->nnodal_position_type();
494 
495  // Set up memory for the shape/test functions
496  Shape psi(n_node);
497 
498  // Loop over all nodes belonging to the element
499  for(unsigned l=0;l<n_node;l++)
500  {
501  // Loop of types of dofs
502  for(unsigned k=0; k<n_position_type; k++)
503  {
504  // Loop over the force components
505  for(unsigned i=0; i<DIM; i++)
506  {
507  // Get the local equation
508  int local_eqn = this->position_local_eqn(l,k,i);
509  // IF it's not a boundary condition and the interpolated force is nonzero
510  if (local_eqn >= 0)
511  {
512  // add the nodal coupling residual to the global residual vector
513  residuals[local_eqn] += coupling_residual[l][i];
514  }
515  // Get Jacobian too?
516  if (flag==1)
517  {
518  // No further action if no coupling forces
519  if (coupling_jacobian.empty()) continue;
520 
521  // Loop over the nodes of the element again
522  for(unsigned ll=0;ll<n_node;ll++)
523  {
524  // Loop of types of dofs again
525  for(unsigned kk=0;kk<n_position_type;kk++)
526  {
527  // Loop over the force components again
528  for(unsigned ii=0;ii<DIM;ii++)
529  {
530  // Get the number of the unknown
531  int local_unknown = this->position_local_eqn(ll,kk,ii);
532  // IF it's not a boundary condition
533  if ((local_unknown >= 0) && (i==ii) && (local_unknown >= local_eqn))
534  {
535  // add the nodal coupling jacobian to the global jacobian matrix
536  jacobian(local_eqn,local_unknown) += coupling_jacobian[l][ll];
537  if(local_eqn != local_unknown)
538  {
539  jacobian(local_unknown,local_eqn) += coupling_jacobian[l][ll];
540  }
541  }
542  }
543  }
544  }
545  }
546  }
547  }
548  }
549  }
Vector< Vector< double > > coupling_jacobian
Contribution of coupling to jacobian (jacobian(local_eqn,local_unknown) += coupling_jacobian[node0][n...
Definition: ScaleCoupledElement.h:560
Vector< Vector< double > > coupling_residual
Contribution of coupling to residual (residuals[local_eqn] += coupling_residual[node][dim];)
Definition: ScaleCoupledElement.h:555
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51

References oomph::ScaleCoupledElement< ELEMENT >::coupling_jacobian, oomph::ScaleCoupledElement< ELEMENT >::coupling_residual, and constants::i.

Referenced by oomph::ScaleCoupledElement< ELEMENT >::fill_in_contribution_to_jacobian(), and oomph::ScaleCoupledElement< ELEMENT >::fill_in_contribution_to_residuals().

◆ clear_coupling_jacobian()

template<class ELEMENT >
void oomph::ScaleCoupledElement< ELEMENT >::clear_coupling_jacobian ( )
inline

Empties the coupling Jacobian.

57  {
58  coupling_jacobian.clear();
59  }

References oomph::ScaleCoupledElement< ELEMENT >::coupling_jacobian.

◆ clear_coupling_residual()

template<class ELEMENT >
void oomph::ScaleCoupledElement< ELEMENT >::clear_coupling_residual ( )
inline

Empties the coupling residual.

41  {
42  coupling_residual.clear();
43  }

References oomph::ScaleCoupledElement< ELEMENT >::coupling_residual.

◆ fill_in_contribution_to_jacobian()

template<class ELEMENT >
void oomph::ScaleCoupledElement< ELEMENT >::fill_in_contribution_to_jacobian ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian 
)
inlineoverride

Add the elements contribution to its residual vector and Jacobian matrix (wrapper)

90  {
91  //Call the generic routine with the flag set to 1
93 
94  // Add nodal coupling force to the residuals
95  add_internal_coupling_forces_to_residuals(residuals,jacobian,1);
96 
97  //if (!coupling_residual.empty()) logger(INFO,"fill_in_contribution_to_jacobian, element %", this);
98  }
void add_internal_coupling_forces_to_residuals(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add the point source contribution to the residual vector.
Definition: ScaleCoupledElement.h:481
void fill_in_generic_contribution_to_residuals_pvd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag) override
\TW I think this just addes the coupling weight
Definition: ScaleCoupledElement.h:101

References oomph::ScaleCoupledElement< ELEMENT >::add_internal_coupling_forces_to_residuals(), and oomph::ScaleCoupledElement< ELEMENT >::fill_in_generic_contribution_to_residuals_pvd().

◆ fill_in_contribution_to_residuals()

template<class ELEMENT >
void oomph::ScaleCoupledElement< ELEMENT >::fill_in_contribution_to_residuals ( Vector< double > &  residuals)
inlineoverride

Adds the elements contribution to its residual vector (wrapper)

78  {
79  // Call the generic residuals function with flag set to 0 using a dummy matrix argument
80  fill_in_generic_contribution_to_residuals_pvd(residuals,GeneralisedElement::Dummy_matrix,0);
81 
82  // Add nodal coupling force to the residuals
83  add_internal_coupling_forces_to_residuals(residuals,GeneralisedElement::Dummy_matrix,0);
84 
85  //if (!coupling_residual.empty()) logger(INFO,"fill_in_contribution_to_residuals, element %", this);
86  }

References oomph::ScaleCoupledElement< ELEMENT >::add_internal_coupling_forces_to_residuals(), and oomph::ScaleCoupledElement< ELEMENT >::fill_in_generic_contribution_to_residuals_pvd().

◆ fill_in_generic_contribution_to_residuals_pvd()

template<class ELEMENT >
void oomph::ScaleCoupledElement< ELEMENT >::fill_in_generic_contribution_to_residuals_pvd ( Vector< double > &  residuals,
DenseMatrix< double > &  jacobian,
const unsigned &  flag 
)
inlineoverride

\TW I think this just addes the coupling weight

104  {
105 
106  // Get problem dimension
107  const unsigned DIM = this->dim();
108 
109  // Simply set up initial condition?
110  if (this->Solid_ic_pt!=0)
111  {
112  this->fill_in_residuals_for_solid_ic(residuals);
113  return;
114  }
115 
116  //Find out how many nodes there are
117  const unsigned n_node = this->nnode();
118 
119  //Find out how many positional dofs there are
120  const unsigned n_position_type = this->nnodal_position_type();
121 
122  //Set up memory for the shape functions
123  Shape psi(n_node,n_position_type);
124  DShape dpsidxi(n_node,n_position_type,DIM);
125 
126  //Set the value of Nintpt -- the number of integration points
127  const unsigned n_intpt = this->integral_pt()->nweight();
128 
129  //Set the vector to hold the local coordinates in the element
130  Vector<double> s(DIM);
131 
132  // Timescale ratio (non-dim density)
133  double lambda_sq = this->lambda_sq();
134 
135  // Time factor
136  double time_factor=0.0;
137  if (lambda_sq>0)
138  {
139  time_factor=this->node_pt(0)->position_time_stepper_pt()->weight(2,0);
140  }
141 
142  //Integer to store the local equation number
143  int local_eqn=0;
144 
145  //Loop over the integration points
146  for(unsigned ipt=0;ipt<n_intpt;ipt++)
147  {
148  //Assign the values of s
149  for(unsigned i=0;i<DIM;++i) {s[i] = this->integral_pt()->knot(ipt,i);}
150 
151  //Get the integral weight
152  double w = this->integral_pt()->weight(ipt);
153 
154  //Call the derivatives of the shape functions (and get Jacobian)
155  double J = this->dshape_lagrangian_at_knot(ipt,psi,dpsidxi);
156 
157  //Get the coupling weight at the integration point
158  double coupling_w = 0;
159  if (!coupling_weight.empty()) {
160  for (unsigned l = 0; l < n_node; l++) {
161  double psi_ = psi(l);
162  coupling_w += psi_ * coupling_weight[l];
163  }
164  }
165 
166  //Calculate interpolated values of the derivative of global position
167  //wrt lagrangian coordinates
168  DenseMatrix<double> interpolated_G(DIM);
169 
170  // Setup memory for accelerations
171  Vector<double> accel(DIM);
172 
173  //Initialise to zero
174  for(unsigned i=0;i<DIM;i++)
175  {
176  // Initialise acclerations
177  accel[i]=0.0;
178  for(unsigned j=0;j<DIM;j++)
179  {
180  interpolated_G(i,j) = 0.0;
181  }
182  }
183 
184  //Storage for Lagrangian coordinates (initialised to zero)
185  Vector<double> interpolated_xi(DIM,0.0);
186 
187  //Calculate displacements and derivatives and lagrangian coordinates
188  for(unsigned l=0;l<n_node;l++)
189  {
190  //Loop over positional dofs
191  for(unsigned k=0;k<n_position_type;k++)
192  {
193  double psi_ = psi(l,k);
194  //Loop over displacement components (deformed position)
195  for(unsigned i=0;i<DIM;i++)
196  {
197  //Calculate the Lagrangian coordinates and the accelerations
198  interpolated_xi[i] += this->lagrangian_position_gen(l,k,i)*psi_;
199 
200  // Only compute accelerations if inertia is switched on
201  if ((lambda_sq>0.0)&&(this->Unsteady))
202  {
203  accel[i] += this->dnodal_position_gen_dt(2,l,k,i)*psi_;
204  }
205 
206  //Loop over derivative directions
207  for(unsigned j=0;j<DIM;j++)
208  {
209  interpolated_G(j,i) +=
210  this->nodal_position_gen(l,k,i)*dpsidxi(l,k,j);
211  }
212  }
213  }
214  }
215 
216  //Get isotropic growth factor
217  double gamma=1.0;
218  this->get_isotropic_growth(ipt,s,interpolated_xi,gamma);
219 
220 
221  //Get body force at current time
222  Vector<double> b(DIM);
223  this->body_force(interpolated_xi,b);
224 
225  // We use Cartesian coordinates as the reference coordinate
226  // system. In this case the undeformed metric tensor is always
227  // the identity matrix -- stretched by the isotropic growth
228  double diag_entry=pow(gamma,2.0/double(DIM));
229  DenseMatrix<double> g(DIM);
230  for(unsigned i=0;i<DIM;i++)
231  {
232  for(unsigned j=0;j<DIM;j++)
233  {
234  if(i==j) {g(i,j) = diag_entry;}
235  else {g(i,j) = 0.0;}
236  }
237  }
238 
239  //Premultiply the undeformed volume ratio (from the isotropic
240  // growth), the integral weights, the coupling weights, and the Jacobian
241  double W = gamma*w*(1.0-coupling_w)*J;
242 
243  //Declare and calculate the deformed metric tensor
244  DenseMatrix<double> G(DIM);
245 
246  //Assign values of G
247  for(unsigned i=0;i<DIM;i++)
248  {
249  //Do upper half of matrix
250  for(unsigned j=i;j<DIM;j++)
251  {
252  //Initialise G(i,j) to zero
253  G(i,j) = 0.0;
254  //Now calculate the dot product
255  for(unsigned k=0;k<DIM;k++)
256  {
257  G(i,j) += interpolated_G(i,k)*interpolated_G(j,k);
258  }
259  }
260  //Matrix is symmetric so just copy lower half
261  for(unsigned j=0;j<i;j++)
262  {
263  G(i,j) = G(j,i);
264  }
265  }
266 
267  //Now calculate the stress tensor from the constitutive law
268  DenseMatrix<double> sigma(DIM);
269  ELEMENT::get_stress(g,G,sigma);
270 
271  // Add pre-stress
272  for(unsigned i=0;i<DIM;i++)
273  {
274  for(unsigned j=0;j<DIM;j++)
275  {
276  sigma(i,j) += this->prestress(i,j,interpolated_xi);
277  }
278  }
279 
280  // Get stress derivative by FD only needed for Jacobian
281  //-----------------------------------------------------
282 
283  // Stress derivative
284  RankFourTensor<double> d_stress_dG(DIM,DIM,DIM,DIM,0.0);
285  // Derivative of metric tensor w.r.t. to nodal coords
286  RankFiveTensor<double> d_G_dX(n_node,n_position_type,DIM,DIM,DIM,0.0);
287 
288  // Get Jacobian too?
289  if (flag==1)
290  {
291  // Derivative of metric tensor w.r.t. to discrete positional dofs
292  // NOTE: Since G is symmetric we only compute the upper triangle
293  // and DO NOT copy the entries across. Subsequent computations
294  // must (and, in fact, do) therefore only operate with upper
295  // triangular entries
296  for (unsigned ll=0;ll<n_node;ll++)
297  {
298  for (unsigned kk=0;kk<n_position_type;kk++)
299  {
300  for (unsigned ii=0;ii<DIM;ii++)
301  {
302  for (unsigned aa=0;aa<DIM;aa++)
303  {
304  for (unsigned bb=aa;bb<DIM;bb++)
305  {
306  d_G_dX(ll,kk,ii,aa,bb)=
307  interpolated_G(aa,ii)*dpsidxi(ll,kk,bb)+
308  interpolated_G(bb,ii)*dpsidxi(ll,kk,aa);
309  }
310  }
311  }
312  }
313  }
314 
315  //Get the "upper triangular" entries of the derivatives of the stress
316  //tensor with respect to G
317  this->get_d_stress_dG_upper(g,G,sigma,d_stress_dG);
318  }
319 
320 //=====EQUATIONS OF ELASTICITY FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
321 
322  //Loop over the test functions, nodes of the element
323  for(unsigned l=0;l<n_node;l++)
324  {
325  //Loop of types of dofs
326  for(unsigned k=0;k<n_position_type;k++)
327  {
328  // Offset for faster access
329  const unsigned offset5=dpsidxi.offset(l ,k);
330 
331  //Loop over the displacement components
332  for(unsigned i=0;i<DIM;i++)
333  {
334  //Get the equation number
335  local_eqn = this->position_local_eqn(l,k,i);
336 
337  /*IF it's not a boundary condition*/
338  if(local_eqn >= 0)
339  {
340  //Initialise contribution to sum
341  double sum=0.0;
342 
343  // Acceleration and body force
344  sum+=(lambda_sq*accel[i]-b[i])*psi(l,k);
345 
346  // Stress term
347  for(unsigned a=0;a<DIM;a++)
348  {
349  unsigned count=offset5;
350  for(unsigned b=0;b<DIM;b++)
351  {
352  //Add the stress terms to the residuals
353  sum+=sigma(a,b)*interpolated_G(a,i)*
354  dpsidxi.raw_direct_access(count);
355  ++count;
356  }
357  }
358  residuals[local_eqn] += W*sum;
359 
360  // Get Jacobian too?
361  if (flag==1)
362  {
363  // Offset for faster access in general stress loop
364  const unsigned offset1=d_G_dX.offset( l, k, i);
365 
366  //Loop over the nodes of the element again
367  for(unsigned ll=0;ll<n_node;ll++)
368  {
369  //Loop of types of dofs again
370  for(unsigned kk=0;kk<n_position_type;kk++)
371  {
372  //Loop over the displacement components again
373  for(unsigned ii=0;ii<DIM;ii++)
374  {
375  //Get the number of the unknown
376  int local_unknown = this->position_local_eqn(ll,kk,ii);
377 
378  /*IF it's not a boundary condition*/
379  if(local_unknown >= 0)
380  {
381  // Offset for faster access in general stress loop
382  const unsigned offset2=d_G_dX.offset( ll, kk, ii);
383  const unsigned offset4=dpsidxi.offset(ll, kk);
384 
385  // General stress term
386  //--------------------
387  double sum=0.0;
388  unsigned count1=offset1;
389  for(unsigned a=0;a<DIM;a++)
390  {
391  // Bump up direct access because we're only
392  // accessing upper triangle
393  count1+=a;
394  for(unsigned b=a;b<DIM;b++)
395  {
396  double factor=d_G_dX.raw_direct_access(count1);
397  if (a==b) factor*=0.5;
398 
399  // Offset for faster access
400  unsigned offset3=d_stress_dG.offset(a,b);
401  unsigned count2=offset2;
402  unsigned count3=offset3;
403 
404  for(unsigned aa=0;aa<DIM;aa++)
405  {
406  // Bump up direct access because we're only
407  // accessing upper triangle
408  count2+=aa;
409  count3+=aa;
410 
411  // Only upper half of derivatives w.r.t. symm tensor
412  for(unsigned bb=aa;bb<DIM;bb++)
413  {
414  sum+=factor*
415  d_stress_dG.raw_direct_access(count3)*
416  d_G_dX.raw_direct_access(count2);
417  ++count2;
418  ++count3;
419  }
420  }
421  ++count1;
422  }
423 
424  }
425 
426  // Multiply by weight and add contribution
427  // (Add directly because this bit is nonsymmetric)
428  jacobian(local_eqn,local_unknown)+=sum*W;
429 
430  // Only upper triangle (no separate test for bc as
431  // local_eqn is already nonnegative)
432  if((i==ii) && (local_unknown >= local_eqn))
433  {
434  //Initialise contribution
435  double sum=0.0;
436 
437  // Inertia term
438  sum+=lambda_sq*time_factor*psi(ll,kk)*psi(l,k);
439 
440  // Stress term
441  unsigned count4=offset4;
442  for(unsigned a=0;a<DIM;a++)
443  {
444  //Cache term
445  const double factor=
446  dpsidxi.raw_direct_access(count4);// ll ,kk
447  ++count4;
448 
449  unsigned count5=offset5;
450  for(unsigned b=0;b<DIM;b++)
451  {
452  sum+=sigma(a,b)*factor*
453  dpsidxi.raw_direct_access(count5); // l ,k
454  ++count5;
455  }
456  }
457  //Add contribution to jacobian
458  jacobian(local_eqn,local_unknown) += sum*W;
459  //Add to lower triangular section
460  if(local_eqn != local_unknown)
461  {
462  jacobian(local_unknown,local_eqn) += sum*W;
463  }
464  }
465 
466  } //End of if not boundary condition
467  }
468  }
469  }
470  }
471 
472  } //End of if not boundary condition
473 
474  } //End of loop over coordinate directions
475  } //End of loop over type of dof
476  } //End of loop over shape functions
477  } //End of loop over integration points
478  }
Vector< double > coupling_weight
Nodal coupling weight (interpolated to compute coupling weight at integration points)
Definition: ScaleCoupledElement.h:563
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
Definition: ExtendedMath.cc:137

References oomph::ScaleCoupledElement< ELEMENT >::coupling_weight, mathsFunc::gamma(), and constants::i.

Referenced by oomph::ScaleCoupledElement< ELEMENT >::fill_in_contribution_to_jacobian(), and oomph::ScaleCoupledElement< ELEMENT >::fill_in_contribution_to_residuals().

◆ get_coupling_residual()

template<class ELEMENT >
double oomph::ScaleCoupledElement< ELEMENT >::get_coupling_residual ( const unsigned &  n,
const unsigned &  d 
)
inline

Returns the coupling residual.

52  {
53  return coupling_residual.empty()?0.0:coupling_residual[n][d];
54  }
const unsigned n
Definition: CG3DPackingUnitTest.cpp:32

References oomph::ScaleCoupledElement< ELEMENT >::coupling_residual, and n.

◆ get_coupling_weight()

template<class ELEMENT >
double oomph::ScaleCoupledElement< ELEMENT >::get_coupling_weight ( const unsigned &  n)
inline

Sets the coupling weight.

72  {
73  return coupling_weight.empty()?0.0:coupling_weight[n];
74  }

References oomph::ScaleCoupledElement< ELEMENT >::coupling_weight, and n.

◆ set_coupling_jacobian()

template<class ELEMENT >
void oomph::ScaleCoupledElement< ELEMENT >::set_coupling_jacobian ( Vector< Vector< double >> &  jacobian)
inline

Sets the coupling residual.

62  {
63  coupling_jacobian = jacobian;
64  }

References oomph::ScaleCoupledElement< ELEMENT >::coupling_jacobian.

◆ set_coupling_residual()

template<class ELEMENT >
void oomph::ScaleCoupledElement< ELEMENT >::set_coupling_residual ( Vector< Vector< double >> &  residual)
inline

Sets the coupling residual.

46  {
47  coupling_residual = residual;
48  }

References oomph::ScaleCoupledElement< ELEMENT >::coupling_residual.

◆ set_coupling_weight()

template<class ELEMENT >
void oomph::ScaleCoupledElement< ELEMENT >::set_coupling_weight ( Vector< double > &  weight)
inline

Sets the coupling weight.

67  {
68  coupling_weight = weight;
69  }

References oomph::ScaleCoupledElement< ELEMENT >::coupling_weight.

Member Data Documentation

◆ coupling_jacobian

template<class ELEMENT >
Vector<Vector<double> > oomph::ScaleCoupledElement< ELEMENT >::coupling_jacobian
private

◆ coupling_residual

◆ coupling_weight

template<class ELEMENT >
Vector<double> oomph::ScaleCoupledElement< ELEMENT >::coupling_weight
private

The documentation for this class was generated from the following file: