1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s1308.c,v 1.1 1994-04-21 12:10:42 boh Exp $
45 *
46 */
47
48
49 #define S1308
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1308(double ep[],int idim,double eimpli[],int ideg,double enorm[],int * jstat)55 s1308(double ep[],int idim,double eimpli[],int ideg,double enorm[],int *jstat)
56 #else
57 void s1308(ep,idim,eimpli,ideg,enorm,jstat)
58 double ep[];
59 int idim;
60 double eimpli[];
61 int ideg;
62 double enorm[];
63 int *jstat;
64 #endif
65 /*
66 *********************************************************************
67 *
68 *********************************************************************
69 *
70 * PURPOSE : To calculate the normal vector at point in an implicit
71 * function.
72 *
73 * INPUT : ep - The coordinates of the point
74 * idim - The dimension of the space
75 * eimpli - Description of the implicit surface
76 * ideg - The degree of the implicit surface
77 * ideg=1: Plane
78 * ideg=2; Quadric surface
79 *
80 *
81 * OUTPUT: enorm - The normal vector
82 * jstat - status messages
83 * > 0 : warning
84 * = 0 : ok
85 * < 0 : error
86 *
87 *
88 * METHOD : For degree=1 the 3 first component of eimpli is the
89 * normal vector. For ideg=2: T
90 * The surface can be represented as (X ,1)A(X,1) where A is
91 * a matrix (idim+1xidim+1) representing eimpli, and X a point.
92 * Assume that X is parameterized with respect to some parameters
93 * X = X(s,...). If we take the derivative of X with respect
94 * to s, we get a tangent in the surface.
95 *
96 * dX/ds A X = 0. Thus the expression AX is a normal to
97 * the implicit function. The idim first components of
98 * AX is the acutal normal
99 *
100 *
101 * REFERENCES :
102 *
103 *-
104 * CALLS :
105 * WRITTEN BY : Tor Dokken, SI, Oslo, Norway, 3. July 1988
106 *
107 *********************************************************************
108 */
109 {
110 int ki,kj,kl; /* Variables in loop */
111 int kdimp1=idim+1; /* Dimension + 1 */
112 int kpos=0; /* Position of error */
113 int kstat=0; /* Local error */
114 double tsum; /* Dummy variable */
115
116 if (ideg != 1 && ideg !=2 && ideg != 1001) goto err175;
117
118 if (ideg == 1)
119 {
120 /* First degree implicit surface normal vector is eimpli[0:idim-1] */
121 memcopy(enorm,eimpli,idim,DOUBLE);
122 }
123 else if (ideg==2)
124 {
125
126 /* Calculate the matrix product */
127
128 for (ki=0;ki<idim;ki++)
129 {
130 tsum = eimpli[idim*kdimp1+ki];
131 for (kj=0,kl=ki ; kj<idim ; kj++,kl+=kdimp1)
132 {
133 tsum +=(eimpli[kl]*ep[kj]);
134 }
135 enorm[ki] = tsum;
136 }
137 }
138 else if (ideg==1001)
139 {
140 /* Torus surface */
141
142 double *scentr; /* The center of the torus */
143 double *snorm; /* The normal of the torus symmetry plane */
144 double tbigr; /* The big radius of the torus */
145 double tsmalr; /* The small radius of the torus */
146 double sdum1[3]; /* Temporary storage for point */
147 double sdum2[3]; /* Temporary storage for point */
148 double tproj; /* Projection of vector onto snorm */
149
150
151 scentr = eimpli;
152 snorm = eimpli+3;
153 tbigr = *(eimpli+6);
154 tsmalr = *(eimpli+7);
155
156 /* Find projection of vector from torus center on to torus axis */
157 s6diff(ep,scentr,3,sdum1);
158 tproj = s6scpr(sdum1,snorm,3);
159
160 /* Project vector from torus center to ep onto torus plane */
161 for (ki=0;ki<3;ki++)
162 sdum2[ki] = sdum1[ki] - tproj*snorm[ki];
163 (void)s6norm(sdum2,3,sdum2,&kstat);
164 if (kstat<0) goto error;
165
166 /* Find vector from torus circle to ep */
167 for (ki=0;ki<3;ki++)
168 sdum1[ki] = sdum1[ki] - tbigr*sdum2[ki];
169
170 /* Normalize this vector */
171 (void)s6norm(sdum1,3,enorm,&kstat);
172 if (kstat<0) goto error;
173 }
174
175 *jstat = 0;
176 goto out;
177
178 /* IDEG NOT 1 OR 2 */
179 err175:
180 *jstat = -175;
181 s6err("s1308",*jstat,kpos);
182 goto out;
183
184
185 /* Error in lower leve function */
186 error:
187 *jstat = kstat;
188 s6err("s1308",*jstat,kpos);
189 goto out;
190
191 out:
192 return;
193 }
194
195