/[escript]/trunk/cusplibrary/cusp/gallery/diffusion.h
ViewVC logotype

Contents of /trunk/cusplibrary/cusp/gallery/diffusion.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5148 - (show annotations)
Mon Sep 15 01:25:23 2014 UTC (5 years, 9 months ago) by caltinay
File MIME type: text/plain
File size: 3484 byte(s)
Merging ripley diagonal storage + CUDA support into trunk.
Options file version has been incremented due to new options
'cuda' and 'nvccflags'.

1 /*
2 * Copyright 2008-2009 NVIDIA Corporation
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17 /*! \file diffusion.h
18 * \brief Anisotropic diffusion matrix generators
19 */
20
21 #pragma once
22
23 #include <cusp/detail/config.h>
24 #include <cusp/gallery/stencil.h>
25
26 #ifdef _WIN32
27 #define _USE_MATH_DEFINES 1 // make sure M_PI is defined
28 #endif
29 #include <math.h>
30
31 namespace cusp
32 {
33
34 namespace gallery
35 {
36
37 /*! \addtogroup gallery Matrix Gallery
38 * \ingroup gallery
39 * \{
40 */
41
42 struct disc_type {};
43
44 struct FD : public disc_type {};
45 struct FE : public disc_type {};
46
47 /*! \p diffusion: Create a matrix representing an anisotropic
48 * Poisson problem discretized on an \p m by \p n grid with
49 * the a given direction.
50 *
51 * \param matrix output
52 * \param m number of grid rows
53 * \param n number of grid columns
54 * \param eps magnitude of anisotropy
55 * \param theta angle of anisotropy in radians
56 *
57 * \tparam MatrixType matrix container
58 *
59 */
60 template <typename Method, typename MatrixType>
61 void diffusion( MatrixType& matrix, size_t m, size_t n,
62 const double eps = 1e-5,
63 const double theta = M_PI/4.0)
64 {
65 typedef typename MatrixType::index_type IndexType;
66 typedef typename MatrixType::value_type ValueType;
67 typedef thrust::tuple<IndexType,IndexType> StencilIndex;
68 typedef thrust::tuple<StencilIndex,ValueType> StencilPoint;
69
70 ValueType C = cos(theta);
71 ValueType S = sin(theta);
72 ValueType CC = C*C;
73 ValueType SS = S*S;
74 ValueType CS = C*S;
75
76 ValueType a;
77 ValueType b;
78 ValueType c;
79 ValueType d;
80 ValueType e;
81
82 if( thrust::detail::is_same<Method, FE>::value )
83 {
84 a = (-1.0*eps - 1.0)*CC + (-1.0*eps - 1.0)*SS + ( 3.0*eps - 3.0)*CS;
85 b = ( 2.0*eps - 4.0)*CC + (-4.0*eps + 2.0)*SS;
86 c = (-1.0*eps - 1.0)*CC + (-1.0*eps - 1.0)*SS + (-3.0*eps + 3.0)*CS;
87 d = (-4.0*eps + 2.0)*CC + ( 2.0*eps - 4.0)*SS;
88 e = ( 8.0*eps + 8.0)*CC + ( 8.0*eps + 8.0)*SS;
89
90 a /= 6.0;
91 b /= 6.0;
92 c /= 6.0;
93 d /= 6.0;
94 e /= 6.0;
95 }
96 else if( thrust::detail::is_same<Method, FD>::value )
97 {
98 a = 0.5 * (eps-1.0) * CS;
99 b = -(eps*SS + CC);
100 c = -a;
101 d = -(eps*CC + SS);
102 e = 2.0 * (eps+1.0);
103 }
104 else
105 {
106 throw cusp::invalid_input_exception("unrecognized discretization method");
107 }
108
109 cusp::array1d<StencilPoint, cusp::host_memory> stencil;
110
111 stencil.push_back(StencilPoint(StencilIndex( -1, -1), a));
112 stencil.push_back(StencilPoint(StencilIndex( 0, -1), b));
113 stencil.push_back(StencilPoint(StencilIndex( 1, -1), c));
114 stencil.push_back(StencilPoint(StencilIndex( -1, 0), d));
115 stencil.push_back(StencilPoint(StencilIndex( 0, 0), e));
116 stencil.push_back(StencilPoint(StencilIndex( 1, 0), d));
117 stencil.push_back(StencilPoint(StencilIndex( -1, 1), c));
118 stencil.push_back(StencilPoint(StencilIndex( 0, 1), b));
119 stencil.push_back(StencilPoint(StencilIndex( 1, 1), a));
120
121 cusp::gallery::generate_matrix_from_stencil(matrix, stencil, StencilIndex(m,n));
122 }
123 /*! \}
124 */
125
126 } // end namespace gallery
127 } // end namespace cusp
128

  ViewVC Help
Powered by ViewVC 1.1.26