00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static char const Id[] = "$Id: tettri.c 2681 2009-07-31 04:30:53Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "util.h"
00028
00029 PetscReal triarea(PetscReal *v0,PetscReal *v1,PetscReal *v2,PetscReal *n)
00030 {
00031 PetscReal a[ND],b[ND],outer[ND];
00032 PetscReal result;
00033
00034
00035 my_dcopy(ND,v1,1,a,1);
00036 my_daxpy(ND,-1.0,v0,1,a,1);
00037
00038
00039 my_dcopy(ND,v2,1,b,1);
00040 my_daxpy(ND,-1.0,v0,1,b,1);
00041
00042
00043 douter(ND,a,b,outer);
00044 result=0.5*my_dnrm2(ND,outer,1);
00045
00046 if (n!=PETSC_NULL) {
00047 my_dscal(ND,1.0/(2.0*result),outer,1);
00048 my_dcopy(ND,outer,1,n,1);
00049 }
00050
00051 return(result);
00052 }
00053
00054
00055 PetscReal tetvol(PetscReal *x1,PetscReal *x2,PetscReal *x3,PetscReal *x4)
00056 {
00057 PetscReal result;
00058
00059
00060 result=(1.0/6.0)*
00061 (x3[0]*x2[1]*x1[2] - x4[0]*x2[1]*x1[2] - x2[0]*x3[1]*x1[2] + x4[0]*x3[1]*x1[2] +
00062 x2[0]*x4[1]*x1[2] - x3[0]*x4[1]*x1[2] - x3[0]*x1[1]*x2[2] + x4[0]*x1[1]*x2[2] +
00063 x1[0]*x3[1]*x2[2] - x4[0]*x3[1]*x2[2] - x1[0]*x4[1]*x2[2] + x3[0]*x4[1]*x2[2] +
00064 x2[0]*x1[1]*x3[2] - x4[0]*x1[1]*x3[2] - x1[0]*x2[1]*x3[2] + x4[0]*x2[1]*x3[2] +
00065 x1[0]*x4[1]*x3[2] - x2[0]*x4[1]*x3[2] - x2[0]*x1[1]*x4[2] + x3[0]*x1[1]*x4[2] +
00066 x1[0]*x2[1]*x4[2] - x3[0]*x2[1]*x4[2] - x1[0]*x3[1]*x4[2] + x2[0]*x3[1]*x4[2]
00067 );
00068
00069 return(result);
00070 }
00071
00072 int tetnb(int *tet1, int *tet2, int *vshared)
00073 {
00074 int t1[NV],t2[NV];
00075
00076 for (int i=0;i<NV;i++) {
00077 t1[i]=tet1[i];
00078 t2[i]=tet2[i];
00079 }
00080
00081
00082 PetscSortInt(NV,(PetscInt*)t1);
00083 PetscSortInt(NV,(PetscInt*)t2);
00084
00085 int cntshared=0;
00086
00087 for (int i=0;i<NV;i++) {
00088 for (int j=0;j<NV;j++) {
00089 if (t1[i]==t2[j]) {
00090 if (vshared!=PETSC_NULL) vshared[cntshared]=t1[i];
00091 cntshared++;
00092 }
00093 }
00094 }
00095
00096 assert(cntshared<=NV);
00097
00098 return(cntshared);
00099 }
00100
00101 int TetSharedVertices(int *v1,int *v2,int *s)
00102 {
00103 int k;
00104 k=0;
00105 for (int i=0;i<NV;i++) {
00106 for (int j=0;j<NV;j++) {
00107 if (v1[i]==v2[j]) {
00108 s[k]=v1[i];
00109 k++;
00110 continue;
00111 }
00112 }
00113 }
00114 return(k);
00115 }