v0.14.0
SphereSurfaceIntegration.cpp
Go to the documentation of this file.
1 /**
2  * \fiele SphereSurfaceIntegration.cpp
3  *
4  * \brief Integration points are taken from OOFEM
5  * <http://www.oofem.org/resources/doc/oofemrefman/microplanematerial_8C_source.html>
6  *
7  */
8 
9 /*
10 * This file is part of MoFEM.
11 * MoFEM is free software: you can redistribute it and/or modify it under
12 * the terms of the GNU Lesser General Public License as published by the
13 * Free Software Foundation, either version 3 of the License, or (at your
14 * option) any later version.
15 *
16 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19 * License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #include <MoFEM.hpp>
25 using namespace MoFEM;
26 
28 using namespace BoneRemodeling;
29 
31  MatrixDouble &gauss_pts,VectorDouble &weights
32 ) {
33  PetscFunctionBegin;
34  gauss_pts.resize(28,3,false);
35  weights.resize(28,false);
36  for(int i = 0; i != 4; i++ ) {
37  weights[ i ] = 0.0160714276E0;
38  weights[ i + 4 ] = 0.0204744730E0;
39  weights[ i + 8 ] = 0.0204744730E0;
40  weights[ i + 12 ] = 0.0204744730E0;
41  weights[ i + 16 ] = 0.0158350505E0;
42  weights[ i + 20 ] = 0.0158350505E0;
43  weights[ i + 24 ] = 0.0158350505E0;
44  }
45  double help [ 7 ] [ 3 ] = {
46  { .577350259E0, .577350259E0, .577350259E0 }, { .935113132E0, .250562787E0, .250562787E0 },
47  { .250562787E0, .935113132E0, .250562787E0 }, { .250562787E0, .250562787E0, .935113132E0 },
48  { .186156720E0, .694746614E0, .694746614E0 }, { .694746614E0, .186156720E0, .694746614E0 },
49  { .694746614E0, .694746614E0, .186156720E0 }
50  };
51  for(int i = 0; i != 7; i++ ) {
52  i4 = i * 4;
53  gauss_pts(i4,0) = help [ i ] [ 0 ];
54  gauss_pts(i4,1) = help [ i ] [ 1 ];
55  gauss_pts(i4,2) = help [ i ] [ 2 ];
56  gauss_pts(i4 + 1,0) = help [ i ] [ 0 ];
57  gauss_pts(i4 + 1,1) = help [ i ] [ 1 ];
58  gauss_pts(i4 + 1,2) = -help [ i ] [ 2 ];
59  gauss_pts(i4 + 2,0) = help [ i ] [ 0 ];
60  gauss_pts(i4 + 2,1) = -help [ i ] [ 1 ];
61  gauss_pts(i4 + 2,2) = help [ i ] [ 2 ];
62  gauss_pts(i4 + 3,0) = help [ i ] [ 0 ];
63  gauss_pts(i4 + 3,1) = -help [ i ] [ 1 ];
64  gauss_pts(i4 + 3,2) = -help [ i ] [ 2 ];
65  }
66  PetscFunctionReturn(0);
67 }
68 
70  MatrixDouble &gauss_pts,VectorDouble &weights
71 ) {
72  PetscFunctionBegin;
73  weights.resize(21,false);
74  gauss_pts.resize(21,3,false);
75  for(int i = 0; i != 3; i++ ) {
76  weights[ i ] = 0.02652141274;
77  }
78  for(int i = 3; i != 9; i++ ) {
79  weights[ i ] = 0.01993014153;
80  }
81  for(int i = 9; i != 21; i++ ) {
82  weights[ i ] = 0.02507124272;
83  }
84  gauss_pts(0,0) = gauss_pts(1,1) = gauss_pts(2,2) = 1.;
85  gauss_pts(0,1) = gauss_pts(0,2) = gauss_pts(1,0) = gauss_pts(1,2) = gauss_pts(2,0) = gauss_pts(2,1) = 0.;
86  gauss_pts(3,0) = gauss_pts(3,1) = gauss_pts(4,0) = gauss_pts(5,0) = gauss_pts(5,2) = gauss_pts(6,0) = gauss_pts(7,1) = gauss_pts(7,2) = gauss_pts(8,1) = 0.7071067812;
87  gauss_pts(4,1) = gauss_pts(6,2) = gauss_pts(8,2) = -0.7071067812;
88  gauss_pts(3,2) = gauss_pts(4,2) = gauss_pts(5,1) = gauss_pts(6,1) = gauss_pts(7,0) = gauss_pts(8,0) = 0.;
89  double help [ 3 ] [ 3 ] =
90  { { 0.3879072746, 0.3879072746, 0.8360956240 },
91  { 0.3879072746, 0.8360956240, 0.3879072746 },
92  { 0.8360956240, 0.3879072746, 0.3879072746 } };
93 
94  for ( i = 0; i < 3; i++ ) {
95  i4 = i * 4;
96  gauss_pts(9 + i4,0) = help [ i ] [ 0 ];
97  gauss_pts(9 + i4,1) = help [ i ] [ 1 ];
98  gauss_pts(9 + i4,2) = help [ i ] [ 2 ];
99 
100  gauss_pts(10 + i4,0) = help [ i ] [ 0 ];
101  gauss_pts(10 + i4,1) = help [ i ] [ 1 ];
102  gauss_pts(10 + i4,2) = -help [ i ] [ 2 ];
103 
104  gauss_pts(11 + i4,0) = help [ i ] [ 0 ];
105  gauss_pts(11 + i4,1) = -help [ i ] [ 1 ];
106  gauss_pts(11 + i4,2) = help [ i ] [ 2 ];
107 
108  gauss_pts(12 + i4,0) = help [ i ] [ 0 ];
109  gauss_pts(12 + i4,1) = -help [ i ] [ 1 ];
110  gauss_pts(12 + i4,2) = -help [ i ] [ 2 ];
111  }
112  PetscFunctionReturn(0);
113 }
114 
116  MatrixDouble &gauss_pts,VectorDouble &weights
117 ) {
118  PetscFunctionBegin;
119  double help [ 61 ] [ 4 ] = {
120  { 1.000000000000, 0.000000000000, 0.000000000000, 0.00795844204678 },
121  { 0.745355992500, 0.0, 0.666666666667, 0.00795844204678 },
122  { 0.745355992500, -0.577350269190, -0.333333333333, 0.00795844204678 },
123  { 0.745355992500, 0.577350269190, -0.333333333333, 0.00795844204678 },
124  { 0.333333333333, 0.577350269190, 0.745355992500, 0.00795844204678 },
125  { 0.333333333333, -0.577350269190, 0.745355992500, 0.00795844204678 },
126  { 0.333333333333, -0.934172358963, 0.127322003750, 0.00795844204678 },
127  { 0.333333333333, -0.356822089773, -0.872677996250, 0.00795844204678 },
128  { 0.333333333333, 0.356822089773, -0.872677996250, 0.00795844204678 },
129  { 0.333333333333, 0.934172358963, 0.127322003750, 0.00795844204678 },
130  { 0.794654472292, -0.525731112119, 0.303530999103, 0.0105155242892 },
131  { 0.794654472292, 0.0, -0.607061998207, 0.0105155242892 },
132  { 0.794654472292, 0.525731112119, 0.303530999103, 0.0105155242892 },
133  { 0.187592474085, 0.0, 0.982246946377, 0.0105155242892 },
134  { 0.187592474085, -0.850650808352, -0.491123473188, 0.0105155242892 },
135  { 0.187592474085, 0.850650808352, -0.491123473188, 0.0105155242892 },
136  { 0.934172358963, 0.0, 0.356822089773, 0.0100119364272 },
137  { 0.934172358963, -0.309016994375, -0.178411044887, 0.0100119364272 },
138  { 0.934172358963, 0.309016994375, -0.178411044887, 0.0100119364272 },
139  { 0.577350269190, 0.309016994375, 0.755761314076, 0.0100119364272 },
140  { 0.577350269190, -0.309016994375, 0.755761314076, 0.0100119364272 },
141  { 0.577350269190, -0.809016994375, -0.110264089708, 0.0100119364272 },
142  { 0.577350269190, -0.5, -0.645497224368, 0.0100119364272 },
143  { 0.577350269190, 0.5, -0.645497224368, 0.0100119364272 },
144  { 0.577350269190, 0.809016994375, -0.110264089708, 0.0100119364272 },
145  { 0.356822089773, -0.809016994375, 0.467086179481, 0.0100119364272 },
146  { 0.356822089773, 0.0, -0.934172358963, 0.0100119364272 },
147  { 0.356822089773, 0.809016994375, 0.467086179481, 0.0100119364272 },
148  { 0.0, 0.5, 0.866025403784, 0.0100119364272 },
149  { 0.0, -1.0, 0.0, 0.0100119364272 },
150  { 0.0, 0.5, -0.866025403784, 0.0100119364272 },
151  { 0.947273580412, -0.277496978165, 0.160212955043, 0.00690477957966 },
152  { 0.812864676392, -0.277496978165, 0.512100034157, 0.00690477957966 },
153  { 0.595386501297, -0.582240127941, 0.553634669695, 0.00690477957966 },
154  { 0.595386501297, -0.770581752342, 0.227417407053, 0.00690477957966 },
155  { 0.812864676392, -0.582240127941, -0.015730584514, 0.00690477957966 },
156  { 0.492438766306, -0.753742692223, -0.435173546254, 0.00690477957966 },
157  { 0.274960591212, -0.942084316623, -0.192025554687, 0.00690477957966 },
158  { -0.076926487903, -0.942084316623, -0.326434458707, 0.00690477957966 },
159  { -0.076926487903, -0.753742692223, -0.652651721349, 0.00690477957966 },
160  { 0.274960591212, -0.637341166847, -0.719856173359, 0.00690477957966 },
161  { 0.947273580412, 0.0, -0.320425910085, 0.00690477957966 },
162  { 0.812864676392, -0.304743149777, -0.496369449643, 0.00690477957966 },
163  { 0.595386501297, -0.188341624401, -0.781052076747, 0.00690477957966 },
164  { 0.595386501297, 0.188341624401, -0.781052076747, 0.00690477957966 },
165  { 0.812864676392, 0.304743149777, -0.496369449643, 0.00690477957966 },
166  { 0.492438766306, 0.753742692223, -0.435173546254, 0.00690477957966 },
167  { 0.274960591212, 0.637341166847, -0.719856173359, 0.00690477957966 },
168  { -0.076926487903, 0.753742692223, -0.652651721349, 0.00690477957966 },
169  { -0.076926487903, 0.942084316623, -0.326434458707, 0.00690477957966 },
170  { 0.274960591212, 0.942084316623, -0.192025554687, 0.00690477957966 },
171  { 0.947273580412, 0.277496978165, 0.160212955043, 0.00690477957966 },
172  { 0.812864676392, 0.582240127941, -0.015730584514, 0.00690477957966 },
173  { 0.595386501297, 0.770581752342, 0.227417407053, 0.00690477957966 },
174  { 0.595386501297, 0.582240127941, 0.553634669695, 0.00690477957966 },
175  { 0.812864676392, 0.277496978165, 0.512100034157, 0.00690477957966 },
176  { 0.492438766306, 0.0, 0.870347092509, 0.00690477957966 },
177  { 0.274960591212, 0.304743149777, 0.911881728046, 0.00690477957966 },
178  { -0.076926487903, 0.188341624401, 0.979086180056, 0.00690477957966 },
179  { -0.076926487903, -0.188341624401, 0.979086180056, 0.00690477957966 },
180  { 0.274960591212, -0.304743149777, 0.911881728046, 0.00690477957966 }
181  };
182  for ( i = 0; i < numberOfMicroplanes; i++ ) {
183  gauss_pts(i,0 ) = help [ i ] [ 0 ];
184  gauss_pts(i,1 ) = help [ i ] [ 1 ];
185  gauss_pts(i,2 ) = help [ i ] [ 2 ];
186  weights[i] = help [ i ] [ 3 ];
187  }
188  PetscFunctionReturn(0);
189 }
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM.hpp
BoneRemodeling::sphereSurfaceIntegration61
PetscErrorCode sphereSurfaceIntegration61(MatrixDouble &gauss_pts, VectorDouble &weights)
61 integration points on sphere
Definition: SphereSurfaceIntegration.cpp:115
SphereSurfaceIntegration.hpp
Quadrature points on sphere surface.
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
BoneRemodeling
Definition: DensityMaps.hpp:27
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
BoneRemodeling::sphereSurfaceIntegration21
PetscErrorCode sphereSurfaceIntegration21(MatrixDouble &gauss_pts, VectorDouble &weights)
21 integration points on sphere
Definition: SphereSurfaceIntegration.cpp:69
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
BoneRemodeling::sphereSurfaceIntegration28
PetscErrorCode sphereSurfaceIntegration28(MatrixDouble &gauss_pts, VectorDouble &weights)
28 integration points on sphere
Definition: SphereSurfaceIntegration.cpp:30