| 1 | /*!\file: GetGlobalDofList.cpp
|
|---|
| 2 | * \brief create transform matrix for different coordinate systems
|
|---|
| 3 | */
|
|---|
| 4 | #include "./elements.h"
|
|---|
| 5 | #include <math.h>
|
|---|
| 6 |
|
|---|
| 7 | void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array){
|
|---|
| 8 |
|
|---|
| 9 | int i,counter;
|
|---|
| 10 | int numdofs = 0;
|
|---|
| 11 | IssmDouble norm;
|
|---|
| 12 | IssmDouble *transform = NULL;
|
|---|
| 13 | IssmDouble *values = NULL;
|
|---|
| 14 | IssmDouble coord_system[3][3];
|
|---|
| 15 |
|
|---|
| 16 | /*Some checks in debugging mode*/
|
|---|
| 17 | _assert_(numnodes && nodes);
|
|---|
| 18 |
|
|---|
| 19 | /*Get total number of dofs*/
|
|---|
| 20 | for(i=0;i<numnodes;i++){
|
|---|
| 21 | switch(cs_array[i]){
|
|---|
| 22 | case XYEnum: numdofs+=2; break;
|
|---|
| 23 | case XYZPEnum: numdofs+=4; break;
|
|---|
| 24 | default: _error_("Coordinate system %s not supported yet",EnumToStringx(cs_array[i]));
|
|---|
| 25 | }
|
|---|
| 26 | }
|
|---|
| 27 |
|
|---|
| 28 | /*Allocate and initialize transform matrix*/
|
|---|
| 29 | transform=xNew<IssmDouble>(numdofs*numdofs);
|
|---|
| 30 | for(i=0;i<numdofs*numdofs;i++) transform[i]=0.0;
|
|---|
| 31 |
|
|---|
| 32 | /*Create transform matrix for all nodes (x,y for 2d and x,y,z for 3d). It is a block matrix
|
|---|
| 33 | *for 3 nodes:
|
|---|
| 34 |
|
|---|
| 35 | * | T1 0 0 |
|
|---|
| 36 | * Q = | 0 T2 0 |
|
|---|
| 37 | * | 0 0 T3|
|
|---|
| 38 | *
|
|---|
| 39 | * Where T1 is the transform matrix for node 1. It is a simple copy of the coordinate system
|
|---|
| 40 | * associated to this node*/
|
|---|
| 41 | counter=0;
|
|---|
| 42 | for(i=0;i<numnodes;i++){
|
|---|
| 43 | nodes[i]->GetCoordinateSystem(&coord_system[0][0]);
|
|---|
| 44 | switch(cs_array[i]){
|
|---|
| 45 | case XYEnum:
|
|---|
| 46 | /*We remove the z component, we need to renormalize x and y: x=[x1 x2 0] y=[-x2 x1 0]*/
|
|---|
| 47 | norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]); _assert_(norm>1.e-4);
|
|---|
| 48 | transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0]/norm;
|
|---|
| 49 | transform[(numdofs)*(counter+0) + counter+1] = - coord_system[1][0]/norm;
|
|---|
| 50 | transform[(numdofs)*(counter+1) + counter+0] = coord_system[1][0]/norm;
|
|---|
| 51 | transform[(numdofs)*(counter+1) + counter+1] = coord_system[0][0]/norm;
|
|---|
| 52 | counter+=2;
|
|---|
| 53 | break;
|
|---|
| 54 | case XYZPEnum:
|
|---|
| 55 | /*Only the first 3 coordinates are changed (x,y,z), leave the others (P) unchanged*/
|
|---|
| 56 | transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0];
|
|---|
| 57 | transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1];
|
|---|
| 58 | transform[(numdofs)*(counter+0) + counter+2] = coord_system[0][2];
|
|---|
| 59 | transform[(numdofs)*(counter+1) + counter+0] = coord_system[1][0];
|
|---|
| 60 | transform[(numdofs)*(counter+1) + counter+1] = coord_system[1][1];
|
|---|
| 61 | transform[(numdofs)*(counter+1) + counter+2] = coord_system[1][2];
|
|---|
| 62 | transform[(numdofs)*(counter+2) + counter+0] = coord_system[2][0];
|
|---|
| 63 | transform[(numdofs)*(counter+2) + counter+1] = coord_system[2][1];
|
|---|
| 64 | transform[(numdofs)*(counter+2) + counter+2] = coord_system[2][2];
|
|---|
| 65 | transform[(numdofs)*(counter+3) + counter+3] = 1.0;
|
|---|
| 66 | counter+=4;
|
|---|
| 67 | break;
|
|---|
| 68 | default:
|
|---|
| 69 | _error_("Coordinate system %s not supported yet",EnumToStringx(cs_array[i]));
|
|---|
| 70 | }
|
|---|
| 71 | }
|
|---|
| 72 |
|
|---|
| 73 | /*Assign output pointer*/
|
|---|
| 74 | *ptransform=transform;
|
|---|
| 75 | }
|
|---|