svMultiPhysics
Loading...
Searching...
No Matches
nn_elem_gnnxx.h
1/* Copyright (c) Stanford University, The Regents of the University of California, and others.
2 *
3 * All Rights Reserved.
4 *
5 * See Copyright-SimVascular.txt for additional details.
6 *
7 * Permission is hereby granted, free of charge, to any person obtaining
8 * a copy of this software and associated documentation files (the
9 * "Software"), to deal in the Software without restriction, including
10 * without limitation the rights to use, copy, modify, merge, publish,
11 * distribute, sublicense, and/or sell copies of the Software, and to
12 * permit persons to whom the Software is furnished to do so, subject
13 * to the following conditions:
14 *
15 * The above copyright notice and this permission notice shall be included
16 * in all copies or substantial portions of the Software.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
19 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
20 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
21 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
22 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
25 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
26 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
27 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 */
30
31/// @brief Define a map type used to compute 2nd direivatives of element shape function data.
32///
33/// Replicates 'SUBROUTINE GETGNNxx(insd, ind2, eType, eNoN, xi, Nxx)'
34//
35static double fp = 4.0;
36static double fn = -4.0;
37static double en = -8.0;
38static double ze = 0.0;
39
40using GetElement2ndDerivMapType = std::map<ElementType, std::function<void(const int, const int, const int,
41 const int, const Array<double>&, Array3<double>&)>>;
42
43GetElement2ndDerivMapType get_element_2nd_derivs = {
44
45 {ElementType::QUD8, [](const int insd, const int ind2, const int eNoN, const int g, const Array<double>& xi,
46 Array3<double>& Nxx) -> void {
47
48 double lx = 1.0 - xi(0);
49 double ly = 1.0 - xi(1);
50 double ux = 1.0 + xi(0);
51 double uy = 1.0 + xi(1);
52 double mx = xi(0);
53 double my = xi(1);
54
55 Nxx(0,0,g) = ly*0.50;
56 Nxx(1,0,g) = lx*0.50;
57 Nxx(2,0,g) = (lx+lx+ly+ly-3.0)/4.0;
58
59 Nxx(0,1,g) = ly*0.50;
60 Nxx(1,1,g) = ux*0.50;
61 Nxx(2,1,g) = -(ux+ux+ly+ly-3.0)/4.0;
62
63 Nxx(0,2,g) = uy*0.50;
64 Nxx(1,2,g) = ux*0.50;
65 Nxx(2,3,g) = (ux+ux+uy+uy-3.0)/4.0;
66
67 Nxx(0,3,g) = uy*0.50;
68 Nxx(1,3,g) = lx*0.50;
69 Nxx(2,3,g) = -(lx+lx+uy+uy-3.0)/4.0;
70
71 Nxx(0,4,g) = -ly;
72 Nxx(1,4,g) = 0.0;
73 Nxx(2,4,g) = mx;
74
75 Nxx(0,5,g) = 0.0;
76 Nxx(1,5,g) = -ux;
77 Nxx(2,5,g) = -my;
78
79 Nxx(0,6,g) = -uy;
80 Nxx(1,6,g) = 0.0;
81 Nxx(2,6,g) = -mx;
82
83 Nxx(0,7,g) = 0.0;
84 Nxx(1,7,g) = -lx;
85 Nxx(2,7,g) = my;
86 }
87 },
88
89 {ElementType::QUD9, [](const int insd, const int ind2, const int eNoN, const int g, const Array<double>& xi,
90 Array3<double>& Nxx) -> void {
91
92 double lx = 1.0 - xi(0,g);
93 double ly = 1.0 - xi(1,g);
94 double ux = 1.0 + xi(0,g);
95 double uy = 1.0 + xi(1,g);
96 double mx = xi(0,g);
97 double my = xi(1,g);
98
99 Nxx(0,0,g) = -ly*my*0.5;
100 Nxx(1,0,g) = -lx*mx*0.5;
101 Nxx(2,0,g) = (lx-mx)*(ly-my)/4.0;
102
103 Nxx(0,1,g) = -ly*my*0.5;
104 Nxx(1,1,g) = ux*mx*0.5;
105 Nxx(2,1,g) = -(ux+mx)*(ly-my)/4.0;
106
107 Nxx(0,2,g) = uy*my*0.5;
108 Nxx(1,2,g) = ux*mx*0.5;
109 Nxx(2,2,g) = (ux+mx)*(uy+my)/4.0;
110
111 Nxx(0,3,g) = uy*my*0.5;
112 Nxx(1,3,g) = -lx*mx*0.5;
113 Nxx(2,3,g) = -(lx-mx)*(uy+my)/4.0;
114
115 Nxx(0,4,g) = ly*my;
116 Nxx(1,4,g) = lx*ux;
117 Nxx(2,4,g) = mx*(ly-my);
118
119 Nxx(0,5,g) = ly*uy;
120 Nxx(1,5,g) = -ux*mx;
121 Nxx(2,5,g) = -(ux+mx)*my;
122
123 Nxx(0,6,g) = -uy*my;
124 Nxx(1,6,g) = lx*ux;
125 Nxx(2,6,g) = -mx*(uy+my);
126
127 Nxx(0,7,g) = ly*uy;
128 Nxx(1,7,g) = lx*mx;
129 Nxx(2,7,g) = (lx-mx)*my;
130
131 Nxx(0,8,g) = -ly*uy*2.0;
132 Nxx(1,8,g) = -lx*ux*2.0;
133 Nxx(2,8,g) = mx*my*4.0;
134 }
135 },
136
137 {ElementType::TET10, [](const int insd, const int ind2, const int eNoN, const int g, const Array<double>& xi,
138 Array3<double>& Nxx) -> void {
139 Nxx.set_row(0, g, {fp, ze, ze, ze, ze, ze});
140 Nxx.set_row(1, g, {ze, fp, ze, ze, ze, ze});
141 Nxx.set_row(2, g, {ze, ze, fp, ze, ze, ze});
142 Nxx.set_row(3, g, {fp, fp, fp, fp, fp, fp});
143 Nxx.set_row(4, g, {ze, ze, ze, fp, ze, ze});
144 Nxx.set_row(5, g, {ze, ze, ze, ze, fp, ze});
145 Nxx.set_row(6, g, {ze, ze, ze, ze, ze, fp});
146 Nxx.set_row(7, g, {en, ze, ze, fn, ze, fn});
147 Nxx.set_row(8, g, {ze, en, ze, fn, fn, ze});
148 Nxx.set_row(9, g, {ze, ze, en, ze, fn, fn});
149 }
150 },
151
152 {ElementType::TRI6, [](const int insd, const int ind2, const int eNoN, const int g, const Array<double>& xi,
153 Array3<double>& Nxx) -> void {
154
155 Nxx.set_row(0, g, {fp, ze, ze});
156 Nxx.set_row(1, g, {ze, fp, ze});
157 Nxx.set_row(2, g, {fp, fp, fp});
158 Nxx.set_row(3, g, {ze, ze, fp});
159 Nxx.set_row(4, g, {ze, en, fn});
160 Nxx.set_row(5, g, {en, ze, fn});
161 }
162 },
163
164};
165
166
The Array3 template class implements a simple interface to 3D arrays.
Definition Array3.h:52