Skip to content

Commit c241be8

Browse files
jacobmersonjacobmerson
committed
compute the derivative of a matrix field
closes #201
1 parent 872aa0a commit c241be8

File tree

6 files changed

+72
-5
lines changed

6 files changed

+72
-5
lines changed

‎apf/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ set(SOURCES
1616
apfDynamicMatrix.cc
1717
apfDynamicVector.cc
1818
apfMatrixField.cc
19+
apfMatrixElement.cc
1920
apfMesh.cc
2021
apfMesh2.cc
2122
apfMigrate.cc

‎apf/apf.cc‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -278,6 +278,13 @@ void getMatrix(Element* e, Vector3 const& param, Matrix3x3& value)
278278
value = element->getValue(param);
279279
}
280280

281+
282+
void getMatrixGrad(Element* e, Vector3 const& param, Vector<27>& deriv)
283+
{
284+
MatrixElement* element = static_cast<MatrixElement*>(e);
285+
return element->grad(param,deriv);
286+
}
287+
281288
void getComponents(Element* e, Vector3 const& param, double* components)
282289
{
283290
e->getComponents(param,components);

‎apf/apf.h‎

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -339,8 +339,8 @@ void getCurl(Element* e, Vector3 const& param, Vector3& curl);
339339
* \param param The local coordinates in the element.
340340
* \param deriv The gradient matrix at that point.
341341
* \details Note: The return parameter component deriv[j][i] stores the
342-
* value (grad u)_{i,j} = \frac{\partial u_i} / \frac{\partial x_j},
343-
* where u is the vector field of interest, i is the vector index,
342+
* value $(grad u)_{i,j} = \frac{\partial u_i} / \frac{\partial x_j}$,
343+
* where $u$ is the vector field of interest, i is the vector index,
344344
* and j is the derivative index.
345345
*/
346346
void getVectorGrad(Element* e, Vector3 const& param, Matrix3x3& deriv);
@@ -352,6 +352,13 @@ void getVectorGrad(Element* e, Vector3 const& param, Matrix3x3& deriv);
352352
*/
353353
void getMatrix(Element* e, Vector3 const& param, Matrix3x3& value);
354354

355+
/** \brief get the gradient of a matrix field
356+
* \param param The local coordinate in the element.
357+
* \param value the gradient of the field at a point. If the matrix is defined by A_{ij},
358+
* the gradient $\frac{\partial A_{ij}}{\partial X_k}=value[i*3+j+9*k]$ .
359+
*/
360+
void getMatrixGrad(Element* e, Vector3 const& param, Vector<27>& value);
361+
355362
/** \brief Evaluate a field into an array of component values. */
356363
void getComponents(Element* e, Vector3 const& param, double* components);
357364

‎apf/apfMatrixElement.cc‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
/*
2+
* Copyright 2018 Scientific Computation Research Center
3+
*
4+
* This work is open source software, licensed under the terms of the
5+
* BSD license as described in the LICENSE file in the top-level directory.
6+
*/
7+
8+
#include "apfMatrixElement.h"
9+
#include "apfVectorElement.h"
10+
11+
namespace apf {
12+
// laid out in array as F_i*3+j+9*d
13+
void MatrixElement::grad(Vector3 const& xi, Vector<27>& g)
14+
{
15+
Matrix3x3* nodeValues = getNodeValues();
16+
NewArray<Vector3> globalGradients;
17+
getGlobalGradients(xi, globalGradients);
18+
// for the first time through g, set the values of g
19+
for(int i=0; i<3; ++i) {
20+
for(int j=0; j<3; ++j) {
21+
for(int d=0; d<3; ++d) {
22+
g[i*3+j+d*9]= nodeValues[0][i][j]*globalGradients[0][d];
23+
}
24+
}
25+
}
26+
for(int nd=1; nd<nen; ++nd) {
27+
for(int i=0; i<3; ++i) {
28+
for(int j=0; j<3; ++j) {
29+
for(int d=0; d<3; ++d) {
30+
g[i*3+j+d*9]+= nodeValues[nd][i][j]*globalGradients[nd][d];
31+
}
32+
}
33+
}
34+
}
35+
}
36+
}//namespace apf

‎apf/apfMatrixElement.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ class MatrixElement : public ElementOf<Matrix3x3>
1818
{
1919
public:
2020
MatrixElement(MatrixField* f, VectorElement* e);
21+
void grad(Vector3 const& xi, Vector<27>& g);
2122
virtual ~MatrixElement();
2223
};
2324

‎spr/sprGetGradIPField.cc‎

Lines changed: 18 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77

88
#include "spr.h"
99
#include "apfMesh.h"
10+
#include "apfShape.h"
1011
#include <pcu_util.h>
1112

1213
namespace spr {
@@ -16,8 +17,16 @@ apf::Field* getGradIPField(apf::Field* f, const char* name, int order)
1617
PCU_ALWAYS_ASSERT(f);
1718
apf::Mesh* m = getMesh(f);
1819
int vt = apf::getValueType(f);
19-
PCU_ALWAYS_ASSERT(vt == apf::SCALAR || vt == apf::VECTOR);
20-
apf::Field* ip_field = apf::createIPField(m,name,vt+1,order);
20+
PCU_ALWAYS_ASSERT(vt == apf::SCALAR || vt == apf::VECTOR || vt == apf::MATRIX);
21+
apf::Field* ip_field = NULL;
22+
if (vt == apf::MATRIX)
23+
{
24+
ip_field = apf::createPackedField(m, name, 27, apf::getIPShape(m->getDimension(), order));
25+
}
26+
else
27+
{
28+
ip_field = apf::createIPField(m,name,vt+1,order);
29+
}
2130
apf::MeshIterator* it = m->begin(m->getDimension());
2231
apf::MeshEntity* e;
2332
while ((e = m->iterate(it)))
@@ -35,12 +44,18 @@ apf::Field* getGradIPField(apf::Field* f, const char* name, int order)
3544
apf::getGrad(fe,xi,value);
3645
apf::setVector(ip_field,e,p,value);
3746
}
38-
else
47+
else if (vt == apf::VECTOR)
3948
{
4049
apf::Matrix3x3 value;
4150
apf::getVectorGrad(fe,xi,value);
4251
apf::setMatrix(ip_field,e,p,value);
4352
}
53+
else if (vt == apf::MATRIX)
54+
{
55+
apf::Vector<27> value;
56+
apf::getMatrixGrad(fe, xi, value);
57+
apf::setComponents(ip_field, e, p, &value[0]);
58+
}
4459
}
4560
apf::destroyElement(fe);
4661
apf::destroyMeshElement(me);

0 commit comments

Comments
 (0)