Newer
Older
#ifndef AVERAGE_INTERFACE_HH
#define AVERAGE_INTERFACE_HH
#include <dune/common/fmatrix.hh>
#include "svd.hh"
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
// Given a resultant force and torque (from a rod problem), this method computes the corresponding
// Neumann data for a 3d elasticity problem.
template <class GridType>
void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>& resultantForce,
const Dune::FieldVector<double,GridType::dimension>& resultantTorque,
const BoundaryPatch<GridType>& interface,
const Configuration& crossSection,
Dune::BlockVector<Dune::FieldVector<double, GridType::dimension> >& pressure)
{
const GridType& grid = interface.getGrid();
const int level = interface.level();
const typename GridType::Traits::LevelIndexSet& indexSet = grid.levelIndexSet(level);
const int dim = GridType::dimension;
typedef typename GridType::ctype ctype;
typedef typename GridType::template Codim<dim>::LevelIterator VertexIterator;
// set up output array
pressure.resize(indexSet.size(dim));
pressure = 0;
ctype area = interface.area();
VertexIterator vIt = indexSet.template begin<dim, Dune::All_Partition>();
VertexIterator vEndIt = indexSet.template end<dim, Dune::All_Partition>();
for (; vIt!=vEndIt; ++vIt) {
int vIdx = indexSet.index(*vIt);
if (interface.containsVertex(vIdx)) {
// force part
pressure[vIdx] = resultantForce;
pressure[vIdx] /= area;
// torque part
double x = (vIt->geometry()[0] - crossSection.r) * crossSection.q.director(0);
double y = (vIt->geometry()[0] - crossSection.r) * crossSection.q.director(1);
Dune::FieldVector<double,3> localTorque;
for (int i=0; i<3; i++)
localTorque[i] = resultantTorque * crossSection.q.director(i);
// add it up
pressure[vIdx][0] += -2 * M_PI * localTorque[2] * y / (area*area);
pressure[vIdx][1] += 2 * M_PI * localTorque[2] * x / (area*area);
pressure[vIdx][2] += 4 * M_PI * localTorque[0] * y / (area*area);
pressure[vIdx][2] += -4 * M_PI * localTorque[1] * x / (area*area);
}
}
}
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
template <class GridType>
void computeAverageInterface(const BoundaryPatch<GridType>& interface,
const Dune::BlockVector<Dune::FieldVector<double,GridType::dimension> > deformation,
Configuration& average)
{
using namespace Dune;
typedef typename GridType::template Codim<0>::LevelIterator ElementIterator;
typedef typename GridType::template Codim<0>::Entity EntityType;
typedef typename EntityType::LevelIntersectionIterator NeighborIterator;
const GridType& grid = interface.getGrid();
const int level = interface.level();
const typename GridType::Traits::LevelIndexSet& indexSet = grid.levelIndexSet(level);
const int dim = GridType::dimension;
// ///////////////////////////////////////////
// Initialize output configuration
// ///////////////////////////////////////////
average.r = 0;
double interfaceArea = 0;
FieldMatrix<double,dim,dim> deformationGradient(0);
// ///////////////////////////////////////////
// Loop and integrate over the interface
// ///////////////////////////////////////////
ElementIterator eIt = grid.template lbegin<0>(level);
ElementIterator eEndIt = grid.template lend<0>(level);
for (; eIt!=eEndIt; ++eIt) {
NeighborIterator nIt = eIt->ilevelbegin();
NeighborIterator nEndIt = eIt->ilevelend();
for (; nIt!=nEndIt; ++nIt) {
if (!interface.contains(*eIt, nIt))
continue;
const typename NeighborIterator::Geometry& segmentGeometry = nIt.intersectionGlobal();
const ReferenceElement<double,dim>& refElement = ReferenceElements<double, dim>::general(eIt->geometry().type());
int nDofs = refElement.size(nIt.numberInSelf(),1,dim);
// Get quadrature rule
const QuadratureRule<double, dim-1>& quad = QuadratureRules<double, dim-1>::rule(segmentGeometry.type(), dim-1);
// Get set of shape functions on this segment
const typename LagrangeShapeFunctionSetContainer<double,double,dim>::value_type& sfs
= LagrangeShapeFunctions<double,double,dim>::general(eIt->geometry().type(),1);
/* Loop over all integration points */
for (int ip=0; ip<quad.size(); ip++) {
// Local position of the quadrature point
//const FieldVector<double,dim-1>& quadPos = quad[ip].position();
const FieldVector<double,dim> quadPos = nIt.intersectionSelfLocal().global(quad[ip].position());
const double integrationElement = segmentGeometry.integrationElement(quad[ip].position());
// Evaluate base functions
FieldVector<double,dim> posAtQuadPoint(0);
for(int i=0; i<eIt->geometry().corners(); i++) {
int idx = indexSet.template subIndex<dim>(*eIt, i);
// Deformation at the quadrature point
posAtQuadPoint.axpy(sfs[i].evaluateFunction(0,quadPos), deformation[idx]);
}
const FieldMatrix<double,dim,dim>& inv = eIt->geometry().jacobianInverseTransposed(quadPos);
/* Compute the weight of the current integration point */
double weight = quad[ip].weight() * integrationElement;
/**********************************************/
/* compute gradients of the shape functions */
/**********************************************/
std::vector<FieldVector<double, dim> > shapeGrads(eIt->geometry().corners());
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
for (int dof=0; dof<eIt->geometry().corners(); dof++) {
for (int i=0; i<dim; i++)
shapeGrads[dof][i] = sfs[dof].evaluateDerivative(0, i, quadPos);
// multiply with jacobian inverse
FieldVector<double,dim> tmp(0);
inv.umv(shapeGrads[dof], tmp);
shapeGrads[dof] = tmp;
//std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
}
/****************************************************/
// The deformation gradient of the deformation
// in formulas: F(i,j) = $\partial \phi_i / \partial x_j$
// or F(i,j) = Id + $\partial u_i / \partial x_j$
/****************************************************/
FieldMatrix<double, dim, dim> F;
for (int i=0; i<dim; i++) {
for (int j=0; j<dim; j++) {
F[i][j] = (i==j) ? 1 : 0;
for (int k=0; k<eIt->geometry().corners(); k++)
F[i][j] += deformation[indexSet.template subIndex<dim>(*eIt, k)][i]*shapeGrads[k][j];
}
}
interfaceArea += quad[ip].weight() * integrationElement;
average.r.axpy(quad[ip].weight() * integrationElement, posAtQuadPoint);
F *= quad[ip].weight();
deformationGradient += F;
}
}
}
// average deformation of the interface is the integral over its
// deformation divided by its area
average.r /= interfaceArea;
// average deformation is the integral over the deformation gradient
// divided by its area
deformationGradient /= interfaceArea;
std::cout << "deformationGradient: " << std::endl << deformationGradient << std::endl;
// Get the rotational part of the deformation gradient by performing a
// polar composition.
FieldVector<double,dim> W;
FieldMatrix<double,dim,dim> VT;
svdcmp<double,dim,dim>(deformationGradient, W, VT);
deformationGradient.rightmultiply(VT);
assert( std::abs(1-deformationGradient.determinant()) < 1e-3);
//std::cout << "determinant: " << deformationGradient.determinant() << " (should be 1)\n";
// average orientation not implemented yet
average.q = Quaternion<double>::identity();
}
#endif