2531 {
2532 int debugFn = 0;
2533
2534
2535 if (debugFn) {
2536 for (
int ele = 1; ele <=
NbElements; ++ele) {
2537 printf("ele, Assigned charge: %d, %lg\n", ele,
2538 (
EleArr + ele - 1)->Assigned);
2539 }
2540 }
2541
2542
2543
2544
2545
2546 {
2547 FILE *ChargingUpInpFile = fopen("neBEMChargingUp.inp", "r");
2548 if (ChargingUpInpFile == NULL) {
2549 printf(
2550 "neBEMChargingUp.inp absent ... assuming no charging up effect "
2551 "...\n");
2552
2553 } else {
2554 fscanf(ChargingUpInpFile,
"OptChargingUp: %d\n", &
OptChargingUp);
2556 printf("OptChargingUp = 0 ... assuming no charging up effect ...\n");
2558
2560 char ChargingUpEFile[256];
2561 char ChargingUpIFile[256];
2562 double ChUpFactor;
2563 int *NbChUpEonEle, *NbChUpIonEle;
2564
2565 fscanf(ChargingUpInpFile, "ChargingUpEFile: %255s\n", ChargingUpEFile);
2566 fscanf(ChargingUpInpFile, "ChargingUpIFile: %255s\n", ChargingUpIFile);
2567 fscanf(ChargingUpInpFile, "ChUpFactor: %lg\n", &ChUpFactor);
2568 if (1) {
2569 printf("ChargingUpEFile: %s\n", ChargingUpEFile);
2570 printf("ChargingUpIFile: %s\n", ChargingUpIFile);
2571 printf("ChUpFactor: %lg\n", ChUpFactor);
2572 }
2573
2574 {
2575 FILE *fptrChargingUpEFile = fopen(ChargingUpEFile, "r");
2576 if (fptrChargingUpEFile == NULL) {
2577 neBEMMessage(
"ChargingUpE file absent ... returning\n");
2578 return -10;
2579 }
2581 if (NbOfE <= 1) {
2582 neBEMMessage(
"Too few lines in ChargingUpE ... returning\n");
2583 fclose(fptrChargingUpEFile);
2584 return -11;
2585 }
2586 NbChUpEonEle = (
int *)malloc((
NbElements + 1) *
sizeof(int));
2587 if (NbChUpEonEle == NULL) {
2588 neBEMMessage(
"Memory allocation failed ... returning\n");
2589 fclose(fptrChargingUpEFile);
2590 return -11;
2591 }
2592 for (
int ele = 0; ele <=
NbElements; ++ele) {
2593
2594 NbChUpEonEle[ele] = 0;
2595 }
2596
2597
2598 char header[256];
2599 fgets(header, 256, fptrChargingUpEFile);
2600
2601 --NbOfE;
2602 if (debugFn) printf("No. of electrons: %d\n", NbOfE);
2603 const char *tmpEFile = "/tmp/ElectronTempFile.out";
2604 FILE *ftmpEF = fopen(tmpEFile, "w");
2605 if (ftmpEF == NULL) {
2606 printf("cannot open temporary output file ... returning ...\n");
2607 free(NbChUpEonEle);
2608 return -100;
2609 }
2610 FILE *fPtEChUpMap = fopen("PtEChUpMap.out", "w");
2611 if (fPtEChUpMap == NULL) {
2612 printf("cannot open PtEChUpMap.out file for writing ...\n");
2613 free(NbChUpEonEle);
2614 fclose(ftmpEF);
2615 return 110;
2616 }
2617
2618 char label;
2619 int vol, enb;
2620 double xlbend, ylbend, zlbend;
2621 double xend, yend, zend;
2623 ptintsct;
2624 for (int electron = 1; electron <= NbOfE; ++electron) {
2625 fscanf(fptrChargingUpEFile, "%c %d %d %lg %lg %lg %lg %lg %lg\n",
2626 &label, &vol, &enb, &xlbend, &xend, &ylbend, ¥d, &zlbend,
2627 &zend);
2628 xlbend *= 0.01;
2629 xend *= 0.01;
2630 ylbend *= 0.01;
2631 yend *= 0.01;
2632 zlbend *= 0.01;
2633 zend *= 0.01;
2637
2638
2639
2640
2641
2642
2643 if (xend < xlbend)
2644 {
2645 }
2646 double lseg = (xend - xlbend) * (xend - xlbend) +
2647 (yend - ylbend) * (yend - ylbend) +
2648 (zend - zlbend) * (zend - zlbend);
2650 double xgrd =
2651 (xend - xlbend) / lseg;
2652 double ygrd = (yend - ylbend) / lseg;
2653 double zgrd = (zend - zlbend) / lseg;
2654 if (debugFn) {
2655 printf("\nelectron: %d\n", electron);
2656 printf("xlbend: %lg, ylbend: %lg, zlbend: %lg\n", xlbend, ylbend,
2657 zlbend);
2658 printf("xend: %lg, yend: %lg, zend: %lg\n", xend, yend, zend);
2659 printf("xgrd: %lg, ygrd: %lg, zgrd: %lg\n", xgrd, ygrd, zgrd);
2660 fprintf(ftmpEF, "#e: %d, label: %c, vol: %d\n", electron, label,
2661 vol);
2662 fprintf(ftmpEF, "%lg %lg %lg\n", xlbend, ylbend, zlbend);
2663 fprintf(ftmpEF, "%lg %lg %lg\n", xend, yend, zend);
2664 fprintf(ftmpEF, "#xgrd: %lg, ygrd: %lg, zgrd: %lg\n", xgrd, ygrd,
2665 zgrd);
2666 fprintf(ftmpEF, "\n");
2667 }
2668
2669
2670
2671
2672
2673
2674
2675
2676
2677
2678
2679
2680
2681
2682 double SumOfAngles;
2683 int PrimIntsctd = -1,
2684 EleIntsctd = -1;
2685 int nearestprim = -1;
2686 double dist = 1.0e6, mindist = 1.0e6;
2687
2690 continue;
2691
2692 int intersect = 0, extrasect = 1;
2693 int InPrim = 0, InEle = 0;
2694 if (debugFn)
2695 printf("prim: %d, mindist: %lg, nearestprim: %d\n", prim,
2696 mindist, nearestprim);
2697
2698
2699
2700
2701
2702
2704 if (debugFn) {
2705 printf("prim: %d\n", prim);
2706 printf(
"vertex0: %lg, %lg, %lg\n",
XVertex[prim][0],
2708 printf(
"vertex1: %lg, %lg, %lg\n",
XVertex[prim][1],
2710 printf(
"vertex2: %lg, %lg, %lg\n",
XVertex[prim][2],
2713 printf(
"vertex3: %lg, %lg, %lg\n",
XVertex[prim][3],
2715 }
2716 fprintf(ftmpEF, "#prim: %d\n", prim);
2717 fprintf(ftmpEF,
"%lg %lg %lg\n",
XVertex[prim][0],
2719 fprintf(ftmpEF,
"%lg %lg %lg\n",
XVertex[prim][1],
2721 fprintf(ftmpEF,
"%lg %lg %lg\n",
XVertex[prim][2],
2724 fprintf(ftmpEF,
"%lg %lg %lg\n",
XVertex[prim][3],
2726 }
2727 fprintf(ftmpEF,
"%lg %lg %lg\n",
XVertex[prim][0],
2729 fprintf(ftmpEF, "\n");
2730 fflush(stdout);
2731 }
2732
2733
2734
2735 double a =
XNorm[prim];
2736 double b =
YNorm[prim];
2737 double c =
ZNorm[prim];
2740
2741
2742 dist = (xend * a + yend * b + zend * c + d) /
2743 sqrt(a * a + b * b + c * c);
2745 if (prim == 1) {
2746 mindist = dist;
2747 nearestprim = prim;
2748 }
2749 if ((prim == 1) && debugFn)
2750 printf(
2751 "after prim == 1 check mindist: %lg, nearestprim: %d\n",
2752 mindist, nearestprim);
2753
2754
2755
2756
2757
2758
2759
2760
2761 double norm1 =
sqrt(a * a + b * b + c * c);
2762 double norm2 =
sqrt(xgrd * xgrd + ygrd * ygrd + zgrd * zgrd);
2763 double denom =
2764 a * xgrd + b * ygrd + c * zgrd;
2765 double tol =
2766 1.0e-16;
2767 intersect = extrasect = 0;
2768
2769 if (debugFn) {
2770 printf("a, b, c, d, dist: %lg, %lg, %lg, %lg, %lg\n", a, b, c,
2771 d, dist);
2772 printf("vector n: ai + bj + ck\n");
2773 printf("vector v: xgrd, ygrd, zgrd: %lg, %lg, %lg\n", xgrd,
2774 ygrd, zgrd);
2775 printf("norm1, norm2, (vec n . vec v) denom: %lg, %lg, %lg\n",
2776 norm1, norm2, denom);
2777 printf("if vec n . vec v == 0, line and plane parallel\n");
2778 fflush(stdout);
2779 }
2780
2781 if (
fabs(denom) < tol * norm1 * norm2) {
2782
2783 if (
fabs(a * xlbend + b * ylbend + c * zlbend + d) <=
2784 1.0e-16) {
2785 intersect = 1;
2786 extrasect = 0;
2787 ptintsct.
X = xlbend;
2788 ptintsct.
Y = ylbend;
2789 ptintsct.
Z = zlbend;
2790 } else {
2791 intersect = 0;
2792 extrasect = 1;
2795 0.0;
2797 }
2798 if (debugFn) {
2799 printf("line and plane parallel ...\n");
2800 printf("intersect: %d, extrasect: %d\n", intersect,
2801 extrasect);
2802 printf(
"intersection point: %lg, %lg, %lg\n", ptintsct.
X,
2803 ptintsct.
Y, ptintsct.
Z);
2804 }
2805 } else {
2806 intersect = 1;
2807 double t =
2808 -(a * xlbend + b * ylbend + c * zlbend + d) / denom;
2809
2810
2811
2812
2813
2814 if ((t < 0.0) ||
2816 {
2817 extrasect = 1;
2818 ptintsct.
X = xlbend + t * xgrd;
2819 ptintsct.
Y = ylbend + t * ygrd;
2820 ptintsct.
Z = zlbend + t * zgrd;
2821 } else {
2822 extrasect = 0;
2823 ptintsct.
X = xlbend + t * xgrd;
2824 ptintsct.
Y = ylbend + t * ygrd;
2825 ptintsct.
Z = zlbend + t * zgrd;
2826 }
2827 if (debugFn) {
2828 printf("line and plane NOT parallel ...\n");
2829 printf("intersect: %d, extrasect: %d\n", intersect,
2830 extrasect);
2831 printf(
"intersection point: %lg, %lg, %lg\n", ptintsct.
X,
2832 ptintsct.
Y, ptintsct.
Z);
2833 printf("t, lseg: %lg, %lg\n", t, lseg);
2834 printf(
2835 "for an interesting intersection, lseg > t > 0.0 "
2836 "...\n\n");
2837 fflush(stdout);
2838 }
2839 }
2840 }
2841 else {
2842 dist = -1.0;
2843 intersect = 0;
2844 extrasect = 0;
2845 continue;
2846 }
2847
2848 if (dist < mindist) {
2849 mindist = dist;
2850 nearestprim = prim;
2851 }
2852 if (debugFn)
2853 printf("nearestprim: %d, mindist: %lg\n\n", nearestprim,
2854 mindist);
2855
2856
2857
2858
2859
2860 if ((intersect == 1) && (extrasect == 0)) {
2861
2862
2878 }
2879
2882 InPrim = 1;
2883 PrimIntsctd = prim;
2884 }
2885 if (debugFn) {
2886
2887 printf("Prim: %d\n", prim);
2888 printf(
"ptintsct: %lg, %lg, %lg\n", ptintsct.
X, ptintsct.
Y,
2890 printf("nvert: %d\n", nvert);
2891 printf("polynode0: %lg, %lg, %lg\n", polynode[0].X,
2892 polynode[0].Y, polynode[0].Z);
2893 printf("polynode1: %lg, %lg, %lg\n", polynode[1].X,
2894 polynode[1].Y, polynode[1].Z);
2895 printf("polynode2: %lg, %lg, %lg\n", polynode[2].X,
2896 polynode[2].Y, polynode[2].Z);
2897 if (nvert == 4) {
2898 printf("polynode3: %lg, %lg, %lg\n", polynode[3].X,
2899 polynode[3].Y, polynode[3].Z);
2900 }
2901 printf("SumOfAngles: %lg, InPrim: %d\n", SumOfAngles, InPrim);
2902 fflush(stdout);
2903 }
2904
2906 continue;
2907 }
2908
2909
2910
2911 if (InPrim) {
2912 InEle = 0;
2914 ++ele) {
2915 nvert = 0;
2916 if ((
EleArr + ele - 1)->G.Type == 3) nvert = 3;
2917 if ((
EleArr + ele - 1)->G.Type == 4) nvert = 4;
2918 if (!nvert) {
2920 "no vertex in element! ... neBEMKnownCharges ...\n");
2921 fclose(fPtEChUpMap);
2922 return -20;
2923 }
2924
2925 polynode[0].
X = (
EleArr + ele - 1)->G.Vertex[0].X;
2926 polynode[0].Y = (
EleArr + ele - 1)->G.Vertex[0].Y;
2927 polynode[0].
Z = (
EleArr + ele - 1)->G.Vertex[0].Z;
2928 polynode[1].X = (
EleArr + ele - 1)->G.Vertex[1].X;
2929 polynode[1].
Y = (
EleArr + ele - 1)->G.Vertex[1].Y;
2930 polynode[1].Z = (
EleArr + ele - 1)->G.Vertex[1].Z;
2931 polynode[2].
X = (
EleArr + ele - 1)->G.Vertex[2].X;
2932 polynode[2].Y = (
EleArr + ele - 1)->G.Vertex[2].Y;
2933 polynode[2].
Z = (
EleArr + ele - 1)->G.Vertex[2].Z;
2934 if (nvert == 4) {
2935 polynode[3].
X = (
EleArr + ele - 1)->G.Vertex[3].X;
2936 polynode[3].Y = (
EleArr + ele - 1)->G.Vertex[3].Y;
2937 polynode[3].
Z = (
EleArr + ele - 1)->G.Vertex[3].Z;
2938 }
2939
2940
2943 InEle = 1;
2944 if (debugFn) {
2945
2946 printf("Ele: %d\n", ele);
2947 printf(
"ptintsct: %lg, %lg, %lg\n", ptintsct.
X,
2948 ptintsct.
Y, ptintsct.
Z);
2949 printf("nvert: %d\n", nvert);
2950 printf("polynode0: %lg, %lg, %lg\n", polynode[0].X,
2951 polynode[0].Y, polynode[0].Z);
2952 printf("polynode1: %lg, %lg, %lg\n", polynode[1].X,
2953 polynode[1].Y, polynode[1].Z);
2954 printf("polynode2: %lg, %lg, %lg\n", polynode[2].X,
2955 polynode[2].Y, polynode[2].Z);
2956 if (nvert == 4) {
2957 printf("polynode3: %lg, %lg, %lg\n", polynode[3].X,
2958 polynode[3].Y, polynode[3].Z);
2959 }
2960 printf("SumOfAngles: %lg, InEle: %d\n", SumOfAngles,
2961 InEle);
2962 fflush(stdout);
2963 }
2964 if (InEle) {
2965 ptintsct.
X = (
EleArr + ele - 1)->G.Origin.X;
2966 ptintsct.
Y = (
EleArr + ele - 1)->G.Origin.Y;
2967 ptintsct.
Z = (
EleArr + ele - 1)->G.Origin.Z;
2968
2969 EleIntsctd = ele;
2970 NbChUpEonEle[ele]++;
2971 fprintf(fPtEChUpMap, "%d %lg %lg %lg %d %d %d %d\n",
2972 electron, ptintsct.
X, ptintsct.
Y, ptintsct.
Z,
2973 prim, InPrim, ele, InEle);
2974
2975 if (debugFn) {
2976 printf("# electron: %d\n", electron);
2977 printf("%lg %lg %lg\n", xlbend, ylbend, zlbend);
2978 printf("%lg %lg %lg\n", xend, yend, zend);
2979 printf(
"%lg, %lg, %lg\n", ptintsct.
X, ptintsct.
Y,
2981 printf("# Associated primitive: %d\n", prim);
2982 printf(
2983 "# Associated element and origin: %d, %lg, %lg, "
2984 "%lg\n",
2985 ele, (
EleArr + ele - 1)->G.Origin.X,
2986 (
EleArr + ele - 1)->G.Origin.Y,
2987 (
EleArr + ele - 1)->G.Origin.Z);
2988 printf("#NbChUpEonEle on element: %d\n",
2989 NbChUpEonEle[ele]);
2990 fprintf(ftmpEF, "#Element: %d\n", ele);
2991 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[0].X,
2992 polynode[0].Y, polynode[0].Z);
2993 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[1].X,
2994 polynode[1].Y, polynode[1].Z);
2995 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[2].X,
2996 polynode[2].Y, polynode[2].Z);
2997 if (nvert == 4) {
2998 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[3].X,
2999 polynode[3].Y, polynode[3].Z);
3000 }
3001 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[0].X,
3002 polynode[0].Y, polynode[0].Z);
3003 fprintf(ftmpEF, "\n");
3004 fflush(stdout);
3005 }
3006 break;
3007 }
3008 }
3009
3010 if (InEle)
3011 break;
3012 else {
3014 "Element not identified ... neBEMKnownCharges\n");
3015 return -2;
3016 }
3017 }
3018 }
3019
3020 if ((InPrim) && (intersect) && (!extrasect) && (InEle)) {
3021
3022
3023 break;
3024 }
3025
3026
3027
3028 if (prim ==
3030 {
3031 int nvert;
3034 double distele = 1.0e6;
3035 double mindistele = 1.0e6;
3036
3037 if (debugFn) {
3038 printf("prim == (NbPrimitives) ... checking nearest ...\n");
3039 printf("nearestprim: %d, mindist: %lg\n", nearestprim,
3040 mindist);
3041 }
3042
3043 if (mindist <= 10.0e-6) {
3044 PrimIntsctd = nearestprim;
3045 InPrim = 1;
3046 } else {
3047 InPrim = 0;
3048 InEle = 0;
3049 break;
3050 }
3051
3054 nvert = 0;
3055 if ((
EleArr + ele - 1)->G.Type == 3) nvert = 3;
3056 if ((
EleArr + ele - 1)->G.Type == 4) nvert = 4;
3057 if (!nvert) {
3059 "no vertex element! ... neBEMKnownCharges ...\n");
3060 return -20;
3061 }
3062
3063
3064
3065
3066
3067
3068
3069
3070
3071
3072
3073
3074
3075
3076
3077
3078
3079
3080
3081
3082
3083
3084
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
3102
3103
3104
3105
3106
3107
3108
3109
3110
3111
3112
3113
3114
3115
3116
3117
3118
3119
3120
3121
3122
3123
3124
3125
3126
3127
3128
3129
3130
3131
3133 eleOrigin.
X = (
EleArr + ele - 1)->G.Origin.X;
3134 eleOrigin.
Y = (
EleArr + ele - 1)->G.Origin.Y;
3135 eleOrigin.
Z = (
EleArr + ele - 1)->G.Origin.Z;
3136 distele = (eleOrigin.
X - xend) * (eleOrigin.
X - xend) +
3137 (eleOrigin.
Y - yend) * (eleOrigin.
Y - yend) +
3138 (eleOrigin.
Z - zend) * (eleOrigin.
Z - zend);
3139 distele =
pow(distele, 0.5);
3140
3142 mindistele = distele;
3143 nearestele = ele;
3144 }
3145 if (distele < mindistele) {
3146 mindistele = distele;
3147 nearestele = ele;
3148 }
3149
3150 if (debugFn) {
3151
3152
3153
3154
3155
3156 printf(
3157 "distele: %lg, mindistele: %lg,from nearest ele "
3158 "origin: %d\n",
3159 distele, mindistele, nearestele);
3160 fflush(stdout);
3161 }
3162
3163
3164 }
3165
3166
3167
3168 EleIntsctd = nearestele;
3169 InEle = 1;
3170 ptintsct.
X = (
EleArr + EleIntsctd - 1)->G.Origin.X;
3171 ptintsct.
Y = (
EleArr + EleIntsctd - 1)->G.Origin.Y;
3172 ptintsct.
Z = (
EleArr + EleIntsctd - 1)->G.Origin.Z;
3173 NbChUpEonEle[EleIntsctd]++;
3174
3175 fprintf(fPtEChUpMap, "%d %lg %lg %lg %d %d %d %d\n", electron,
3176 ptintsct.
X, ptintsct.
Y, ptintsct.
Z, PrimIntsctd, InPrim,
3177 EleIntsctd, InEle);
3178
3179
3180 if (debugFn) {
3181 printf("# electron: %d\n", electron);
3182 printf("%lg %lg %lg\n", xlbend, ylbend, zlbend);
3183 printf("%lg %lg %lg\n", xend, yend, zend);
3184 printf(
"%lg, %lg, %lg\n", ptintsct.
X, ptintsct.
Y, ptintsct.
Z);
3185 printf("# Associated primitive: %d\n", PrimIntsctd);
3186 printf("# Associated element and origin: %d, %lg, %lg, %lg\n",
3187 EleIntsctd, (
EleArr + EleIntsctd - 1)->G.Origin.X,
3188 (
EleArr + EleIntsctd - 1)->G.Origin.Y,
3189 (
EleArr + EleIntsctd - 1)->G.Origin.Z);
3190 printf("#NbChUpEonEle on element: %d\n",
3191 NbChUpEonEle[EleIntsctd]);
3192 fflush(stdout);
3193
3194 fprintf(ftmpEF, "#Element: %d\n", EleIntsctd);
3195 polynode[0].
X = (
EleArr + EleIntsctd - 1)->G.Vertex[0].X;
3196 polynode[0].Y = (
EleArr + EleIntsctd - 1)->G.Vertex[0].Y;
3197 polynode[0].
Z = (
EleArr + EleIntsctd - 1)->G.Vertex[0].Z;
3198 polynode[1].X = (
EleArr + EleIntsctd - 1)->G.Vertex[1].X;
3199 polynode[1].
Y = (
EleArr + EleIntsctd - 1)->G.Vertex[1].Y;
3200 polynode[1].Z = (
EleArr + EleIntsctd - 1)->G.Vertex[1].Z;
3201 polynode[2].
X = (
EleArr + EleIntsctd - 1)->G.Vertex[2].X;
3202 polynode[2].Y = (
EleArr + EleIntsctd - 1)->G.Vertex[2].Y;
3203 polynode[2].
Z = (
EleArr + EleIntsctd - 1)->G.Vertex[2].Z;
3204 if (nvert == 4) {
3205 polynode[3].
X = (
EleArr + EleIntsctd - 1)->G.Vertex[3].X;
3206 polynode[3].Y = (
EleArr + EleIntsctd - 1)->G.Vertex[3].Y;
3207 polynode[3].
Z = (
EleArr + EleIntsctd - 1)->G.Vertex[3].Z;
3208 }
3209 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[0].X, polynode[0].Y,
3210 polynode[0].Z);
3211 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[1].X, polynode[1].Y,
3212 polynode[1].Z);
3213 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[2].X, polynode[2].Y,
3214 polynode[2].Z);
3215 if (nvert == 4) {
3216 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[3].X,
3217 polynode[3].Y, polynode[3].Z);
3218 }
3219 fprintf(ftmpEF, "%lg %lg %lg\n", polynode[0].X, polynode[0].Y,
3220 polynode[0].Z);
3221 fprintf(ftmpEF, "\n");
3222 }
3223 }
3224
3225 }
3226
3227 if (debugFn)
3228 printf("writing file for checking electron positions\n");
3229
3230 if (debugFn)
3231
3232 {
3233 char elecposdbg[256], enbstr[10];
3234 sprintf(enbstr, "%d", electron);
3235 strcpy(elecposdbg, "/tmp/Electron");
3236 strcat(elecposdbg, enbstr);
3237 strcat(elecposdbg, ".out");
3238 FILE *fepd = fopen(elecposdbg, "w");
3239 if (fepd == NULL) {
3240 printf(
3241 "cannot open writable file to debug electron positions "
3242 "...\n");
3243 printf("returning ...\n");
3244 return -111;
3245 }
3246
3247
3248 fprintf(fepd, "#electron: %d %d\n", enb,
3249 electron);
3250 fprintf(fepd, "#last but end position:\n");
3251 fprintf(fepd, "%lg %lg %lg\n", xlbend, ylbend, zlbend);
3252 fprintf(fepd, "#end position:\n");
3253 fprintf(fepd, "%lg %lg %lg\n\n", xend, yend, zend);
3254 fprintf(fepd, "#intersected primitive number: %d\n", PrimIntsctd);
3255 if (PrimIntsctd >= 1) {
3256 fprintf(fepd,
"#PrimType: %d\n",
PrimType[PrimIntsctd]);
3257 fprintf(fepd, "#prim vertices:\n");
3258 fprintf(fepd,
"%lg %lg %lg\n",
XVertex[PrimIntsctd][0],
3260 fprintf(fepd,
"%lg %lg %lg\n",
XVertex[PrimIntsctd][1],
3262 fprintf(fepd,
"%lg %lg %lg\n",
XVertex[PrimIntsctd][2],
3265 fprintf(fepd,
"%lg %lg %lg\n",
XVertex[PrimIntsctd][3],
3267 }
3268 fprintf(fepd,
"%lg %lg %lg\n",
XVertex[PrimIntsctd][0],
3270 fprintf(fepd, "\n");
3271
3272 fprintf(fepd, "#ptintsct:\n");
3273 fprintf(fepd,
"%lg %lg %lg\n", ptintsct.
X, ptintsct.
Y,
3275 fprintf(fepd, "\n");
3276 }
3277 fprintf(fepd, "#intersected element number: %d\n", EleIntsctd);
3278 if (EleIntsctd >= 1) {
3279 int gtype = (
EleArr + EleIntsctd - 1)->G.Type;
3280 fprintf(fepd, "#EleType: %d\n", gtype);
3281 fprintf(fepd, "#element vertices:\n");
3282 double x0 = (
EleArr + EleIntsctd - 1)->G.Vertex[0].X;
3283 double y0 = (
EleArr + EleIntsctd - 1)->G.Vertex[0].Y;
3284 double z0 = (
EleArr + EleIntsctd - 1)->G.Vertex[0].Z;
3285 double x1 = (
EleArr + EleIntsctd - 1)->G.Vertex[1].X;
3286 double y1 = (
EleArr + EleIntsctd - 1)->G.Vertex[1].Y;
3287 double z1 = (
EleArr + EleIntsctd - 1)->G.Vertex[1].Z;
3288 double x2 = (
EleArr + EleIntsctd - 1)->G.Vertex[2].X;
3289 double y2 = (
EleArr + EleIntsctd - 1)->G.Vertex[2].Y;
3290 double z2 = (
EleArr + EleIntsctd - 1)->G.Vertex[2].Z;
3291 fprintf(fepd, "%lg %lg %lg\n", x0, y0, z0);
3292 fprintf(fepd, "%lg %lg %lg\n", x1, y1, z1);
3293 fprintf(fepd, "%lg %lg %lg\n", x2, y2, z2);
3294 if (gtype == 4) {
3295 double x3 = (
EleArr + EleIntsctd - 1)->G.Vertex[3].X;
3296 double y3 = (
EleArr + EleIntsctd - 1)->G.Vertex[3].Y;
3297 double z3 = (
EleArr + EleIntsctd - 1)->G.Vertex[3].Z;
3298 fprintf(fepd, "%lg %lg %lg\n", x3, y3, z3);
3299 }
3300 fprintf(fepd, "%lg %lg %lg\n", x0, y0, z0);
3301 fprintf(fepd, "\n");
3302
3303 fprintf(fepd, "#ptintsct:\n");
3304 fprintf(fepd,
"%lg %lg %lg\n", ptintsct.
X, ptintsct.
Y,
3306 fprintf(fepd, "\n");
3307 }
3308
3309 fclose(fepd);
3310 }
3311 if (debugFn)
3312 printf("done writing file for checking electron positions\n");
3313 }
3314 fclose(fPtEChUpMap);
3315 if (debugFn) printf("done for all electrons\n");
3316
3317 FILE *fEleEChUpMap = fopen("EleEChUpMap.out", "w");
3318 if (fEleEChUpMap == NULL) {
3319 printf("cannot open EleEChUpMap.out file for writing ...\n");
3320 return 111;
3321 }
3322 for (
int ele = 1; ele <=
NbElements; ++ele) {
3323 (
EleArr + ele - 1)->Assigned +=
3324 ChUpFactor *
Q_E * NbChUpEonEle[ele] / (
EleArr + ele - 1)->G.dA;
3325 fprintf(fEleEChUpMap, "%d %lg %lg %lg %d %lg\n", ele,
3326 (
EleArr + ele - 1)->G.Origin.X,
3327 (
EleArr + ele - 1)->G.Origin.Y,
3328 (
EleArr + ele - 1)->G.Origin.Z, NbChUpEonEle[ele],
3329 (
EleArr + ele - 1)->Assigned);
3330 }
3331 fclose(fEleEChUpMap);
3332
3333 fclose(ftmpEF);
3334 free(NbChUpEonEle);
3335 }
3336
3337 {
3338 FILE *fptrChargingUpIFile = fopen(ChargingUpIFile, "r");
3339 if (fptrChargingUpIFile == NULL) {
3340 neBEMMessage(
"ChargingUpI file absent ... returning\n");
3341 return -10;
3342 }
3344 if (NbOfI <= 1) {
3345 neBEMMessage(
"Too few lines in ChargingUpI ... returning\n");
3346 fclose(fptrChargingUpIFile);
3347 return -11;
3348 }
3349
3350 NbChUpIonEle = (
int *)malloc((
NbElements + 1) *
sizeof(int));
3351 if (NbChUpIonEle == NULL) {
3352 neBEMMessage(
"Memory allocation failed ... returning\n");
3353 fclose(fptrChargingUpIFile);
3354 return -11;
3355 }
3356 for (
int ele = 0; ele <=
NbElements; ++ele) {
3357
3358 NbChUpIonEle[ele] = 0;
3359 }
3360
3361
3362 char header[256];
3363 fgets(header, 256, fptrChargingUpIFile);
3364
3365 --NbOfI;
3366 if (debugFn) printf("No. of ions: %d\n", NbOfI);
3367 const char *tmpIFile = "/tmp/IonTempFile.out";
3368 FILE *ftmpIF = fopen(tmpIFile, "w");
3369 if (ftmpIF == NULL) {
3370 printf("cannot open temporary ion output file ... returning ...\n");
3371 free(NbChUpEonEle);
3372 return -100;
3373 }
3374 FILE *fPtIChUpMap = fopen("PtIChUpMap.out", "w");
3375 if (fPtIChUpMap == NULL) {
3376 printf("cannot open PtIChUpMap.out file for writing ...\n");
3377 fclose(ftmpIF);
3378 free(NbChUpEonEle);
3379 return 110;
3380 }
3381
3382 char label;
3383 int inb, vol;
3384 double xlbend, ylbend, zlbend;
3385 double xend, yend, zend;
3387 for (int ion = 1; ion <= NbOfI; ++ion) {
3388 fscanf(fptrChargingUpIFile, "%c %d %d %lg %lg %lg %lg %lg %lg\n",
3389 &label, &vol, &inb, &xlbend, &xend, &ylbend, ¥d, &zlbend,
3390 &zend);
3391 xlbend *= 0.01;
3392 xend *= 0.01;
3393 ylbend *= 0.01;
3394 yend *= 0.01;
3395 zlbend *= 0.01;
3396 zend *= 0.01;
3400
3401
3402
3403
3404
3405
3406 if (xend < xlbend)
3407 {
3408 }
3409 double lseg = (xend - xlbend) * (xend - xlbend) +
3410 (yend - ylbend) * (yend - ylbend) +
3411 (zend - zlbend) * (zend - zlbend);
3413 double xgrd =
3414 (xend - xlbend) / lseg;
3415 double ygrd = (yend - ylbend) / lseg;
3416 double zgrd = (zend - zlbend) / lseg;
3417 if (debugFn) {
3418 printf("\nion: %d\n", ion);
3419 printf("xlbend: %lg, ylbend: %lg, zlbend: %lg\n", xlbend, ylbend,
3420 zlbend);
3421 printf("xend: %lg, yend: %lg, zend: %lg\n", xend, yend, zend);
3422 printf("xgrd: %lg, ygrd: %lg, zgrd: %lg\n", xgrd, ygrd, zgrd);
3423 fprintf(ftmpIF, "#e: %d, label: %c, vol: %d\n", ion, label, vol);
3424 fprintf(ftmpIF, "%lg %lg %lg\n", xlbend, ylbend, zlbend);
3425 fprintf(ftmpIF, "%lg %lg %lg\n", xend, yend, zend);
3426 fprintf(ftmpIF, "#xgrd: %lg, ygrd: %lg, zgrd: %lg\n", xgrd, ygrd,
3427 zgrd);
3428 fprintf(ftmpIF, "\n");
3429 }
3430
3431
3432
3433
3434
3435
3436
3437
3438
3439
3440
3441
3442
3443
3444 int PrimIntsctd = -1,
3445 EleIntsctd = -1;
3446 int nearestprim = -1;
3447 double dist = 1.0e6, mindist = 1.0e6;
3448 double SumOfAngles;
3449
3452 continue;
3453
3454 int intersect = 0, extrasect = 1;
3455 int InPrim = 0, InEle = 0;
3456 if (debugFn)
3457 printf("prim: %d, mindist: %lg, nearestprim: %d\n", prim,
3458 mindist, nearestprim);
3459
3460
3461
3462
3463
3464
3465
3466
3468 if (debugFn) {
3469 printf("prim: %d\n", prim);
3470 printf(
"vertex0: %lg, %lg, %lg\n",
XVertex[prim][0],
3472 printf(
"vertex1: %lg, %lg, %lg\n",
XVertex[prim][1],
3474 printf(
"vertex2: %lg, %lg, %lg\n",
XVertex[prim][2],
3477 printf(
"vertex3: %lg, %lg, %lg\n",
XVertex[prim][3],
3479 }
3480 fprintf(ftmpIF, "#prim: %d\n", prim);
3481 fprintf(ftmpIF,
"%lg %lg %lg\n",
XVertex[prim][0],
3483 fprintf(ftmpIF,
"%lg %lg %lg\n",
XVertex[prim][1],
3485 fprintf(ftmpIF,
"%lg %lg %lg\n",
XVertex[prim][2],
3488 fprintf(ftmpIF,
"%lg %lg %lg\n",
XVertex[prim][3],
3490 }
3491 fprintf(ftmpIF,
"%lg %lg %lg\n",
XVertex[prim][0],
3493 fprintf(ftmpIF, "\n");
3494 fflush(stdout);
3495 }
3496
3497
3498
3502
3503
3504 dist =
3506 zend *
ZNorm[prim] + d) /
3510 if (prim == 1) {
3511 mindist = dist;
3512 nearestprim = prim;
3513 }
3514 if ((prim == 1) && debugFn)
3515 printf(
3516 "after prim == 1 check mindist: %lg, nearestprim: %d\n",
3517 mindist, nearestprim);
3518
3519
3520
3521
3522
3523
3524
3525
3526 double a =
XNorm[prim];
3527 double b =
YNorm[prim];
3528 double c =
ZNorm[prim];
3529 double norm1 =
sqrt(a * a + b * b + c * c);
3530 double norm2 =
sqrt(xgrd * xgrd + ygrd * ygrd + zgrd * zgrd);
3531 double denom =
3532 a * xgrd + b * ygrd + c * zgrd;
3533 double tol =
3534 1.0e-12;
3535 intersect = extrasect = 0;
3536
3537 if (debugFn) {
3538 printf("a, b, c, d, dist: %lg, %lg, %lg, %lg, %lg\n", a, b, c,
3539 d, dist);
3540 printf("vector n: ai + bj + ck\n");
3541 printf("vector v: xgrd, ygrd, zgrd: %lg, %lg, %lg\n", xgrd,
3542 ygrd, zgrd);
3543 printf("norm1, norm2, (vec n . vec v) denom: %lg, %lg, %lg\n",
3544 norm1, norm2, denom);
3545 printf("if vec n . vec v == 0, line and plane parallel\n");
3546 fflush(stdout);
3547 }
3548
3549 if (
fabs(denom) < tol * norm1 * norm2) {
3550
3551 if (a * xlbend + b * ylbend + c * zlbend + d ==
3552 0.0)
3553 {
3554 intersect = 1;
3555 extrasect = 0;
3556 ptintsct.
X = xlbend;
3557 ptintsct.
Y = ylbend;
3558 ptintsct.
Z = zlbend;
3559 } else {
3560 intersect = 0;
3561 extrasect = 0;
3564 0.0;
3566 }
3567 if (debugFn) {
3568 printf("line and plane parallel ...\n");
3569 printf("intersect: %d, extrasect: %d\n", intersect,
3570 extrasect);
3571 printf(
"intersection point: %lg, %lg, %lg\n", ptintsct.
X,
3572 ptintsct.
Y, ptintsct.
Z);
3573 }
3574 } else {
3575 intersect = 1;
3576 double t =
3577 -(a * xlbend + b * ylbend + c * zlbend + d) / denom;
3578
3579
3580
3581
3582
3583 if ((t < 0.0) || (
fabs(t) >
fabs(lseg))) {
3584 extrasect = 1;
3585 ptintsct.
X = xlbend + t * xgrd;
3586 ptintsct.
Y = ylbend + t * ygrd;
3587 ptintsct.
Z = zlbend + t * zgrd;
3588 } else {
3589 extrasect = 0;
3590 ptintsct.
X = xlbend + t * xgrd;
3591 ptintsct.
Y = ylbend + t * ygrd;
3592 ptintsct.
Z = zlbend + t * zgrd;
3593 }
3594 if (debugFn) {
3595 printf("line and plane NOT parallel ...\n");
3596 printf("intersect: %d, extrasect: %d\n", intersect,
3597 extrasect);
3598 printf(
"intersection point: %lg, %lg, %lg\n", ptintsct.
X,
3599 ptintsct.
Y, ptintsct.
Z);
3600 printf("t, lseg: %lg, %lg\n", t, lseg);
3601 printf(
3602 "for an interesting intersection, lseg > t > 0.0 "
3603 "...\n\n");
3604 fflush(stdout);
3605 }
3606 }
3607 }
3608 else {
3609 dist = -1.0;
3610 intersect = 0;
3611 extrasect = 0;
3612 }
3613
3614 if (dist < mindist) {
3615 mindist = dist;
3616 nearestprim = prim;
3617 }
3618
3619
3620
3621
3622
3623 if ((intersect == 1) && (extrasect == 0)) {
3624
3625
3641 }
3642
3645 InPrim = 1;
3646 PrimIntsctd = prim;
3647 }
3648 if (debugFn) {
3649
3650 printf("Prim: %d\n", prim);
3651 printf(
"ptintsct: %lg, %lg, %lg\n", ptintsct.
X, ptintsct.
Y,
3653 printf("nvert: %d\n", nvert);
3654 printf("polynode0: %lg, %lg, %lg\n", polynode[0].X,
3655 polynode[0].Y, polynode[0].Z);
3656 printf("polynode1: %lg, %lg, %lg\n", polynode[1].X,
3657 polynode[1].Y, polynode[1].Z);
3658 printf("polynode2: %lg, %lg, %lg\n", polynode[2].X,
3659 polynode[2].Y, polynode[2].Z);
3660 if (nvert == 4) {
3661 printf("polynode3: %lg, %lg, %lg\n", polynode[3].X,
3662 polynode[3].Y, polynode[3].Z);
3663 }
3664 printf("SumOfAngles: %lg, InPrim: %d\n", SumOfAngles, InPrim);
3665 fflush(stdout);
3666 }
3667 if (!InPrim) continue;
3668
3669
3670
3671 InEle = 0;
3673 ++ele) {
3674 nvert = 0;
3675 if ((
EleArr + ele - 1)->G.Type == 3) nvert = 3;
3676 if ((
EleArr + ele - 1)->G.Type == 4) nvert = 4;
3677 if (!nvert) {
3679 "no vertex in element! ... neBEMKnownCharges ...\n");
3680 if (fPtIChUpMap) fclose(fPtIChUpMap);
3681 return -20;
3682 }
3683
3684 polynode[0].
X = (
EleArr + ele - 1)->G.Vertex[0].X;
3685 polynode[0].Y = (
EleArr + ele - 1)->G.Vertex[0].Y;
3686 polynode[0].
Z = (
EleArr + ele - 1)->G.Vertex[0].Z;
3687 polynode[1].X = (
EleArr + ele - 1)->G.Vertex[1].X;
3688 polynode[1].
Y = (
EleArr + ele - 1)->G.Vertex[1].Y;
3689 polynode[1].Z = (
EleArr + ele - 1)->G.Vertex[1].Z;
3690 polynode[2].
X = (
EleArr + ele - 1)->G.Vertex[2].X;
3691 polynode[2].Y = (
EleArr + ele - 1)->G.Vertex[2].Y;
3692 polynode[2].
Z = (
EleArr + ele - 1)->G.Vertex[2].Z;
3693 if (nvert == 4) {
3694 polynode[3].
X = (
EleArr + ele - 1)->G.Vertex[3].X;
3695 polynode[3].Y = (
EleArr + ele - 1)->G.Vertex[3].Y;
3696 polynode[3].
Z = (
EleArr + ele - 1)->G.Vertex[3].Z;
3697 }
3698
3699
3702 if (debugFn) {
3703
3704 printf("Ele: %d\n", ele);
3705 printf(
"ptintsct: %lg, %lg, %lg\n", ptintsct.
X, ptintsct.
Y,
3707 printf("nvert: %d\n", nvert);
3708 printf("polynode0: %lg, %lg, %lg\n", polynode[0].X,
3709 polynode[0].Y, polynode[0].Z);
3710 printf("polynode1: %lg, %lg, %lg\n", polynode[1].X,
3711 polynode[1].Y, polynode[1].Z);
3712 printf("polynode2: %lg, %lg, %lg\n", polynode[2].X,
3713 polynode[2].Y, polynode[2].Z);
3714 if (nvert == 4) {
3715 printf("polynode3: %lg, %lg, %lg\n", polynode[3].X,
3716 polynode[3].Y, polynode[3].Z);
3717 }
3718 printf("SumOfAngles: %lg, InEle: %d\n", SumOfAngles, InEle);
3719 fflush(stdout);
3720 }
3721 if (InEle) {
3722 ptintsct.
X = (
EleArr + ele - 1)->G.Origin.X;
3723 ptintsct.
Y = (
EleArr + ele - 1)->G.Origin.Y;
3724 ptintsct.
Z = (
EleArr + ele - 1)->G.Origin.Z;
3725 EleIntsctd = ele;
3726
3727 NbChUpIonEle[ele]++;
3728 fprintf(fPtIChUpMap, "%d %lg %lg %lg %d %d %d %d\n", ion,
3729 ptintsct.
X, ptintsct.
Y, ptintsct.
Z, prim, InPrim,
3730 ele, InEle);
3731
3732 if (debugFn) {
3733 printf("# ion: %d\n", ion);
3734 printf("%lg %lg %lg\n", xlbend, ylbend, zlbend);
3735 printf("%lg %lg %lg\n", xend, yend, zend);
3736 printf(
"%lg, %lg, %lg\n", ptintsct.
X, ptintsct.
Y,
3738 printf("# Associated primitive: %d\n", prim);
3739 printf(
3740 "# Associated element and origin: %d, %lg, %lg, "
3741 "%lg\n",
3742 ele, (
EleArr + ele - 1)->G.Origin.X,
3743 (
EleArr + ele - 1)->G.Origin.Y,
3744 (
EleArr + ele - 1)->G.Origin.Z);
3745 printf("#NbChUpIonEle on element: %d\n",
3746 NbChUpIonEle[ele]);
3747 fprintf(ftmpIF, "#Element: %d\n", ele);
3748 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[0].X,
3749 polynode[0].Y, polynode[0].Z);
3750 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[1].X,
3751 polynode[1].Y, polynode[1].Z);
3752 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[2].X,
3753 polynode[2].Y, polynode[2].Z);
3754 if (nvert == 4) {
3755 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[3].X,
3756 polynode[3].Y, polynode[3].Z);
3757 }
3758 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[0].X,
3759 polynode[0].Y, polynode[0].Z);
3760 fprintf(ftmpIF, "\n");
3761 fflush(stdout);
3762 }
3763 break;
3764 }
3765 }
3766
3767 if (InEle) break;
3769 "Element cannot be identified ... neBEMKnownCharges\n");
3770 return -2;
3771 }
3772
3773 if ((InPrim) && (intersect) && (!extrasect) && (InEle)) {
3774
3775
3776 break;
3777 }
3778
3779
3780
3782
3783 int nvert;
3786 double distele = 1.0e6;
3787 double mindistele = 1.0e6;
3788
3789 if (debugFn) {
3790 printf("prim == (NbPrimitives) ... checking nearest ...\n");
3791 printf("nearestprim: %d, mindist: %lg\n", nearestprim,
3792 mindist);
3793 }
3794
3795 if (mindist <= 10.0e-6) {
3796 PrimIntsctd = nearestprim;
3797 InPrim = 1;
3798 } else {
3799 InPrim = 0;
3800 InEle = 0;
3801 break;
3802 }
3803
3806 nvert = 0;
3807 if ((
EleArr + ele - 1)->G.Type == 3) nvert = 3;
3808 if ((
EleArr + ele - 1)->G.Type == 4) nvert = 4;
3809 if (!nvert) {
3811 "no vertex element! ... neBEMKnownCharges ...\n");
3812 return -20;
3813 }
3814
3815
3816
3817
3818
3819
3820
3821
3822
3823
3824
3825
3826
3827
3828
3829
3830
3831
3832
3833
3834
3835
3836
3837
3838
3839
3840
3841
3842
3843
3844
3845
3846
3847
3848
3849
3850
3851
3852
3853
3854
3855
3856
3857
3858
3859
3860
3861
3862
3863
3864
3865
3866
3867
3868
3869
3870
3871
3872
3873
3874
3875
3876
3877
3878
3879
3880
3881
3882
3883
3884
3885
3887 eleOrigin.
X = (
EleArr + ele - 1)->G.Origin.X;
3888 eleOrigin.
Y = (
EleArr + ele - 1)->G.Origin.Y;
3889 eleOrigin.
Z = (
EleArr + ele - 1)->G.Origin.Z;
3890 distele = (eleOrigin.
X - xend) * (eleOrigin.
X - xend) +
3891 (eleOrigin.
Y - yend) * (eleOrigin.
Y - yend) +
3892 (eleOrigin.
Z - zend) * (eleOrigin.
Z - zend);
3893 distele =
pow(distele, 0.5);
3894
3896 mindistele = distele;
3897 nearestele = ele;
3898 }
3899 if (distele < mindistele) {
3900 mindistele = distele;
3901 nearestele = ele;
3902 }
3903
3904 if (debugFn) {
3905
3906
3907
3908
3909
3910 printf(
3911 "distele: %lg, mindist: %lg, from nearest ele: %d\n",
3912 distele, mindistele, nearestele);
3913 fflush(stdout);
3914 }
3915
3916
3917 }
3918
3919
3920
3921 EleIntsctd = nearestele;
3922 InEle = 1;
3923 ptintsct.
X = (
EleArr + EleIntsctd - 1)->G.Origin.X;
3924 ptintsct.
Y = (
EleArr + EleIntsctd - 1)->G.Origin.Y;
3925 ptintsct.
Z = (
EleArr + EleIntsctd - 1)->G.Origin.Z;
3926 NbChUpIonEle[EleIntsctd]++;
3927
3928 fprintf(fPtIChUpMap, "%d %lg %lg %lg %d %d %d %d\n", ion,
3929 ptintsct.
X, ptintsct.
Y, ptintsct.
Z, PrimIntsctd, InPrim,
3930 EleIntsctd, InEle);
3931
3932
3933 if (debugFn) {
3934 printf("# ion: %d\n", ion);
3935 printf("%lg %lg %lg\n", xlbend, ylbend, zlbend);
3936 printf("%lg %lg %lg\n", xend, yend, zend);
3937 printf(
"%lg, %lg, %lg\n", ptintsct.
X, ptintsct.
Y, ptintsct.
Z);
3938 printf("# Associated primitive: %d\n", PrimIntsctd);
3939 printf("# Associated element and origin: %d, %lg, %lg, %lg\n",
3940 EleIntsctd, (
EleArr + EleIntsctd - 1)->G.Origin.X,
3941 (
EleArr + EleIntsctd - 1)->G.Origin.Y,
3942 (
EleArr + EleIntsctd - 1)->G.Origin.Z);
3943 printf("#NbChUpIonEle on element: %d\n",
3944 NbChUpIonEle[EleIntsctd]);
3945 fprintf(ftmpIF, "#Element: %d\n", EleIntsctd);
3946 polynode[0].
X = (
EleArr + EleIntsctd - 1)->G.Vertex[0].X;
3947 polynode[0].Y = (
EleArr + EleIntsctd - 1)->G.Vertex[0].Y;
3948 polynode[0].
Z = (
EleArr + EleIntsctd - 1)->G.Vertex[0].Z;
3949 polynode[1].X = (
EleArr + EleIntsctd - 1)->G.Vertex[1].X;
3950 polynode[1].
Y = (
EleArr + EleIntsctd - 1)->G.Vertex[1].Y;
3951 polynode[1].Z = (
EleArr + EleIntsctd - 1)->G.Vertex[1].Z;
3952 polynode[2].
X = (
EleArr + EleIntsctd - 1)->G.Vertex[2].X;
3953 polynode[2].Y = (
EleArr + EleIntsctd - 1)->G.Vertex[2].Y;
3954 polynode[2].
Z = (
EleArr + EleIntsctd - 1)->G.Vertex[2].Z;
3955 if (nvert == 4) {
3956 polynode[3].
X = (
EleArr + EleIntsctd - 1)->G.Vertex[3].X;
3957 polynode[3].Y = (
EleArr + EleIntsctd - 1)->G.Vertex[3].Y;
3958 polynode[3].
Z = (
EleArr + EleIntsctd - 1)->G.Vertex[3].Z;
3959 }
3960 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[0].X, polynode[0].Y,
3961 polynode[0].Z);
3962 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[1].X, polynode[1].Y,
3963 polynode[1].Z);
3964 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[2].X, polynode[2].Y,
3965 polynode[2].Z);
3966 if (nvert == 4) {
3967 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[3].X,
3968 polynode[3].Y, polynode[3].Z);
3969 }
3970 fprintf(ftmpIF, "%lg %lg %lg\n", polynode[0].X, polynode[0].Y,
3971 polynode[0].Z);
3972 fprintf(ftmpIF, "\n");
3973 fflush(stdout);
3974 }
3975 }
3976
3977 }
3978
3979 if (debugFn)
3980 {
3981 char ionposdbg[256], inbstr[10];
3982 sprintf(inbstr, "%d", ion);
3983 strcpy(ionposdbg, "/tmp/Ion");
3984 strcat(ionposdbg, inbstr);
3985 strcat(ionposdbg, ".out");
3986 FILE *fipd = fopen(ionposdbg, "w");
3987 if (fipd == NULL) {
3988 printf(
3989 "cannot open writable file to debug ion positions ...\n");
3990 printf("returning ...\n");
3991 return -111;
3992 }
3993
3994
3995 fprintf(fipd, "#ion: %d %d\n", inb, ion);
3996 fprintf(fipd, "#last but end position:\n");
3997 fprintf(fipd, "%lg %lg %lg\n", xlbend, ylbend, zlbend);
3998 fprintf(fipd, "#end position:\n");
3999 fprintf(fipd, "%lg %lg %lg\n\n", xend, yend, zend);
4000
4001 fprintf(fipd, "#intersected primitive number: %d\n", PrimIntsctd);
4002 if (PrimIntsctd >= 1) {
4003 fprintf(fipd,
"#PrimType: %d\n",
PrimType[PrimIntsctd]);
4004 fprintf(fipd, "#prim vertices:\n");
4005 fprintf(fipd,
"%lg %lg %lg\n",
XVertex[PrimIntsctd][0],
4007 fprintf(fipd,
"%lg %lg %lg\n",
XVertex[PrimIntsctd][1],
4009 fprintf(fipd,
"%lg %lg %lg\n",
XVertex[PrimIntsctd][2],
4012 fprintf(fipd,
"%lg %lg %lg\n",
XVertex[PrimIntsctd][3],
4014 }
4015 fprintf(fipd,
"%lg %lg %lg\n",
XVertex[PrimIntsctd][0],
4017 fprintf(fipd, "\n");
4018
4019 fprintf(fipd, "#ptintsct:\n");
4020 fprintf(fipd,
"%lg %lg %lg\n", ptintsct.
X, ptintsct.
Y,
4022 fprintf(fipd, "\n");
4023 }
4024
4025 fprintf(fipd, "#intersected element number: %d\n", EleIntsctd);
4026 if (EleIntsctd >= 1) {
4027 int gtype = (
EleArr + EleIntsctd - 1)->G.Type;
4028 fprintf(fipd, "#EleType: %d\n", gtype);
4029 fprintf(fipd, "#element vertices:\n");
4030 double x0 = (
EleArr + EleIntsctd - 1)->G.Vertex[0].X;
4031 double y0 = (
EleArr + EleIntsctd - 1)->G.Vertex[0].Y;
4032 double z0 = (
EleArr + EleIntsctd - 1)->G.Vertex[0].Z;
4033 double x1 = (
EleArr + EleIntsctd - 1)->G.Vertex[1].X;
4034 double y1 = (
EleArr + EleIntsctd - 1)->G.Vertex[1].Y;
4035 double z1 = (
EleArr + EleIntsctd - 1)->G.Vertex[1].Z;
4036 double x2 = (
EleArr + EleIntsctd - 1)->G.Vertex[2].X;
4037 double y2 = (
EleArr + EleIntsctd - 1)->G.Vertex[2].Y;
4038 double z2 = (
EleArr + EleIntsctd - 1)->G.Vertex[2].Z;
4039 fprintf(fipd, "%lg %lg %lg\n", x0, y0, z0);
4040 fprintf(fipd, "%lg %lg %lg\n", x1, y1, z1);
4041 fprintf(fipd, "%lg %lg %lg\n", x2, y2, z2);
4042 if (gtype == 4) {
4043 double x3 = (
EleArr + EleIntsctd - 1)->G.Vertex[3].X;
4044 double y3 = (
EleArr + EleIntsctd - 1)->G.Vertex[3].Y;
4045 double z3 = (
EleArr + EleIntsctd - 1)->G.Vertex[3].Z;
4046 fprintf(fipd, "%lg %lg %lg\n", x3, y3, z3);
4047 }
4048 fprintf(fipd, "%lg %lg %lg\n", x0, y0, z0);
4049 fprintf(fipd, "\n");
4050
4051 fprintf(fipd, "#ptintsct:\n");
4052 fprintf(fipd,
"%lg %lg %lg\n", ptintsct.
X, ptintsct.
Y,
4054 fprintf(fipd, "\n");
4055 }
4056 fclose(fipd);
4057 }
4058 }
4059 fclose(fPtIChUpMap);
4060
4061
4062
4063 FILE *fEleEIChUpMap = fopen("EleE+IChUpMap.out", "w");
4064 if (fEleEIChUpMap == NULL) {
4065 printf("cannot open EleE+IChUpMap.out file for writing ...\n");
4066 return 111;
4067 }
4068 for (
int ele = 1; ele <=
NbElements; ++ele) {
4069 (
EleArr + ele - 1)->Assigned +=
4070 ChUpFactor *
Q_I * NbChUpIonEle[ele] / (
EleArr + ele - 1)->G.dA;
4071 fprintf(fEleEIChUpMap, "%d %lg %lg %lg %d %lg\n", ele,
4072 (
EleArr + ele - 1)->G.Origin.X,
4073 (
EleArr + ele - 1)->G.Origin.Y,
4074 (
EleArr + ele - 1)->G.Origin.Z, NbChUpIonEle[ele],
4075 (
EleArr + ele - 1)->Assigned);
4076 }
4077 fclose(fEleEIChUpMap);
4078
4079 fclose(ftmpIF);
4080 free(NbChUpIonEle);
4081 }
4082
4083 }
4084
4085 fclose(ChargingUpInpFile);
4086 }
4087
4088 if (debugFn) {
4089
4090 }
4091 }
4092
4093 return (0);
4094}
DoubleAc pow(const DoubleAc &f, double p)
DoubleAc fabs(const DoubleAc &f)
DoubleAc sqrt(const DoubleAc &f)
int neBEMGetNbOfLines(const char fname[])
double neBEMChkInPoly(int n, Point3D *p, Point3D q)
neBEMGLOBAL Element * EleArr
neBEMGLOBAL double * ZNorm
neBEMGLOBAL double ** ZVertex
neBEMGLOBAL int OptChargingUp
neBEMGLOBAL double * XNorm
neBEMGLOBAL int NbPrimitives
neBEMGLOBAL int NbElements
neBEMGLOBAL double ** YVertex
neBEMGLOBAL double * YNorm
neBEMGLOBAL double ** XVertex
neBEMGLOBAL int * ElementEnd
neBEMGLOBAL int * InterfaceType
neBEMGLOBAL int * PrimType
neBEMGLOBAL int * ElementBgn